结果示意图:

GlobeLand30是30米空间分辨率全球地表覆盖数据,目前可供下载的有3年的数据:2000-2010-2020。本文主要讲解在ArcGIS10.6平台下进行GlobeLand30的预处理操作,主要预处理步骤包括:批量分幅投影转换批量分幅无效值处理批量图幅拼接数据集掩膜提取

目录

  • 一、GlobeLand30 ArcGIS批处理完整流程
  • 1. 批量投影转换
  • 2. 批量处理无效值
  • 3. 批量图幅拼接
  • 4. 添加分类字段并赋值
  • 5. 掩膜提取
  • 二、GlobeLand30成品数据下载

一、GlobeLand30 ArcGIS批处理完整流程

GlobeLand30的预处理步骤主要有:批量分幅投影转换批量分幅无效值处理批量图幅拼接数据集掩膜提取

1. 批量投影转换

由于原数据集使用的是6°带UTM投影,图幅从43-53跨了11个带,因此需要先对所有的图幅进行投影转换,然后进行后续操作。为了便于精确计算面积,本文将其转换为Albers等积投影,地理坐标系为WGS1984。

Python源代码:

由于原始数据集分幅是保存在不同的文件夹,因此第一步批量投影变换需要遍历根目录下的所有文件夹,过程稍微有些复杂,从第二步开始,将所有的处理结果保存在一个文件夹内,减少了代码量。

arcpy.ProjectRaster_management(raster,out, "PROJCS['MyAlbers',GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Albers'],PARAMETER['False_Easting',0.0],PARAMETER['False_Northing',0.0],PARAMETER['central_meridian',105.0],PARAMETER['Standard_Parallel_1',25.0],PARAMETER['Standard_Parallel_2',47.0],PARAMETER['latitude_of_origin',0.0],UNIT['Meter',1.0]]", "NEAREST", cellsize , "", "", "GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]")
                    n=n+1
    print str(n)+" rasters are processed !!!"

查看投影坐标系:

全球10米土地覆盖产品 ESRI 全球30米地表覆盖数据(globeland30)_全球10米土地覆盖产品 ESRI

2. 批量处理无效值

第一步投影结束后,影像的黑边依然存在,这将影像后面的拼接处理结果,黑边对应的像元值为0,只需要将其设置为NoData即可,对应的工具是栅格计算器中的SetNull函数,我们继续使用Python批处理。

Python源码:

for inRaster in inRasters:
    outRaster=outpath+inRaster
    ##print outRaster
    print "Processing "+inRaster +" ......"
    outSetNull = SetNull(inRaster,inRaster,WhereClause)
    outSetNull.save(outRaster)
    n=n+1
print str(n)+" rasters are processed!"

处理过程:

全球10米土地覆盖产品 ESRI 全球30米地表覆盖数据(globeland30)_数据集_02

处理结果:

全球10米土地覆盖产品 ESRI 全球30米地表覆盖数据(globeland30)_Python_03

3. 批量图幅拼接

在ArcGIS中拼接栅格图幅有两个工具:【镶嵌】和【镶嵌至新栅格】,当然了需要手动添加数据,好在我们的数据在同一个文件夹下,可以直接全选拖动,当然了,也可以写Python代码来完成。

方法一:使用【镶嵌至新栅格】工具

全球10米土地覆盖产品 ESRI 全球30米地表覆盖数据(globeland30)_d3_04

方法二:Python代码

arcpy.MosaicToNewRaster_management(str, outPath, outName, "", "8_BIT_UNSIGNED", "", "1", "LAST", "FIRST")

处理过程:

全球10米土地覆盖产品 ESRI 全球30米地表覆盖数据(globeland30)_d3_05

需要注意的是,原始数据集是tif格式,且每个图幅都是有金字塔的,加载和显示速度都比较快,但是SetNull后会丢失金字塔,加载速度很慢。

4. 添加分类字段并赋值

打开属性表,添加一个文本型字段type,打开编辑器,输入分类名称并保存。

全球10米土地覆盖产品 ESRI 全球30米地表覆盖数据(globeland30)_数据集_06

5. 掩膜提取

使用【按掩膜提取】工具,按研究区范围进行提取。

全球10米土地覆盖产品 ESRI 全球30米地表覆盖数据(globeland30)_全球10米土地覆盖产品 ESRI_07


配色参考值:

全球10米土地覆盖产品 ESRI 全球30米地表覆盖数据(globeland30)_d3_08

注意:拼接后的数据集是有金字塔的,但是使用【按掩膜提取】工具后又会丢失金字塔,必须重新构建金字塔,否则加载速度异常慢。

二、GlobeLand30成品数据下载