GEE:批量处理和下载SoilGrids 250m v2.0
01 说明
1.1 任务和要求
任务:利用GEE
分别输出soc
(有机碳含量)和nitrogen
(全碳含量)在0-5cm、5-15cm、15-30cm的单波段影像。具体要求包括:
全球范围,地理坐标系WGS84,输出分辨率0.05°,应用缩放(基于比例因子scaler_factor
)还原单位,无效值处理;
1.2 参考
比例因子scaler_factor参考官网(https://gee-community-catalog.org/projects/isric/
),具体如下:
参考官方代码示例:
https://git.wur.nl/isric/soilgrids/soilgrids.notebooks/-/blob/master/markdown/access_on_gee.md
https://code.earthengine.google.com/08f74bc377b253523cd4cffe4268f66d
ps
:本文所下载的SoilGrids 250产品与GEE公开数据集检索的SoilGrids产品(ee.ImageCollection("ISRIC/SoilGrids250m/v2_0")
, 见下图)不同,GEE
收录的是旧版本(但也是ISRIC生产的但不是最新的,是静态不更新的),而本文所下载的SoilGrids产品是ISRIC托管在GEE
上的数据集是由ISRIC自己直接维护和更新的,是最新的。
02 代码说明
完整代码如下:
/*
转换因子scaler_factor见: https://gee-community-catalog.org/projects/isric/#citation
soc、nitrogen的转换因子分别为10、100,转换后单位分别为dg/kg, cg/kg
*/// 准备
var out_res = 0.05;
var soc_bands_name = ['soc_0-5cm_mean', 'soc_5-15cm_mean', 'soc_15-30cm_mean'];
var nitrogen_bands_name = ['nitrogen_0-5cm_mean', 'nitrogen_0-5cm_mean', 'nitrogen_0-5cm_mean'];
var geo_transform = [out_res, 0, -180, 0, -out_res, 90];
var no_data_value = -9999; // 无效值设置// 获取soc和nitrogen
var soc = ee.Image("projects/soilgrids-isric/soc_mean"); // 有机碳含量
var nitrogen = ee.Image("projects/soilgrids-isric/nitrogen_mean"); // 全碳含量
print(soc)
print(nitrogen)// 选取0-30cm的波段
soc = soc.select(soc_bands_name);
nitrogen = nitrogen.select(nitrogen_bands_name)// 无效值处理(等于0的掩膜掉)
soc = soc.updateMask(soc.neq(0));
nitrogen = soc.updateMask(nitrogen.neq(0));// 缩放
soc = soc.divide(10); // 缩放后单位为dg/kg
nitrogen = nitrogen.divide(100); // 缩放后单位为cg/kg// 定义地理范围
var global_region = ee.Geometry.Rectangle([-180, -90, 180, 90], 'EPSG:4326', false);// 定义导出函数
function export_band(img_band, band_name) {var file_name = 'isric_' + band_name;Export.image.toDrive({image: img_band,description: file_name,folder: 'Global_ISRIC',fileNamePrefix: file_name,region: global_region,crs: 'EPSG:4326',crsTransform: geo_transform,maxPixels: 1e13,formatOptions: {'noData': no_data_value}
});
}// 分别导出soc各个波段
soc_bands_name.forEach(function(band_name){var cur_soc = soc.select(band_name);export_band(cur_soc, band_name);
})
// 分别导出nitrogen各个波段
nitrogen_bands_name.forEach(function(band_name){var cur_nitrogen = nitrogen.select(band_name);export_band(cur_nitrogen, band_name);
})
这里有一些细节进行说明:
- 由于要求里需要输出地理坐标系WGS84且输出分辨率为0.05°。如果按照一般的输出使用
sclae
参数,则只能输出为以WGS84椭球体为基准的投影坐标系,单位为m.这是指定scale
参数之后gee会依据scale
参数自动计算crsTransform
参数(仿射系数:[x_res,0,left_right_corner_x,0, -y_res,left_right_corner_x]
,x_res
表示x轴上的分辨率例如经度分辨率,y_res
表示y轴上的分辨率例如纬度分辨率,left_right_corner
涉及到了投影原点的定义,这里是指栅格矩阵左上角像元(第一行第一列的像元)的左上角位置,这里不详细讨论)。(具体可查看:https://developers.google.com/earth-engine/guides/exporting_images?hl=zh-cn&authuser=1
);
只需要明白一点,如果使用地理坐标系输出且分辨率单位为°/度,就不要指定scale
而是定义好仿射系数crsTransform
,此外还需要指定输出范围region
参数(否则GEE
无法知晓输出的范围因为scale
未指定部分地理参数都没有进行计算);
- 对于掩膜后的数据集或者本身存在无效值的数据集,直接输出为tiff文件后
ArcGIS
软件能正确识别nan但是渲染时存在异常显示为黑色背景等问题,需要在输出时指定formatOptions
参数将掩膜的无效值定义为例如-9999,这样GIS软件可以正确识别其为无效值nan并且渲染其为透明色;(亦或者使用栅格计算器SetNull(IsNull("temp.tif"), "temp.tif")
,即将识别到的无效值设置为ArcGIS
中的无效值定义); - 对于输出分辨率过高,一张tiff文件可能存在的像元数超过一亿,那么就会超出默认的限制
maxPixels
而输出报错,这里将默认限制提高; crsTransform
、scale
和dimensions
是互斥的,即指定了其中某一个参数,则其他两个参数就不应指定;