Python GDAL实战：ASTER GDEM V003全国DEM高程数据处理（1123幅分幅合并→34省裁切，30m精度）

摘要：本文详细记录了基于 Python（rasterio + GDAL）对 ASTER GDEM V003 原始数据进行全国34省级行政区DEM裁切的完整工程实践。涵盖数据下载、分幅合并（mosaic）、矢量裁切（mask）、质量检验全流程，附完整可复现代码。数据覆盖N18°~N53°、E73°~E134°，成果34个GeoTIFF文件，总计10.61GB，124.8亿有效像素。

关键词：DEM高程数据、ASTER GDEM、Python rasterio、GDAL、遥感数据处理、GeoTIFF、空间数据裁切、数字高程模型

标签：GIS 遥感 空间数据库 数据处理 地理信息系统

---

一、数据背景与选型

1.1 ASTER GDEM V003 vs V002 vs SRTM

选型结论：全国DEM处理选 ASTER GDEM V003，兼顾覆盖范围（含高纬度）和精度。SRTM适合低纬度区域专项研究。

1.2 数据获取来源

---

二、原始数据技术规格

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.tif

2.2 全国覆盖统计

分幅总数：1,123幅
覆盖范围：N18°~N53°, E73°~E134°
总数据量：~49 GB
单幅示例：ASTGTMV003_N39E116_dem.tif → 北京
         ASTGTMV003_N30E114_dem.tif → 武汉

2.3 分幅筛选——根据省份bbox确定所需图幅

import numpy as np

def get_tile_range(min_lat, max_lat, min_lon, max_lon):
    """根据经纬度范围获取需要的ASTER GDEM分幅编号"""
    lats = range(int(np.floor(min_lat)), int(np.ceil(max_lat)) + 1)
    lons = range(int(np.floor(min_lon)), int(np.ceil(max_lon)) + 1)
    tiles = [f"N{lat}E{lon}" for lat in lats for lon in lons]
    return tiles

# 示例：河南省范围约 N31°~N36°, E110°~E117°
tiles = get_tile_range(31, 36, 110, 117)
print(f"河南需 {len(tiles)} 幅: {tiles}")
# 输出: 河南需 42 幅: ['N31E110', 'N31E111', ..., 'N36E117']

---

三、处理流程与核心代码

┌──────────┐    ┌──────────┐    ┌──────────┐    ┌──────────┐
│  数据下载  │ -> │  分幅合并  │ -> │  省界裁切  │ -> │  质量检验  │
│ 1123幅   │    │ 按省镶嵌  │    │ 矢量裁切  │    │  统计+核查 │
└──────────┘    └──────────┘    └──────────┘    └──────────┘
   ~49 GB          临时大文件       34个TIF       _summary.json

3.1 分幅合并（Mosaic）

import rasterio
from rasterio.merge import merge
from rasterio.transform import array_bounds

def mosaic_tiles(tile_paths, output_path):
    """合并多个DEM分幅为单个TIFF（LZW压缩）"""
    src_files = [rasterio.open(p) for p in tile_paths]
    mosaic, transform = merge(src_files)

    profile = src_files[0].profile.copy()
    profile.update({
        'height': mosaic.shape[1],
        'width': mosaic.shape[2],
        'transform': transform,
        'compress': 'lzw',
        'tiled': True,        # 启用分块写入，提升读取性能
        'blockxsize': 256,
        'blockysize': 256,
    })

    with rasterio.open(output_path, 'w', **profile) as dst:
        dst.write(mosaic)

    for src in src_files:
        src.close()

    print(f"[OK] 合并完成: {output_path} ({mosaic.shape[2]}×{mosaic.shape[1]})")

3.2 省界裁切（Mask）

import rasterio
from rasterio.mask import mask
import json

def clip_to_province(input_tif, geojson_path, output_tif, nodata=-9999):
    """用省级行政区划矢量（GeoJSON）裁切栅格"""
    with open(geojson_path) as f:
        geojson = json.load(f)

    # 处理 FeatureCollection / Feature / Geometry 等不同结构
    if geojson['type'] == 'FeatureCollection':
        shapes = [feat['geometry'] for feat in geojson['features']]
    elif geojson['type'] == 'Feature':
        shapes = [geojson['geometry']]
    else:
        shapes = [geojson]

    with rasterio.open(input_tif) as src:
        out_image, out_transform = mask(
            src, shapes, crop=True, nodata=nodata, all_touched=True
        )
        profile = src.profile.copy()
        profile.update({
            'height': out_image.shape[1],
            'width': out_image.shape[2],
            'transform': out_transform,
            'nodata': nodata,
            'compress': 'lzw',
            'tiled': True,
        })

        with rasterio.open(output_tif, 'w', **profile) as dst:
            dst.write(out_image)

    print(f"[OK] 裁切完成: {output_tif}")

