1. 为什么我要自己动手做ASTER GDEM镶嵌
第一次接触ASTER GDEM是在做一个区域地形分析的小项目,当时需要覆盖三个省份的连续高程数据。我一开始想得很简单,直接去下载现成的分幅数据,用GIS软件拼一下不就完了。结果打开一看,每幅数据大概1度乘1度的范围,光目标区域就涉及几十幅,而且相邻图幅之间还有重叠带和边缘噪声。手动拼了几次,接边处要么出现明显的亮度断层,要么高程值在重叠区对不上,做坡度分析的时候直接出现条带状异常。
后来我才认真去研究ASTER GDEM这个数据源本身。它是由搭载在Terra卫星上的ASTER传感器获取的立体像对生成的全球数字高程模型,第一版在2009年发布,第二版在2011年更新,空间分辨率是1弧秒,也就是大约30米。这个数据覆盖了北纬83度到南纬83度之间的所有陆地,对于没有实测地形数据的区域来说,几乎是唯一能免费拿到的中高分辨率高程来源。但它的原始分发方式就是按1度乘1度的分幅打包,全球一共两万两千多幅,你要用哪块就得自己镶嵌。
镶嵌这件事听起来不复杂,但真正做起来坑不少。首先是数据量的问题,一幅GeoTIFF大概二十多兆,几十幅加起来就是几个G,如果涉及上百幅,普通笔记本的内存和磁盘读写就会吃紧。其次是接边处理,ASTER GDEM不同图幅之间由于成像时间、大气条件、传感器姿态的差异,边缘区域的高程值会有系统性偏差,直接拼接会出现明显的接缝。再就是无效值处理,水体、云覆盖、陡峭阴影区域会出现-9999这样的填充值,如果不做掩膜和插补,后续分析全废。
我踩过的最大的坑是第一次镶嵌完没有做投影统一。ASTER GDEM原始数据是地理坐标系,经纬度以WGS84为基准,但如果你拿到的分幅数据里混入了不同投影版本,或者你自己在镶嵌过程中不小心改了投影参数,拼出来的结果就会出现错位。我当时有一幅数据不知道什么时候被转成了UTM,结果那一块整体偏移了几百米,做等高线的时候直接穿帮。
所以这篇文章我想把整个ASTER GDEM镶嵌的流程从头到尾讲清楚,包括数据准备、投影统一、接边处理、无效值填充、镶嵌后裁剪和质检。适合谁看呢,我觉得三类人比较合适:一是做区域地形分析的学生和研究人员,二是做遥感影像处理的技术人员,三是需要高程数据做三维建模或者水文分析的朋友。哪怕你之前没碰过GDEM,跟着走一遍也能把流程跑通。
2. 镶嵌前的整体设计与核心思路拆解
2.1 为什么不能直接拼接而要设计流程
很多人觉得镶嵌就是把几幅图放到一起,软件里点一下拼接就完事了。但ASTER GDEM的分幅数据有几个特殊性质,决定了你不能这么干。
第一,分幅之间是有重叠的。ASTER GDEM的每一幅标准分幅是1度乘1度,但实际有效数据范围会略大于这个边界,相邻图幅之间通常有几十个像素的重叠带。这个重叠带的存在本来是为了方便接边,但如果你直接拼接,重叠区域就会出现两套高程值,软件默认取其中一个,接缝处就会出现突变。
第二,不同图幅的高程基准有微小差异。虽然ASTER GDEM整体是基于WGS84椭球的高程,但不同图幅在生成过程中使用的控制数据不同,导致图幅之间可能存在米级的系统性偏差。这个偏差在平坦地区不明显,但在坡度分析、水文流向计算中会被放大。
第三,无效值分布不均匀。水体、云、雪、阴影区域都会产生无效值,而且不同图幅的无效值分布不一样。如果你先镶嵌再处理无效值,会发现有些接缝处的无效值连成一片,很难判断是数据本身的问题还是拼接引入的。
所以我的整体思路是:先做单幅数据的预处理,包括投影统一、无效值标记、边缘噪声裁剪;然后做接边处理,通过重叠区的高程统计来调整图幅之间的偏差;最后才做镶嵌和裁剪。这个顺序不能乱,乱了就要返工。
2.2 工具选型:为什么我最终用了GDAL加Python
工具选型这块我试过不少方案。ArcGIS的镶嵌工具确实方便,Mosaic To New Raster点几下就能出结果,但问题是它对大批量数据的处理效率不高,而且接边处理的参数不够透明,你很难控制它到底怎么取重叠区的值。QGIS的虚拟栅格功能适合做快速预览,但真正要输出一个干净的结果,还是得靠命令行工具。
我最终选的是GDAL加Python脚本的方案。GDAL的gdalwarp可以做投影转换和重采样,gdalbuildvrt可以快速构建虚拟镶嵌,gdal_translate可以把虚拟镶嵌转成实体文件。Python的rasterio和numpy可以让我灵活地读取重叠区的高程值,做统计分析和偏差调整。这个组合的好处是每一步都可控,参数都在你手里,出了问题也能定位到具体是哪一步。
提示:如果你只是做一次性的小区域镶嵌,用QGIS的虚拟栅格加导出就够了。但如果涉及几十幅以上的数据,或者需要反复调整接边参数,建议还是走脚本路线。
2.3 投影选择与重采样方法
ASTER GDEM原始数据是地理坐标系,单位是度。做镶嵌的时候,我建议统一转成投影坐标系,单位是米。原因很简单,地理坐标系下不同纬度的1度对应的实际距离不一样,你在做坡度、距离、面积计算的时候会引入纬度相关的误差。转成投影之后,单位统一,计算也方便。
投影的选择要看你的研究区域。如果区域不大,比如一个省的范围,用UTM分区投影就够了。如果区域跨多个UTM带,可以用Albers等面积投影或者Lambert等角投影。我一般用Albers,因为它在中等纬度地区的面积变形小,适合做区域分析。
重采样方法我试过最近邻、双线性、三次卷积。最近邻速度快,但会引入阶梯状的高程突变,做坡度分析的时候会出现明显的块状噪声。双线性平滑一些,但在陡峭地形区域会过度平滑。三次卷积效果最好,但计算量大。我的建议是:如果做坡度、水文分析,用双线性;如果做三维可视化,用三次卷积;如果只是做粗略预览,最近邻也行。
# 投影转换示例:将地理坐标系转为Albers投影 gdalwarp -t_srs "+proj=aea +lat_1=25 +lat_2=47 +lat_0=36 +lon_0=105 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs" \ -tr 30 30 -r bilinear \ input_gdem.tif output_gdem_proj.tif这段命令里,-t_srs指定目标投影,-tr 30 30指定输出分辨率为30米,-r bilinear指定双线性重采样。注意-tr后面的两个30是x和y方向的分辨率,单位是米,和投影单位一致。
3. 核心细节解析与实操要点
3.1 数据下载与命名规范
ASTER GDEM的数据可以从多个渠道获取,我常用的是NASA的EarthData平台。下载的时候要注意,每一幅数据都有唯一的标识码,命名格式一般是ASTGTMV003_NxxExxx这样的结构,其中N后面跟纬度,E后面跟经度。下载下来是压缩包,解压后有两个文件,一个是dem文件,一个是num文件。dem就是高程数据,num是质量评估文件,记录了每个像素的有效性。
我建议在下载阶段就做好命名规范。比如你研究区域是北纬30到33度,东经100到103度,那涉及的分幅就是N30E100到N32E102,一共九幅。你可以按行列号重命名,比如row30_col100.tif,这样后面做镶嵌的时候排序和索引都方便。
注意:下载的时候一定要确认版本号。ASTER GDEM有V1、V2、V3三个版本,V3是最新的,修正了V2中的一些异常值。如果你拿到的数据里混了不同版本,接边处的高程偏差会更大。
3.2 无效值处理与质量掩膜
ASTER GDEM的无效值通常用-9999表示。这些无效值出现在水体、云覆盖、传感器阴影、陡峭地形等区域。如果你不做处理,镶嵌后的结果里就会有大片的-9999,做任何分析都会出错。
我的处理方法是分两步。第一步,用num文件生成质量掩膜。num文件里每个像素的值代表了该像素的有效性等级,一般大于0的值表示有效,等于0表示无效。你可以用gdal_calc把num文件转成二值掩膜,然后应用到dem上。第二步,对于掩膜后的无效区域,用邻域插值或者全局插值填充。邻域插值适合小面积的无效区,全局插值适合大面积的缺失。
import rasterio import numpy as np from scipy import ndimage # 读取dem和num文件 with rasterio.open('dem.tif') as src: dem = src.read(1) profile = src.profile with rasterio.open('num.tif') as src: num = src.read(1) # 生成有效掩膜 valid_mask = (num > 0) & (dem != -9999) # 对无效区域做邻域均值填充 filled_dem = dem.copy() filled_dem[~valid_mask] = np.nan # 用3x3窗口的均值填充 kernel = np.ones((3,3)) / 9 neighbor_mean = ndimage.convolve(np.nan_to_num(filled_dem), kernel, mode='constant') filled_dem[~valid_mask] = neighbor_mean[~valid_mask] # 如果还有剩余的nan,用全局均值填充 global_mean = np.nanmean(filled_dem) filled_dem[np.isnan(filled_dem)] = global_mean # 写出结果 with rasterio.open('dem_filled.tif', 'w', **profile) as dst: dst.write(filled_dem, 1)这段代码的逻辑是先用num文件生成有效掩膜,然后把无效区域设为nan,用3x3窗口的邻域均值填充,最后用全局均值兜底。实际用的时候,窗口大小可以根据无效区的面积调整,大面积缺失可以用更大的窗口或者用反距离加权插值。
3.3 接边处理的核心逻辑
接边处理是镶嵌里最关键的环节。我的做法是:先找到相邻图幅的重叠区,然后在重叠区里统计两幅图的高程差异,如果差异超过阈值,就对其中一幅做整体偏移调整。
具体步骤是这样的。假设你有两幅图A和B,它们在水平方向上有重叠。你先用gdalwarp把两幅图都转到同一个投影和分辨率下,然后用gdalbuildvrt构建一个虚拟镶嵌,再用Python读取重叠区的像素值。计算A和B在重叠区的平均高程差,如果差值大于5米,就把B的整体高程加上这个差值。这个5米的阈值是我根据经验定的,平坦地区可以设小一点,山区可以设大一点。
提示:接边调整不要一次性对所有图幅做,建议按行列顺序逐对处理。先处理第一行,再处理第二行,每一行内部从左到右依次调整。这样可以把累积误差控制住。
3.4 边缘噪声裁剪
ASTER GDEM的每一幅数据边缘都有大约1到2个像素的噪声带,表现为高程值的突变或者无效值。如果你不裁剪,镶嵌后的接缝处就会出现细条状的异常。我的做法是在镶嵌前对每一幅数据做边缘裁剪,裁掉最外圈的2个像素。
# 边缘裁剪示例:裁掉最外圈2个像素 gdal_translate -srcwin 2 2 $(expr $(gdalinfo input.tif | grep "Size is" | awk '{print $3}') - 4) $(expr $(gdalinfo input.tif | grep "Size is" | awk '{print $4}') - 4) input.tif output_trimmed.tif这个命令用gdal_translate的-srcwin参数指定裁剪窗口,从第2行第2列开始,宽度和高度各减去4个像素。实际用的时候,你可以先打开一幅数据看看边缘噪声的宽度,再决定裁多少。
4. 完整实操流程与关键环节实现
4.1 环境准备与依赖安装
我用的环境是Ubuntu 20.04加Python 3.8。GDAL的安装建议用conda,因为pip安装的GDAL经常和系统库版本不匹配。rasterio和numpy也用conda装,scipy用来做插值。
conda create -n gdem python=3.8 conda activate gdem conda install -c conda-forge gdal rasterio numpy scipy装完之后验证一下GDAL版本,确保gdalwarp和gdalbuildvrt都能正常调用。
gdalinfo --version如果输出类似GDAL 3.4.1这样的版本信息,就说明装好了。
4.2 批量投影转换与重采样
假设你已经下载了九幅数据,放在data目录下,命名是N30E100.tif到N32E102.tif。第一步是批量做投影转换。
#!/bin/bash for file in data/*.tif; do basename=$(basename $file .tif) gdalwarp -t_srs "+proj=aea +lat_1=25 +lat_2=47 +lat_0=36 +lon_0=105 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs" \ -tr 30 30 -r bilinear \ $file proj/${basename}_proj.tif done这个脚本会遍历data目录下的所有tif文件,逐个做投影转换,输出到proj目录。注意-tr 30 30指定了输出分辨率,这个值要和你的研究需求匹配。如果你需要更高分辨率,可以设成15或者10,但那样数据量会翻倍。
4.3 无效值填充与边缘裁剪
投影转换完之后,对每一幅数据做无效值填充和边缘裁剪。我把这两步合并到一个Python脚本里。
import rasterio import numpy as np from scipy import ndimage import os def fill_and_trim(input_path, output_path, trim_pixels=2): with rasterio.open(input_path) as src: dem = src.read(1) profile = src.profile # 无效值填充 valid_mask = dem != -9999 filled_dem = dem.copy() filled_dem[~valid_mask] = np.nan # 邻域均值填充 kernel = np.ones((3,3)) / 9 neighbor_mean = ndimage.convolve(np.nan_to_num(filled_dem), kernel, mode='constant') filled_dem[~valid_mask] = neighbor_mean[~valid_mask] # 全局均值兜底 global_mean = np.nanmean(filled_dem) filled_dem[np.isnan(filled_dem)] = global_mean # 边缘裁剪 trimmed_dem = filled_dem[trim_pixels:-trim_pixels, trim_pixels:-trim_pixels] # 更新profile profile.update({ 'height': trimmed_dem.shape[0], 'width': trimmed_dem.shape[1], 'transform': rasterio.transform.from_origin( src.transform.c + trim_pixels * src.transform.a, src.transform.f + trim_pixels * src.transform.e, src.transform.a, -src.transform.e ) }) with rasterio.open(output_path, 'w', **profile) as dst: dst.write(trimmed_dem, 1) # 批量处理 for file in os.listdir('proj'): if file.endswith('_proj.tif'): input_path = os.path.join('proj', file) output_path = os.path.join('filled', file.replace('_proj.tif', '_filled.tif')) fill_and_trim(input_path, output_path)这个脚本里,trim_pixels=2表示裁掉最外圈2个像素。rasterio.transform.from_origin用来更新裁剪后的地理变换参数,确保裁剪后的数据还能正确定位。
4.4 接边偏差统计与调整
接边调整我写了一个单独的脚本,逻辑是遍历所有相邻图幅对,计算重叠区的高程差,然后做偏移调整。
import rasterio import numpy as np import os def get_overlap_diff(file_a, file_b): with rasterio.open(file_a) as src_a: dem_a = src_a.read(1) transform_a = src_a.transform crs_a = src_a.crs with rasterio.open(file_b) as src_b: dem_b = src_b.read(1) transform_b = src_b.transform crs_b = src_b.crs # 确保投影一致 if crs_a != crs_b: return None # 计算重叠区 # 这里简化处理,假设两幅图在水平方向相邻 # 实际用的时候需要用rasterio的window功能精确计算 overlap_width = min(dem_a.shape[1], dem_b.shape[1]) // 10 overlap_a = dem_a[:, -overlap_width:] overlap_b = dem_b[:, :overlap_width] # 计算平均高程差 diff = np.nanmean(overlap_a) - np.nanmean(overlap_b) return diff # 遍历相邻图幅 files = sorted([f for f in os.listdir('filled') if f.endswith('.tif')]) for i in range(len(files) - 1): diff = get_overlap_diff( os.path.join('filled', files[i]), os.path.join('filled', files[i+1]) ) if diff is not None and abs(diff) > 5: print(f"调整 {files[i+1]},偏移量:{diff:.2f}米") # 实际调整逻辑:读取文件,加上偏移量,重新写出这个脚本是个简化版,实际用的时候需要用rasterio的window功能精确计算重叠区,而不是简单取边缘。但核心逻辑是一样的:找到重叠区,算平均差,超过阈值就调整。
4.5 虚拟镶嵌与实体输出
所有单幅数据处理好之后,用gdalbuildvrt构建虚拟镶嵌,再用gdal_translate转成实体文件。
# 构建虚拟镶嵌 gdalbuildvrt mosaic.vrt filled/*.tif # 转成实体GeoTIFF gdal_translate -of GTiff -co COMPRESS=LZW -co TILED=YES mosaic.vrt mosaic.tif-co COMPRESS=LZW指定用LZW压缩,可以显著减小文件体积。-co TILED=YES指定用分块存储,方便后续读取。
4.6 镶嵌后裁剪与质检
镶嵌完的mosaic.tif覆盖了你下载的所有图幅范围,但你可能只需要其中一部分。用gdalwarp加-cutline或者用gdal_translate加-projwin做裁剪。
# 按经纬度范围裁剪 gdal_translate -projwin 100.0 33.0 103.0 30.0 mosaic.tif final.tif-projwin后面跟的是左上角经度、左上角纬度、右下角经度、右下角纬度。注意这个顺序不能错,否则裁出来是反的。
质检我一般做三件事。第一,用gdalinfo看统计信息,确认最大值、最小值、均值在合理范围内。第二,用QGIS打开看接缝处有没有明显的条带。第三,随机取几个点,和已知的高程数据做对比。
提示:质检的时候特别要注意接缝处的高程突变。如果发现某条接缝有明显的亮度差异,说明接边调整没做好,需要回去重新处理那一对图幅。
5. 常见问题与排查技巧实录
5.1 接缝处出现条带状异常怎么办
这是最常见的问题。原因通常有三个:一是接边偏差没调整,二是边缘噪声没裁剪,三是重采样方法不一致。
排查顺序是这样的。先检查所有图幅的投影和分辨率是否一致,用gdalinfo逐个看。如果投影不一致,回去重新做投影转换。如果投影一致,检查边缘裁剪是否都做了,裁了多少像素。如果都没问题,那就是接边偏差的问题,需要重新计算重叠区的高程差,调整阈值。
我遇到过一次,接缝处出现了一条宽度约5个像素的亮带。查了半天发现是其中一幅数据在做投影转换的时候用了最近邻重采样,其他用的是双线性。最近邻在边缘产生了阶梯状的高程突变,和双线性的平滑结果拼在一起就出现了亮带。统一重采样方法之后就解决了。
5.2 镶嵌后数据量太大怎么优化
几十幅数据镶嵌完,文件可能有好几个G。优化方法有几个。第一,用LZW或者DEFLATE压缩,一般能压到原来的三分之一。第二,用分块存储,TILED=YES,这样读取的时候不用加载整个文件。第三,如果研究区域不大,先裁剪再镶嵌,而不是先镶嵌再裁剪。第四,如果不需要原始分辨率,可以重采样到更低的分辨率,比如从30米降到90米,数据量直接降到九分之一。
# 压缩和分块存储 gdal_translate -of GTiff -co COMPRESS=DEFLATE -co PREDICTOR=2 -co TILED=YES -co BIGTIFF=IF_NEEDED mosaic.tif mosaic_compressed.tif-co PREDICTOR=2对高程数据特别有效,因为它利用了相邻像素的相关性。-co BIGTIFF=IF_NEEDED确保文件超过4G的时候自动用BigTIFF格式。
5.3 无效值填充后出现异常高值或低值
这个问题通常出现在大面积无效区填充的时候。邻域均值填充在边缘区域会引入异常值,因为邻域里可能包含无效值或者极端值。我的解决方法是:先用有效掩膜把无效区标记出来,然后只用有效像素做邻域统计。如果无效区面积超过窗口大小,就扩大窗口或者改用反距离加权插值。
还有一个技巧是用中值滤波代替均值滤波。中值滤波对极端值不敏感,填充结果更稳健。但中值滤波的计算量比均值滤波大,大面积填充的时候会慢一些。
5.4 常见问题速查表
| 问题现象 | 可能原因 | 排查方法 | 解决方案 |
|---|---|---|---|
| 接缝处条带状异常 | 接边偏差未调整 | 检查重叠区高程差 | 重新计算偏差并调整 |
| 接缝处细线状异常 | 边缘噪声未裁剪 | 检查边缘像素值 | 裁剪最外圈2-3像素 |
| 镶嵌后文件过大 | 未压缩或未分块 | 查看文件属性 | 用LZW压缩和TILED存储 |
| 无效值填充后异常 | 邻域包含无效值 | 检查填充区域 | 用有效掩膜限制邻域 |
| 投影错位 | 投影未统一 | gdalinfo检查CRS | 统一投影后重新镶嵌 |
| 坡度分析出现块状噪声 | 重采样方法不当 | 检查重采样参数 | 改用双线性或三次卷积 |
5.5 几个我踩过的坑和独家技巧
第一个坑是下载数据的时候没注意版本号。我有一次混用了V2和V3的数据,接边处的高程差达到了十几米,调了半天才发现是版本问题。后来我养成了习惯,下载的时候先看文件名里的版本标识,统一用V3。
第二个坑是投影转换的时候没设分辨率。gdalwarp默认会保持原始分辨率,但地理坐标系下的1度对应的实际距离和投影坐标系下的30米不是一回事。如果你不设-tr参数,转换后的分辨率可能是乱的。我一般显式设成30 30,确保输出分辨率一致。
第三个技巧是接边调整的时候用中位数而不是平均值。平均值容易受极端值影响,中位数更稳健。特别是在山区,重叠区里可能有陡峭地形导致的高程突变,用中位数可以避免这些突变影响整体偏移量。
第四个技巧是镶嵌前先做一个低分辨率的预览。用gdalbuildvrt加-resolution lowest参数,快速生成一个低分辨率镶嵌,看看整体效果。如果预览里就有明显的接缝,那高分辨率的结果肯定也有问题,可以提前发现。
第五个技巧是保留中间文件。投影转换后的文件、填充后的文件、裁剪后的文件,都保留着。万一后面发现某一步有问题,不用从头再来。磁盘空间不够的话,至少保留投影转换后的文件,因为那一步最耗时。
6. 镶嵌后的数据应用与扩展思路
6.1 基于镶嵌结果的地形分析
镶嵌完的GDEM数据可以直接用来做坡度、坡向、地形起伏度分析。我用GDAL的gdaldem工具做坡度分析,一条命令就能出结果。
# 坡度分析 gdaldem slope final.tif slope.tif -p -s 1.0 # 坡向分析 gdaldem aspect final.tif aspect.tif # 地形起伏度(用3x3窗口的极差) gdaldem roughness final.tif roughness.tif-p表示输出以度为单位的坡度,-s 1.0表示垂直方向的缩放系数。如果你的高程数据单位是米,投影单位也是米,这个系数就是1.0。如果投影单位是度,就需要调整这个系数。
坡度分析的结果可以用来做地形分级、水土流失评估、选址分析等。我做过一个区域的地形起伏度分析,用镶嵌后的GDEM算出来的结果和实地考察的数据对比,误差在可接受范围内。
6.2 水文分析中的流向计算
GDEM镶嵌结果在水文分析里用得很多,比如计算流向、汇流累积量、河网提取。但GDEM本身有噪声,直接算流向会出现很多虚假的洼地。我的做法是先做洼地填充,用GDAL的fillnodata或者用Python的richdem库。
import richdem as rd # 读取镶嵌后的DEM dem = rd.LoadGDAL('final.tif') # 洼地填充 dem_filled = rd.FillDepressions(dem, epsilon=True, in_place=False) # 计算流向 flow_dir = rd.FlowDirD8(dem_filled, in_place=False) # 计算汇流累积量 flow_acc = rd.FlowAccumulation(flow_dir, method='D8') # 保存结果 rd.SaveGDAL('flow_acc.tif', flow_acc)richdem的FillDepressions函数可以自动填充洼地,epsilon=True表示用微小坡度填充,避免平坦区域出现流向不确定的问题。汇流累积量可以用来提取河网,设定一个阈值,比如1000个像素,超过这个阈值的区域就是河道。
6.3 三维可视化与场景搭建
镶嵌后的GDEM可以导入Blender或者Unity做三维地形可视化。我一般把GDEM转成PNG的高度图,然后在Blender里用Displace修改器生成地形。具体做法是:先用gdal_translate把GDEM转成16位PNG,然后在Blender里创建一个平面,细分足够多的网格,用Displace修改器加载PNG作为高度图。
提示:转PNG的时候要注意拉伸对比度。GDEM的高程范围可能很大,直接转PNG会丢失细节。我一般先用gdal_calc做线性拉伸,把高程范围映射到0到65535之间,再转PNG。
6.4 后续可以扩展的方向
这个镶嵌流程还可以扩展到其他GDEM数据源,比如SRTM、ALOS World 3D。不同数据源的镶嵌逻辑类似,但接边处理的参数需要调整。SRTM的分辨率是90米,ALOS World 3D是30米,精度和覆盖范围也不一样。如果你需要更高精度的数据,可以考虑用商业数据源,但那个成本就上去了。
另一个扩展方向是自动化。我现在是把整个流程拆成几个脚本,手动按顺序跑。如果你需要频繁做镶嵌,可以把这些脚本串起来,加一个配置文件,指定输入目录、输出目录、投影参数、分辨率、裁剪范围,一键跑完。我用Python的argparse写了一个简单的命令行工具,支持这些参数,用起来方便不少。
还有一个方向是质量评估的自动化。我现在质检是靠人工看,效率不高。后来我写了一个脚本,自动计算镶嵌后数据的统计信息,包括均值、标准差、无效值比例、接缝处的高程梯度。如果某个指标超过阈值,就报警提示。这个脚本帮我省了不少时间,特别是做大批量镶嵌的时候。
我个人在实际操作中的体会是,ASTER GDEM镶嵌这件事,难点不在技术本身,而在细节的把控。投影要统一,分辨率要一致,接边要调整,无效值要填充,边缘要裁剪,每一步都不能省。省了一步,后面就要花更多时间返工。另外,中间文件一定要保留,特别是投影转换后的文件,因为那一步最耗时,万一后面发现问题,不用从头再来。最后再分享一个小技巧:如果你不确定接边调整的阈值设多少合适,可以先设一个保守的值,比如3米,跑一遍看看效果,如果接缝还是明显,再逐步放宽到5米、10米。不要一上来就设很大的阈值,那样可能会把真实的地形差异也抹掉。