最新数据处理 之 最新DSM（AW3D30）数据批量拼接--文末附数据获取方式

本文主要针对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景批量拼接的结果部分显示，可以看到没有拼接痕迹。

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

常用免费DEM数据汇总（含下载使用方法）

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

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