3.3 质量检验

import rasterio
import numpy as np
from pathlib import Path
import json

def validate_dem(tif_path):
    """验证DEM数据质量，返回统计字典"""
    with rasterio.open(str(tif_path)) as src:
        arr = src.read(1)
        valid = arr[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 = {}
for tif in sorted(Path('a_output').glob('*.tif')):
    stats = validate_dem(tif)
    results[stats['code']] = stats
    print(f"{stats['code']}: {stats['valid_pixels']:>12,} px | "
          f"{stats['elev_min']:>5}~{stats['elev_max']:<5}m | "
          f"{stats['size_mb']}MB")

# 导出汇总JSON
with open('_summary.json', 'w') as f:
    json.dump(results, f, ensure_ascii=False, indent=2)

---

四、成果数据总览

4.1 整体统计

4.2 34省完整清单

---

五、工程注意事项

5.1 大文件内存管理

新疆成果文件 650000.tif 像素尺寸达 34260×23418（8亿+像素），直接全量读取会触发内存溢出。

# ❌ 错误：全量读取大文件
with rasterio.open('650000.tif') as src:
    arr = src.read(1)  # 1.6GB内存！

# ✅ 正确：分块读取（windowed read）
with rasterio.open('650000.tif') as src:
    for window in src.block_windows(1):
        block = src.read(1, window=window[1])
        # 逐块处理...

5.2 坐标系一致性检查

裁切前务必确认矢量与栅格坐标系一致，否则裁切结果为空或错位。

import rasterio
import geopandas as gpd

with rasterio.open('merged.tif') as src:
    raster_crs = src.crs.to_epsg()

gdf = gpd.read_file('province_boundary.geojson')
vector_crs = gdf.crs.to_epsg()

if raster_crs != vector_crs:
    print(f"[WARNING] 坐标系不一致! 栅格: EPSG:{raster_crs}, 矢量: EPSG:{vector_crs}")
    gdf = gdf.to_crs(raster_crs)  # 统一到栅格坐标系

5.3 噪声像素过滤

650000.tif（新疆）存在 104 个极端异常像素（-32287~32008m），使用前过滤：

import rasterio
import numpy as np

with rasterio.open('650000.tif') as src:
    profile = src.profile.copy()
    arr = src.read(1).astype(np.float32)
    # 过滤物理不可能的高程值
    arr[(arr < -1000) | (arr > 9000)] = -9999

    with rasterio.open('650000_filtered.tif', 'w', **profile) as dst:
        dst.write(arr.astype(np.int16), 1)

5.4 投影转换注意事项

数据为经纬度坐标（EPSG:4326），进行以下分析前需投影：

# 使用GDAL进行投影转换（命令行）
# gdalwarp -t_srs EPSG:32649 input.tif output_utm49n.tif

# Python方式
import subprocess
subprocess.run([
    'gdalwarp', '-t_srs', 'EPSG:32649',
    '-r', 'bilinear', '-of', 'GTiff',
    '-co', 'COMPRESS=LZW',
    'input.tif', 'output_utm49n.tif'
])

---

六、数据读取与可视化

import rasterio
import numpy as np
import matplotlib.pyplot as plt

with rasterio.open('410000.tif') as src:
    dem = src.read(1)
    dem = np.where(dem == -9999, np.nan, dem).astype(np.float32)

    fig, ax = plt.subplots(figsize=(10, 8))
    im = ax.imshow(dem, cmap='terrain')
    plt.colorbar(im, label='高程 (m)', shrink=0.8)
    ax.set_title('河南省 DEM (ASTER GDEM V003)', fontsize=14)
    ax.set_xlabel('列号')
    ax.set_ylabel('行号')
    plt.tight_layout()
    plt.savefig('henan_dem.png', dpi=150)
    plt.show()

---

七、总结

处理流程全代码开源、可复现，所有34个成果文件技术规格统一。

---

参考文献

1. NASA/METI. ASTER GDEM Version 3. (2019). https://asterweb.jpl.nasa.gov/gdem.asp

2. NASA EarthData Search. https://earthdata.nasa.gov

3. 地理空间数据云. https://www.gscloud.cn

4. rasterio Documentation. https://rasterio.readthedocs.io

5. GDAL Documentation. https://gdal.org

---

数据获取方式见下一篇： 全国DEM高程数据免费下载指南（34省+1123幅原始）

如果觉得有用，点赞收藏关注一键三连，有问题评论区交流 ✌️