Python GDAL实战ASTER GDEM V003全国DEM高程数据处理1123幅分幅合并→34省裁切30m精度摘要本文详细记录了基于 Pythonrasterio GDAL对 ASTER GDEM V003 原始数据进行全国34省级行政区DEM裁切的完整工程实践。涵盖数据下载、分幅合并mosaic、矢量裁切mask、质量检验全流程附完整可复现代码。数据覆盖N18°N53°、E73°E134°成果34个GeoTIFF文件总计10.61GB124.8亿有效像素。关键词DEM高程数据、ASTER GDEM、Python rasterio、GDAL、遥感数据处理、GeoTIFF、空间数据裁切、数字高程模型标签GIS遥感空间数据库数据处理地理信息系统成果免费获取方式数据由WX号YouGIS顽石整理分享完全免费。两种获取方式任选其一方式一关键词推荐关注YouGIS顽石发送DEM-省份简称 获取省份成果数据DEM-分幅号 获取分幅数据DEM-{省份简称} → 获取省份裁切成果数据如 DEM-河南 DEM-{分幅号} → 获取原始分幅数据如 DEM-N30E14方式二YouGIS 数据助手平台入口小程序搜索「YouGIS数据助手」PC 端https://yougis.com.cn/res/home支持在线查看数据元信息分辨率、坐标系、覆盖范围适合批量检索多个数据集。一、数据背景与选型1.1 ASTER GDEM V003 vs V002 vs SRTM对比项SRTM C-bandASTER GDEM V002ASTER GDEM V003分辨率30m30m30m纬度覆盖N60°~S56°N83°~S83°N83°~S83°垂直精度(RMSE)~16m17m8.5m数据源年份2000年2月2000-20102000-2018空洞密度少较多大幅减少中国高纬度覆盖差漠河N53°已覆盖但更北缺失好好选型结论全国DEM处理选 ASTER GDEM V003兼顾覆盖范围含高纬度和精度。SRTM适合低纬度区域专项研究。1.2 数据获取来源来源网址特点NASA EarthDatahttps://earthdata.nasa.gov原始来源需注册地理空间数据云https://www.gscloud.cn国内镜像下载速度快ASTER GDEM官网https://asterweb.jpl.nasa.gov/gdem.asp产品说明文档二、原始数据技术规格2.1 核心参数数据产品ASTER GDEM V003 发布机构NASA METI 发布时间2019年8月 分幅规则1°×1° 经纬度网格 单幅像素3601 × 3601 单幅大小~41 MB 空间分辨率0.000278°赤道约30m 数据类型signed 16-bit integer (int16) 空间参考WGS84 / EPSG:4326 NoData值-9999 文件格式GeoTIFF 命名规则ASTGTMV003_N{纬度}E{经度}_dem.tif2.2 全国覆盖统计分幅总数1,123幅 覆盖范围N18°~N53°, E73°~E134° 总数据量~49 GB 单幅示例ASTGTMV003_N39E116_dem.tif → 北京 ASTGTMV003_N30E114_dem.tif → 武汉2.3 分幅筛选——根据省份bbox确定所需图幅importnumpyasnpdefget_tile_range(min_lat,max_lat,min_lon,max_lon):根据经纬度范围获取需要的ASTER GDEM分幅编号latsrange(int(np.floor(min_lat)),int(np.ceil(max_lat))1)lonsrange(int(np.floor(min_lon)),int(np.ceil(max_lon))1)tiles[fN{lat}E{lon}forlatinlatsforloninlons]returntiles# 示例河南省范围约 N31°~N36°, E110°~E117°tilesget_tile_range(31,36,110,117)print(f河南需{len(tiles)}幅:{tiles})# 输出: 河南需 42 幅: [N31E110, N31E111, ..., N36E117]三、处理流程与核心代码┌──────────┐ ┌──────────┐ ┌──────────┐ ┌──────────┐ │ 数据下载 │ - │ 分幅合并 │ - │ 省界裁切 │ - │ 质量检验 │ │ 1123幅 │ │ 按省镶嵌 │ │ 矢量裁切 │ │ 统计核查 │ └──────────┘ └──────────┘ └──────────┘ └──────────┘ ~49 GB 临时大文件 34个TIF _summary.json3.1 分幅合并Mosaicimportrasteriofromrasterio.mergeimportmergefromrasterio.transformimportarray_boundsdefmosaic_tiles(tile_paths,output_path):合并多个DEM分幅为单个TIFFLZW压缩src_files[rasterio.open(p)forpintile_paths]mosaic,transformmerge(src_files)profilesrc_files[0].profile.copy()profile.update({height:mosaic.shape[1],width:mosaic.shape[2],transform:transform,compress:lzw,tiled:True,# 启用分块写入提升读取性能blockxsize:256,blockysize:256,})withrasterio.open(output_path,w,**profile)asdst:dst.write(mosaic)forsrcinsrc_files:src.close()print(f[OK] 合并完成:{output_path}({mosaic.shape[2]}×{mosaic.shape[1]}))3.2 省界裁切Maskimportrasteriofromrasterio.maskimportmaskimportjsondefclip_to_province(input_tif,geojson_path,output_tif,nodata-9999):用省级行政区划矢量GeoJSON裁切栅格withopen(geojson_path)asf:geojsonjson.load(f)# 处理 FeatureCollection / Feature / Geometry 等不同结构ifgeojson[type]FeatureCollection:shapes[feat[geometry]forfeatingeojson[features]]elifgeojson[type]Feature:shapes[geojson[geometry]]else:shapes[geojson]withrasterio.open(input_tif)assrc:out_image,out_transformmask(src,shapes,cropTrue,nodatanodata,all_touchedTrue)profilesrc.profile.copy()profile.update({height:out_image.shape[1],width:out_image.shape[2],transform:out_transform,nodata:nodata,compress:lzw,tiled:True,})withrasterio.open(output_tif,w,**profile)asdst:dst.write(out_image)print(f[OK] 裁切完成:{output_tif})3.3 质量检验importrasterioimportnumpyasnpfrompathlibimportPathimportjsondefvalidate_dem(tif_path):验证DEM数据质量返回统计字典withrasterio.open(str(tif_path))assrc:arrsrc.read(1)validarr[arr!src.nodata]return{code:Path(tif_path).stem,width:src.width,height:src.height,crs:str(src.crs),dtype:str(src.dtypes[0]),total_pixels:int(arr.size),valid_pixels:int(valid.size),valid_ratio:round(valid.size/arr.size*100,2),elev_min:int(valid.min()),elev_max:int(valid.max()),elev_mean:round(float(valid.mean()),1),size_mb:Path(tif_path).stat().st_size//(1024*1024),}# 批量检验results{}fortifinsorted(Path(a_output).glob(*.tif)):statsvalidate_dem(tif)results[stats[code]]statsprint(f{stats[code]}:{stats[valid_pixels]:12,}px | f{stats[elev_min]:5}~{stats[elev_max]:5}m | f{stats[size_mb]}MB)# 导出汇总JSONwithopen(_summary.json,w)asf:json.dump(results,f,ensure_asciiFalse,indent2)四、成果数据总览4.1 整体统计指标数值成果文件数34个总数据量10,868 MB (10.61 GB)总有效像素12,476,309,771约124.8亿有效像素占比43.42%高程范围-275m ~ 8802m数据类型int16坐标系WGS84 (EPSG:4326)NoData值-9999压缩方式LZW4.2 34省完整清单编码省份图幅数像素尺寸有效像素高程范围大小(MB)110000北京64008×28659,258,6263~2296m19120000天津43605×34918,632,021-5~1097m7130000河北3910845×947065,857,0020~2818m199140000山西279010×719544,388,292148~3055m195150000内蒙古46428848×19948280,552,49091~2840m1230210000辽宁3110846×900162,045,913-275~1337m162220000吉林4312642×893065,444,8910~2630m214230000黑龙江8418054×13394145,335,4350~1671m519310000上海53962×28896,837,004-10~352m4320000江苏309014×690038,332,587-32~685m61330000浙江257937×722634,425,8360~1920m118340000安徽279010×701942,437,091-2~1869m122350000福建238292×829240,449,7600~2160m150360000江西299369×795847,502,717-1~2145m177370000山东4010849×740745,702,055-11~1532m134410000河南4210489×852950,054,54912~2411m156420000湖北4012646×901462,375,665-6~3097m201430000湖南3510489×972856,001,68819~2097m235440000广东3210493×936655,021,5680~1901m185450000广西4812646×1009065,193,980-16~2109m271460000海南309370×630023,508,4200~1839m32500000重庆217585×684536,514,21226~2776m106510000四川7016258×14022125,472,821169~7523m658520000贵州309369×865054,395,708247~2896m223530000云南6014479×11586103,015,27776~6710m506540000西藏15128848×21612260,749,872122~8802m1476610000陕西4011711×901365,016,912166~3762m266620000甘肃18720704×16057172,464,233615~5786m507630000青海9719853×16625205,670,3091682~6751m801640000宁夏126127×541222,346,2511090~3554m56650000新疆22434260×23418357,088,649-154~8199m1833710000台湾306482×793819,908,7800~3886m44810000香港21622×18011,582,8400~958m1820000澳门11082×1441762,668-2~171m0五、工程注意事项5.1 大文件内存管理新疆成果文件 650000.tif 像素尺寸达 34260×234188亿像素直接全量读取会触发内存溢出。# ❌ 错误全量读取大文件withrasterio.open(650000.tif)assrc:arrsrc.read(1)# 1.6GB内存# ✅ 正确分块读取windowed readwithrasterio.open(650000.tif)assrc:forwindowinsrc.block_windows(1):blocksrc.read(1,windowwindow[1])# 逐块处理...5.2 坐标系一致性检查裁切前务必确认矢量与栅格坐标系一致否则裁切结果为空或错位。importrasterioimportgeopandasasgpdwithrasterio.open(merged.tif)assrc:raster_crssrc.crs.to_epsg()gdfgpd.read_file(province_boundary.geojson)vector_crsgdf.crs.to_epsg()ifraster_crs!vector_crs:print(f[WARNING] 坐标系不一致! 栅格: EPSG:{raster_crs}, 矢量: EPSG:{vector_crs})gdfgdf.to_crs(raster_crs)# 统一到栅格坐标系5.3 噪声像素过滤650000.tif新疆存在 104 个极端异常像素-32287~32008m使用前过滤importrasterioimportnumpyasnpwithrasterio.open(650000.tif)assrc:profilesrc.profile.copy()arrsrc.read(1).astype(np.float32)# 过滤物理不可能的高程值arr[(arr-1000)|(arr9000)]-9999withrasterio.open(650000_filtered.tif,w,**profile)asdst:dst.write(arr.astype(np.int16),1)5.4 投影转换注意事项数据为经纬度坐标EPSG:4326进行以下分析前需投影分析类型推荐投影原因面积/体积计算Albers等积投影保持面积不变距离/坡度分析UTM投影局部区域变形小流域提取UTM投影水流方向计算需平面坐标# 使用GDAL进行投影转换命令行# gdalwarp -t_srs EPSG:32649 input.tif output_utm49n.tif# Python方式importsubprocess subprocess.run([gdalwarp,-t_srs,EPSG:32649,-r,bilinear,-of,GTiff,-co,COMPRESSLZW,input.tif,output_utm49n.tif])六、数据读取与可视化importrasterioimportnumpyasnpimportmatplotlib.pyplotaspltwithrasterio.open(410000.tif)assrc:demsrc.read(1)demnp.where(dem-9999,np.nan,dem).astype(np.float32)fig,axplt.subplots(figsize(10,8))imax.imshow(dem,cmapterrain)plt.colorbar(im,label高程 (m),shrink0.8)ax.set_title(河南省 DEM (ASTER GDEM V003),fontsize14)ax.set_xlabel(列号)ax.set_ylabel(行号)plt.tight_layout()plt.savefig(henan_dem.png,dpi150)plt.show()七、总结阶段内容规模原始数据ASTER GDEM V0031123幅~49GB处理下载→合并→裁切→检验rasterio GDAL成果int16, EPSG:4326, LZW压缩34省TIF10.61GB处理流程全代码开源、可复现所有34个成果文件技术规格统一。参考文献NASA/METI. ASTER GDEM Version 3. (2019). https://asterweb.jpl.nasa.gov/gdem.aspNASA EarthData Search. https://earthdata.nasa.gov地理空间数据云. https://www.gscloud.cnrasterio Documentation. https://rasterio.readthedocs.ioGDAL Documentation. https://gdal.org如果觉得有用点赞收藏关注一键三连有问题评论区交流 ✌️