結果示意圖:
GlobeLand30是30米空間分辨率全球地表覆蓋數據,目前可供下載的有3年的數據:2000-2010-2020。本文主要講解在ArcGIS10.6平台下進行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 !!!"
查看投影坐標系:
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. 掩膜提取
使用【按掩膜提取】工具,按研究區范圍進行提取。
配色參考值:
注意:拼接后的數據集是有金字塔的,但是使用【按掩膜提取】工具后又會丟失金字塔,必須重新構建金字塔,否則加載速度異常慢。
二、GlobeLand30成品數據下載
請私信作者(QQ:494535735),僅供學習交流使用。