基于ArcPy的遥感数据拼接——以中国90米DEM为例

最新推荐文章于 2024-04-25 21:44:44 发布

原创最新推荐文章于 2024-04-25 21:44:44 发布 · 1k 阅读

9 ·

CC 4.0 BY-SA版权

文章标签：

#ArcPy #ArcGIS #Python

Arcpy 专栏收录该内容

0 篇文章

订阅专栏

本文详细介绍使用ArcPy脚本批量处理DEM数据的投影转换，以及如何利用MosaicToNewRaster工具实现数据拼接和裁切。通过代码示例，展示了从数据下载到最终成果的全过程。

1.数据下载——DEM

地理空间数据云：网页直通车
地理空间数据云主页

1.1选择DEM数字高程数据：

操作步骤

1.2数据下载

根据研究区的范围下载数据

2.对于遥感数据的投影处理

由于拼接中国整幅DEM所用到的图幅较多，因此采用ArcPy脚本进行批量化的遥感图像投影处理，代码如下：

#Project_Batch.py
import arcpy
def inputstring(list_Rasters):
    #用来转化列表类型，并且添加";"分隔
    inputstr=""
    for instr in range(len(list_Rasters)):
        if instr != len(list_Rasters)-1:
            inputstr+=list_Rasters[instr]+";"
        else:
            inputstr+=list_Rasters[instr]
    return inputstr
arcpy.env.workspace="D:\\遥感解译\\遥感"#工作空间
Input_Rasters=arcpy.ListDatasets()#遍历所有的Datasets()图像
Output_Coordinate_System =open("D:\\project.txt","r").read()#读入坐标文件
Old_Project=open("D:\\Old_Project.txt",',').read()#读入原始坐标文件
ProjectRaster = "D:\\遥感解译\\"
for raster in Input_Rasters:
 raster=str(raster)
 arcpy.ProjectRaster_management(raster,ProjectRaster+raster, Output_Coordinate_System, "NEAREST", /
 "1.01469440773919E-03  1.01469440773919E03","","",Old_Project)
print("定义完成")#提示定义完成
'''#Old_Project.txt
"PROJCS['WGS_1984_UTM_Zone_43N',GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Transverse_Mercator'],PARAMETER['False_Easting',500000.0],PARAMETER['False_Northing',0.0],PARAMETER['Central_Meridian',75.0],PARAMETER['Scale_Factor',0.9996],PARAMETER['Latitude_Of_Origin',0.0],UNIT['Meter',1.0]]"
#Output_Coordinate_System.txt
"GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]"
'''

之所以进行重投影，首先因为所下载数据’WGS_1984_UTM_Zone_43N’——通用横轴墨卡托投影，具有不同的中央纬线，因此将其转化为WGS_1984地理坐标系统，使得数据可以进行裁切。

3.使用MosaicToNewRaster工具进行数据拼接

#MosaicToClip
#批处理前提，首先所有要素在一个投影之下，才能去拼接呀！
import arcpy
def inputstring(list_Rasters):
    #用来转化列表类型，并且添加";"分隔
    inputstr=""
    for instr in range(len(list_Rasters)):
        if instr != len(list_Rasters)-1:
            inputstr+=list_Rasters[instr]+";"
        else:
            inputstr+=list_Rasters[instr]
    return inputstr
arcpy.env.workspace="D:\\遥感解译\\No2-07-20180127\\GDEMV2_30M_Hebei"
Input_Rasters=arcpy.ListDatasets()
Out_Loction="D:\\遥感解译\\RS"
Rater_Extension="Mosaic_China"
inputstr=inputstring(Input_Rasters)
arcpy.MosaicToNewRaster_management(inputstr,Out_Loction,Rater_Extension,"","8_BIT_UNSIGNED","","1","MAXIMUM","FIRST")
#注意所访问要素的格式
Input_Raster="D:\\遥感解译\\Mosaic_China"
Rectangle="-2717252.8552 2207537.7556 2128604.968 6224070.0716"
Output_Extent="D:\\遥感解译\\clip.shp"
Output_Raster_Dataset="D:\\遥感解译\\计算机地图制图\\China.gdb\\China_DEM.tif"
Use_Input_Features_for_Clipping_Geometry=True#使用多边形几何形状进行裁切
arcpy.Clip_management(Input_Raster, Rectangle, Output_Raster_Dataset, Output_Extent,\
                      "", Use_Input_Features_for_Clipping_Geometry, "NO_MAINTAIN_EXTENT")
print("拼接裁切完成")

如果对ArcPy使用较为熟练，可以将以上的两端代码组合起来使用。

4.总结

经过本人的测试，ArcGIS的MosaicToNewRaster，对于图像拼接的算法相较于Erdas的Mosaic工具运行较慢。因此我推荐可以直接定义完投影后直接使用Erdas的Mosaic模块。下面介绍Mosaic的用法：
Data Preparation——Mosaic Images——Mosaic Tool

Erdas裁切模块的使用
效果展示