别再到处找DEM了！手把手教你用ArcGIS Pro + Python脚本，从NASA官网免费下载并拼接出完整的中国90米高程数据

从NASA获取中国90米高程数据的自动化解决方案在GIS和遥感研究领域获取高质量的数字高程模型(DEM)数据是许多项目的基础工作。然而对于中国区域的完整覆盖、高精度且免费可用的DEM数据研究者们常常面临获取困难。本文将介绍如何利用ArcGIS Pro结合Python脚本直接从NASA官网获取SRTM 90米分辨率数据并实现自动化处理和拼接最终生成完整的中国区域DEM数据集。1. 数据源选择与准备NASA提供的SRTM(Shuttle Radar Topography Mission)数据是目前最广泛使用的全球DEM数据之一。90米分辨率的SRTM数据覆盖全球且完全免费开放获取。相比其他来源NASA的数据具有以下优势权威性由NASA官方发布数据质量有保障全球覆盖完整覆盖中国全境免费获取无需付费或特殊权限持续更新数据版本会定期修正和优化要获取这些数据我们首先需要了解NASA Earthdata的访问机制。NASA要求用户在下载数据前先注册账号并获取API密钥。以下是注册步骤访问Earthdata登录页面点击Register创建新账户填写必要信息并完成验证登录后在个人资料页面可以找到API访问凭证提示建议将API凭证保存在安全位置但不要直接硬编码在脚本中考虑使用环境变量或配置文件管理。2. 自动化下载流程设计手动下载SRTM数据不仅耗时而且容易出错。我们将使用Python脚本实现自动化下载流程。核心思路是确定中国区域的经纬度范围根据范围计算需要下载的SRTM图幅编号构建NASA数据服务器的请求URL批量下载所有需要的图幅以下是实现这一流程的关键Python代码片段import requests import os from concurrent.futures import ThreadPoolExecutor # NASA Earthdata认证信息 username os.getenv(EARTHDATA_USER) password os.getenv(EARTHDATA_PASS) # 中国区域边界坐标 china_bbox { min_lon: 73.5, max_lon: 135.0, min_lat: 18.0, max_lat: 53.5 } def get_srtm_tiles(bbox): 根据边界框计算需要的SRTM图幅编号 tiles [] # 计算逻辑略... return tiles def download_srtm_tile(tile): 下载单个SRTM图幅 url fhttps://e4ftl01.cr.usgs.gov/MEASURES/SRTMGL1.003/2000.02.11/{tile}.SRTMGL1.hgt.zip response requests.get(url, auth(username, password), streamTrue) # 保存处理略...为了提高下载效率我们使用多线程并发下载。同时代码中加入了断点续传和错误重试机制确保在大批量下载时的稳定性。3. 数据处理与格式转换下载的原始SRTM数据是HGT格式需要转换为ArcGIS Pro可以直接处理的格式。此外原始数据还存在以下常见问题需要处理数据缺失部分区域可能有数据空洞投影问题原始数据使用地理坐标系(WGS84)边缘效应图幅边缘可能出现不连续我们使用arcpy模块实现自动化处理流程import arcpy from arcpy.sa import * # 设置工作空间 arcpy.env.workspace path/to/workspace arcpy.env.overwriteOutput True def process_srtm_tile(input_hgt, output_raster): 处理单个SRTM图幅 # 转换为ArcGIS格式 arcpy.management.CopyRaster(input_hgt, output_raster) # 修复数据空洞 filled Con(IsNull(output_raster), FocalStatistics(output_raster), output_raster) # 重投影到适合中国区域的投影坐标系 arcpy.management.ProjectRaster(filled, projected_ output_raster, PROJCS[Asia_North_Albers_Equal_Area_Conic]) return projected_ output_raster对于数据缺失问题我们采用邻域填充算法进行修复。投影转换则选择适合中国区域的Albers等面积圆锥投影以减少面积计算时的变形。4. 数据拼接与裁剪当所有图幅都完成下载和处理后需要将它们拼接成一个完整的中国区域DEM数据集。这一过程需要考虑边缘匹配确保相邻图幅间无缝衔接数据融合消除图幅间的亮度差异范围裁剪精确匹配中国国界以下是实现拼接和裁剪的关键步骤创建镶嵌数据集# 创建镶嵌数据集 arcpy.management.CreateMosaicDataset(China_DEM.gdb, China_SRTM, PROJCS[Asia_North_Albers_Equal_Area_Conic]) # 添加所有处理后的图幅 raster_list [tile1.tif, tile2.tif, ...] # 所有图幅列表 arcpy.management.AddRastersToMosaicDataset(China_SRTM, Raster Dataset, raster_list)执行镶嵌操作# 设置镶嵌参数 arcpy.management.SetMosaicDatasetProperties(China_SRTM, mosaic_methodBLEND, blend_width2) # 生成最终镶嵌结果 arcpy.management.BuildOverview(China_SRTM)按中国边界裁剪# 加载中国边界矢量数据 china_boundary path/to/china_boundary.shp # 执行裁剪 arcpy.management.Clip(China_SRTM, in_memory/clipped, china_boundary, , ClippingGeometry)注意中国边界数据应使用官方发布的权威数据源确保边界精确性和合规性。5. 质量控制与优化生成最终DEM后需要进行质量检查以确保数据可用性。主要检查项包括检查项目方法标准数据完整性统计空值比例空值比例0.1%高程范围计算统计值符合中国地形特征边缘连续性可视化检查无明显接缝投影正确性元数据验证与设计一致对于发现的问题可以采取以下优化措施高程异常修正使用邻域统计值替换明显异常值平滑处理应用焦点统计减少数据噪声分辨率优化必要时进行重采样def quality_control(input_dem): DEM质量检查与优化 # 计算空值比例 null_count arcpy.management.GetRasterProperties(input_dem, COUNT)[0] total_cells arcpy.management.GetRasterProperties(input_dem, MEAN)[1] null_ratio float(null_count) / float(total_cells) if null_ratio 0.001: # 执行空值填充 filled Con(IsNull(input_dem), FocalStatistics(input_dem), input_dem) arcpy.management.CopyRaster(filled, input_dem) # 其他优化处理...6. 成果输出与应用完成所有处理后我们可以将最终DEM数据输出为多种格式以适应不同应用场景GeoTIFF通用栅格格式兼容大多数GIS软件File GeodatabaseArcGIS原生格式支持高效存储Cloud Optimized GeoTIFF (COG)适合网络共享和在线应用对于大型项目建议将处理流程封装为ArcGIS Pro工具箱工具方便团队其他成员使用。以下是一个简单的工具封装示例import arcpy class DemProcessor(object): def __init__(self): self.label SRTM DEM Processor self.description Automated SRTM DEM processing tool def getParameterInfo(self): params [] # 定义输入参数... return params def execute(self, parameters, messages): # 执行处理逻辑... arcpy.AddMessage(Processing completed successfully)在实际项目中这套自动化流程可以节省大量时间。以处理整个中国区域90米SRTM数据为例手动处理约需3-5天且容易出错自动化流程首次运行约4-6小时后续更新仅需1-2小时7. 进阶技巧与问题排查即使采用自动化流程在实际操作中仍可能遇到各种问题。以下是一些常见问题及解决方案问题1NASA服务器连接不稳定解决方案使用会话保持和重试机制避开美国白天的高峰时段考虑使用NASA提供的wget脚本批量下载import time from requests.adapters import HTTPAdapter from requests.packages.urllib3.util.retry import Retry def requests_retry_session(retries3, backoff_factor0.3): session requests.Session() retry Retry( totalretries, readretries, connectretries, backoff_factorbackoff_factor, status_forcelist(500, 502, 504), ) adapter HTTPAdapter(max_retriesretry) session.mount(http://, adapter) session.mount(https://, adapter) return session问题2图幅间存在高程差异解决方案使用直方图匹配算法统一色调在镶嵌时设置适当的重叠区域和混合宽度必要时进行手动调整问题3投影转换导致变形解决方案选择适合中国区域的投影如Albers等面积圆锥投影控制投影参数优化变形分布对大区域数据进行分块处理在实际使用中我发现将整个处理流程分解为多个独立步骤并设置检查点非常有用。这样当某一步骤失败时可以从最近的检查点重新开始而不必从头运行整个流程。