ArcGIS风暴GlobeLand30全球数据处理教程（批量投影转换无效值处理拼接）

Posted 2021-06-16 刘一哥GIS

本文讲述GlobeLand30全球数据批处理流程，主要步骤包括：批量分幅投影转换、批量分幅无效值处理、批量图幅拼接和成品展示。由于图幅数目和数据量较大，本文是采用Python代码批量实现的，有关Python的基本操作及处理案例，可以参考实用精品专栏：《ArcGIS遇上Python》

文章目录

• 一、数据简介
• 二、数据准备
• 三、处理流程
• • 1. 批量投影转换
  • 2. 批量处理无效值
  • 3. 批量图幅拼接
  • 4. 添加分类字段并赋值
  • 5. 创建色彩映射表
  • 6. 导入色彩映射表
• 四、成品展示

---

一、数据简介

• GlobeLand30是30米空间分辨率全球地表覆盖数据，目前可供下载使用的有3年的数据：2000-2010-2020。

• GlobeLand30数据采用WGS-84坐标系。南纬85°-北纬85°之间的区域，投影方式采用UTM投影， 6度分带，坐标单位为米，坐标不加带号。南北纬85°-90°之间的区域，投影方式采用极地方位投影，投影面切于地球南北极点。

二、数据准备

下载3年（2000-2010-2020）数据，分别解压到对应的文件夹内。

三、处理流程

1. 批量投影转换

GlobeLand30数据采用WGS-84坐标系。

南纬85°-北纬85°之间的区域，投影方式采用UTM投影， 6度分带，坐标单位为米，坐标不加带号。

南北纬85°-90°之间的区域，投影方式采用极地方位投影，投影面切于地球南北极点。

因此需要先对所有的图幅进行投影转换，然后进行后续拼接等操作。本文将其转换为地理坐标系为WGS1984。

Python批处理核心代码：

def GetRaster(path):
    n=0
    if (os.path.exists(path)):
        files = os.listdir(path)
        for file in files:
            m = os.path.join(path,file)
            if (os.path.isdir(m)):
                h=os.path.split(m)
                ws=root+"/"+h[1]
                arcpy.env.workspace = ws
                rasters = arcpy.ListRasters("*", "tif")
                
                for raster in rasters:
                    print "Finding  "+ raster
                    
                    print "Projecting "+raster
                    out=outPath+"/"+raster
                    arcpy.ProjectRaster_management(raster,out, projecttype, "NEAREST",cellsize , "", "",geogratype)
                    n=n+1
    print str(n)+" rasters are processed !!!"
GetRaster(root)

ArcGIS Python支持多线程处理，3年数据可以同步进行批处理。。。。

由于图幅数目和数据量均较大，该过程将持续很长的时间，可能是几个小时或者几天，因此呢，如果你没有足够强大的机器和心脏，建议不要轻易尝试。

2000年：

2010年：

2020年：

北半球投影结果：

局部图幅展示：

投影转换后的各个图幅，坐标系为WGS 1984。

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!"

3. 批量图幅拼接

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

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

方法二：Python代码

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

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

4. 添加分类字段并赋值

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

分类赋值参照表：

5. 创建色彩映射表

双击2000年数据图层，打开图层属性，切换到【符号系统】，点击【色彩映射表】→【导出色彩映射表】。

选择存放路径，并输入名称，色彩映射表的后缀名是.clr。

6. 导入色彩映射表

有了刚才导出的色彩映射表（GlobeLand30色彩映射表.clr文件），可以快速地对剩余两年的数据做符号化了，方法是：双击2010和2020年数据图层，打开图层属性，切换到【符号系统】，点击【色彩映射表】→【导入色彩映射表】。

四、成品展示