本文主要针对JAXA发布的最新版本3.1&3.2(AW3D30)数据进行数据拼接处理。

JAXA利用ALOS的PRISM立体模式制作了全球的DSM数据,并免费分布了30米版本产品。详细的数据介绍和下载:AW3D30

与DEM数据一样,数据是按照格网大小进行分发的,如 N028E090.zip,就是28N, 90E区域的DSM数据。

数据下载完成后,解压出来,里面包括了 "_DSM.tif","_MSK.tif","_STK.tif" 三种tif格式的数据,这里 DSM 数据为例介绍批量拼接。

数据的拼接,在ArcGIS软件中可以很容易的实现,但如果数据过多的话,在ArcGIS中点击会比较繁琐。因此,考虑用 ArcPy 进行批量拼接。废话不多说,代码如下:

# -*- coding: utf-8 -*-
from osgeo import gdal

arcpy.CheckOutExtension("Spatial")

arcpy.env.workspace = "tiff"  # 存放所有数据的目录
outFolder = "output"  # 输出路径
outName = "DEM.tif"  # 输出名称

base = "ALPSMLC30_N028E090_DSM.tif"  # 基础数据, 主要用来提取基本参数 (随便选一幅就可以)

# 获取基本参数
out_coor_system = arcpy.Describe(base).spatialReference
dataType = arcpy.Describe(base).DataType
piexl_type = arcpy.Describe(base).pixelType
cellwidth = arcpy.Describe(base).meanCellWidth
bandcount = arcpy.Describe(base).bandCount

# 获取待拼接的所有数据名称
rasters = []
for ras in arcpy.ListRasters("*DSM.tif"):
    rasters.append(ras)

ras_list = ";".join(rasters)

print ras_list  # 检查名称

# 批量拼接
arcpy.MosaicToNewRaster_management(ras_list, outFolder, outName, out_coor_system, "16_BIT_SIGNED", cellwidth, bandcount, "LAST", "FIRST")

代码经过了测试,没有问题。下图是60景批量拼接的结果部分显示,可以看到没有拼接痕迹。

gis的esriAddIn怎么打开 arcgis打开dsm_gdal

 对于其他的栅格数据,理论上可以用ArcGIS实现的拼接操作,这个代码稍加修改都可以实现,因为本身就是调用的ArcGIS拼接工具。下图是常用的DEM数据汇总与介绍,感兴趣的可以下载,自行处理一下看看。


gis的esriAddIn怎么打开 arcgis打开dsm_dsm_02

2021.12.24补充

上述代码能够实现同一文件夹内批量数据的拼接,但如果有多个文件夹怎么处理?都放在一个文件夹里,显然操作比较繁琐,而且拼接的时候,处理时间长,需要占用大量的内存。因此,如果数据分布在多个文件夹内,每个文件夹内的数据都需要拼接,那么可以考虑使用下面批量拼接代码:

# -*- coding: utf-8 -*-
import os
import arcpy

dem_path = "tiff"  # 存放数据目录, 如 tiff/N045E090/ALPSMLC30_N049E094_DSM.tif
outFolder = "output"  # 指定输出文件夹

for file in os.listdir(dem_path):
    arcpy.env.workspace = dem_path + '/' + file  # 指定工作目录, 即存放影像的目录

    # 基础数据, 主要用来提取基本参数
    base = os.listdir(dem_path + '/' + file)[1]

    out_coor_system = arcpy.Describe(base).spatialReference
    dataType = arcpy.Describe(base).DataType
    piexl_type = arcpy.Describe(base).pixelType
    cellwidth = arcpy.Describe(base).meanCellWidth
    bandcount = arcpy.Describe(base).bandCount

    arcpy.CheckOutExtension("Spatial")

    # 获取待拼接的所有数据名称
    rasters = []
    for ras in arcpy.ListRasters("*DSM.tif"):
        rasters.append(ras)

    ras_list = ";".join(rasters)

    print(ras_list)  # 检查名称

    outName = file + ".tif"  # 输出名称, 以文件夹命名

    arcpy.MosaicToNewRaster_management(ras_list, outFolder, outName,
                                       out_coor_system,
                                       "16_BIT_SIGNED", cellwidth,
                                       bandcount, "LAST", "FIRST")

这个代码实际上就是上述单个文件夹数据拼接的批量版,都是为了实现批量拼接,减少人工的操作。