本文主要针对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数据汇总与介绍,感兴趣的可以下载,自行处理一下看看。
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")
这个代码实际上就是上述单个文件夹数据拼接的批量版,都是为了实现批量拼接,减少人工的操作。