简介厦门市数字高程模型DEM是一份面向GIS从业者、规划人员及地理信息学习者的地形栅格数据包。该数据以ArcGIS Grid格式存储涵盖厦门市全域地表海拔信息可用于坡度坡向分析、流域划分、洪水风险评估及城市基建选址等场景使用者需具备一定ArcGIS或QGIS操作基础。压缩包共64个文件主要包括27个nit网格索引文件、27个dat属性数据文件以及adf、xml、log等辅助文件完整呈现了GRID栅格数据集的目录结构便于在ArcGIS中直接加载与二次处理。整包仅1.16MB轻量易下载。目前已有579人学习参考适合用于课堂实验、区域地形研究或规划项目初判。借助该DEM数据用户可快速提取厦门市任意区域的高程、坡度与地形起伏信息为环境分析、工程规划或灾害模拟提供基础数据支撑节省自行采集与预处理的时间。1. 「厦门市数字高程模型DEM」到底在解决什么问题把「厦门市数字高程模型DEM」这句话翻译成工程需求大多数人要的不是一份文件而是一条能跑通的数据链路从哪里拿到覆盖厦门全域的高程数据、用哪种分辨率、怎么裁剪到市界、怎么把坑填掉、最后落成坡度坡向或等高线产品。厦门的地形很特殊——本岛加岛外四区城区贴着海岸线北部和西部有低山丘陵高程从0米直接跳到200米以上这种破碎地形对DEM的处理要求比平原城市高得多。反直觉的一点是下载DEM从来不是瓶颈真正的瓶颈在投影、基准面和无效值。同一个厦门Xiamen大学附近的海平面0米、同安山区的负值、以及ALOS和SRTM之间一米多的系统性偏差会把未经预处理的RAW数据变成一张无法使用的「花屏」。这篇文章按「数据源选型 → 下载与预处理 → 裁剪拼接与DSM转DEM → 地形分析与校验」的顺序给出一套可以直接复现的厦门市DEM处理方案兼顾还在找入口的初学者和已经在调参数的老手。2. 厦门DEM数据源怎么选ALOS、SRTM还是ASTER2.1 三种全球DEM的分辨率与适用边界做厦门这种市级范围的DEM常见做法是先从覆盖全球或亚洲区域的开源数据源里选基座。真正在用的主要有三套ALOS World 3D简称AW3D30对外常叫ALOS 12.5米数据、SRTM航天飞机雷达地形测绘、ASTER GDEM。三套数据的分辨率、覆盖年份和高程基准都有差异不能简单说「越细越好」。数据源对外常见分辨率实际网格高程基准厦门适用性ALOS AW3D3030米/12.5米部分付费版本1角秒约30mEGM96大地水准面城区细节最均衡推荐默认SRTM v430米1角秒1角秒EGM96山体阴影质感好但厦门填海区数据旧ASTER GDEM v330米1角秒EGM96多云地区有异常坑点厦门夏季影像易出伪影Copernicus GLO-3030米1角秒大地水准面欧洲中心产的质量稳定可作交叉验证选择排序一般是这样优先ALOS AW3D30因为它在亚洲区域经过大量校正厦门周边的海岛、岸线轮廓比SRTM锐利。SRTM虽然经典但它的采集年代较早厦门港、翔安新城等填海造地路段存在明显的高程「断层」。ASTER的问题则是光学影像立体像对生成的DEM在云层覆盖区域会留下坑洞厦门春季多雾容易出现局部凹陷。2.2 为什么我建议直接下12.5米DEM而不是30米热词里频繁出现「12.5米dem下载」「alos 12.5米dem数据」这背后有一个实际需求厦门岛内的建筑、道路、岸堤在30米网格上会被抹平比如环岛路高架和地面道路在SRTM里几乎分不清而在12.5米数据里能看出结构轮廓。这里要澄清一个常见误解——12.5米数据通常来自ALOS PALSAR雷达数据重采样而成的DEM它的垂直精度RMSE在城区大约4-5米和30米版本相比提升有限但水平方向对地物边界的刻画好一档。提示如果你要算的是土方量、库容这类「体积」指标30米和12.5米的差异在厦门小范围地块上会放大到10%以上建议直接上12.5米。如果只是做全域坡度图或俯瞰可视化30米足够处理速度反而快一个量级。另一个不常被提及的点是文件体积。厦门市域面积约1700平方公里12.5米分辨率下栅格尺寸约1.1万个像素×1.4万个像素单张GeoTIFF约400-600MB而30米数据只有约100-150MB。内存不足的老机器处理这个量级会出现明显的卡顿但现代机器问题不大。考虑到厦门城中村和岛内老城区的地块尺度普遍很小12.5米能保住的边界细节是值得的。2.3 数据源官方入口和本地文件组织常见承载ALOS数据的入口是日本JAXA的AW3D平台以及OpenTopography。OpenTopography最大的好处是支持按多边形区域直接提交任务下载结果自带投影信息省去自己拼接的麻烦。下载之前建议先建立一套文件组织约定/workspace/xiamen_dem/ ├── source/ # 原始下载文件 │ ├── N24E117.tif │ ├── N24E118.tif │ └── N25E117.tif ├── boundary/ # 厦门市行政边界 │ └── xiamen_boundary.shp ├── cut/ # 裁剪后的中间产物 ├── final/ # 最终交付的DEM └── scripts/ # 处理脚本之所以必须按分幅文件管理是因为厦门的经纬度跨度大致在E117.5°-118.6°、N24.3°-24.8°之间不管用哪种数据源都不是单张影像能覆盖的。至少会横跨两到三个分幅而每幅的无效值范围和边缘重叠区域都不一致后期需要统一处理。3. 实操OpenTopography下载厦门12.5米DEM的最小步骤3.1 在OpenTopography上按边界选数据OpenTopography的常规流程是进入网站「Data」页签选择数据集类别里的「Global DEM」先在地图上框出厦门区域的AOIArea of Interest然后在数据源列表里勾选ALOS World 3D - 12.5m。提交后系统会用AWS的Lambda自动处理一般几分钟到十几分钟不等完成后会提供一个S3下载链接。按我的经验手工在网页上拖拽框选厦门全域有个隐蔽问题——很容易漏掉大小金门附近海域的边缘分幅以及厦门最东侧的翔安部分。更可靠的做法是先下载一份厦门市行政边界的GeoJSON或SHP文件然后在OpenTopography的AOI上传功能里直接提交这个边界文件。这样系统会针对多边形范围做瓦片拼接和裁剪下载结果就是「厦门市」形状的DEM而不是一个矩形范围。下载回来的文件通常叫merged_dem.tif之类但要注意它用的是WGS84经纬度坐标高程值单位是米。在厦门本地坐标系CGCS2000 / 3度分带中央经线118.5°E没有建立之前所有处理都先在这个通用坐标里做最后再重投影。3.2 命令行方式批量验证下载文件完整性如果下载的是分幅文件而不是合并结果拿到手后第一件事不是立刻处理而是先看范围、坐标系、无效值。用GDAL的命令行工具几个命令就能检查清楚# 查看栅格基本信息 gdalinfo N24E117.tif # 快速统计高程范围检测负值 gdalinfo -stats N24E117.tif # 输出为ASCII格式方便快速抽样 gdal_translate -of XYZ N24E117.tif sample.xyzgdalinfo输出的关键字段包括Size is 3601, 3601表示像素尺寸、Origin左上角坐标、Pixel Size网格间距ALOS是0.00027778度约30米12.5米版本则不同、TypeInt16高程值按整型存单位米。特别注意NoData Value字段——如果它是默认的-32768后面所有计算里必须显式指定它否则裁剪和叠加时会把海洋区域的高程拉进统计里。检查步骤里最容易忽略的一步是确认无效值的空间分布。海洋区域在ALOS原始数据里是有效的高程0米但OpenTopography输出的产品通常已经把海域设置为-32768或NaN因此厦门本岛周围的海面不会变成一马平川而是显示为无数据。这个特性在做淹没分析或海岸线提取时有用但在计算全区域高程均值时会造成严重偏移。3.3 用Python批量拼接厦门周边分幅瓦片如果不想依赖OpenTopography的在线拼接也需要在本地掌握批量拼接的脚本原因之一是万一最终下载结果只有单幅后续自行处理时用得到。厦门周边可能涉及N24E117、N24E118、N25E117、N25E118等分幅用rasterio.merge是最好的方式import rasterio from rasterio.merge import merge from rasterio.plot import show import glob # 列出所有ALOS分幅文件 tif_files glob.glob(/workspace/xiamen_dem/source/N*.tif) print(f共找到 {len(tif_files)} 个分幅文件) # 执行拼接注意输出路径和输出类型 mosaic, out_transform merge(tif_files, methodfirst) mosaic_meta rasterio.open(tif_files[0]).meta.copy() mosaic_meta.update({ driver: GTiff, height: mosaic.shape[1], width: mosaic.shape[2], transform: out_transform, dtype: int16 }) with rasterio.open(/workspace/xiamen_dem/source/xiamen_mosaic.tif, w, **mosaic_meta) as dest: dest.write(mosaic)这段代码的关键逻辑是merge默认把重叠区域的像素取第一个文件的数值methodfirst对于DEM来说重叠带两边是同一数据源的相邻瓦片取值差异通常在0.1米以内用first没问题。而mosaic_meta要从原始文件继承NoData和坐标信息不能直接新建一个空meta那样会丢坐标系。最后写入时dtype保持int16不要在拼接阶段就转Float32否则文件体积翻倍。4. 厦门DEM的裁剪、重投影与DSM到DEM的修正4.1 用行政边界裁剪DEM并处理边界锯齿拿到合并后的矩形DEM后要用厦门市行政边界把它裁剪出来。这里有个操作细节不能直接用gdalwarp -cutline默认参数因为默认情况下裁剪边界上的像素会做最邻近重采样导致边界线上出现「锯齿状」台阶。厦门海岸线和行政区界有大量曲线比如翔安海岸线需要指定-crop_to_cutline和合适的重采样算法gdalwarp -cutline /workspace/xiamen_dem/boundary/xiamen_boundary.shp \ -crop_to_cutline \ -tr 0.0002777777777777778 0.0002777777777777778 \ -r bilinear \ -of GTiff \ /workspace/xiamen_dem/source/xiamen_mosaic.tif \ /workspace/xiamen_dem/cut/xiamen_cut_30m.tif-tr的数值是目标分辨率的经纬度跨度0.00027777度约等于30米如果你用的是12.5米数据应该改成约0.00011574度1/86400度。用bilinear双线性插值配合crop_to_cutline边界上的高程过渡会比较平滑。gdalwarp后面可以加上--config GDALWARP_IGNORE_BAD_CUTLINE YES这样边界SHP即使有小瑕疵也不会中断整条命令但对厦门这种边界数据质量较高的场景不是必须的。裁剪完成后要验证两个东西一是裁剪后的像素行列数是否需要在后续处理中记录二是裁剪边缘的高程值有没有出现「一圈0米」的异常。厦门市界不是完全贴合海岸线边界外的海域应该保留NoData而不是0这在后续做表面积计算时才不会把海域算进去。4.2 重投影到厦门本地坐标系目标CRS怎么定DEM处理可以全程在WGS84经纬度里做但一旦涉及和已有CAD数据叠加比如热词里提到的「dem dxf 叠加 python」坐标统一就是绕不开的。厦门常见的本地坐标系有厦门94地方坐标系以及基于CGCS2000的3度带高斯投影中央经线118.5°E。在没有明确标注的情况下后台GIS数据通常采用CGCS2000 / 3-degree Gauss-Kruger zone 35或UTM zone 50N。用GDAL做重投影# 从WGS84重投影到CGCS2000 / 3-degree Gauss-Kruger CM 118.5E gdalwarp -t_srs EPSG:4549 \ -r cubic \ -ot Float32 \ /workspace/xiamen_dem/cut/xiamen_cut_30m.tif \ /workspace/xiamen_dem/final/xiamen_dem_cgcs2000.tif这里用EPSG:4549CGCS2000 / 3-degree Gauss-Kruger CM 118.5E是因为厦门正好落在118.5度中央经线上变形最小。重采样选cubic三次卷积因为经过投影转换后网格方向变了最邻近会产生轻微的马赛克感而双线性会对地形细节有一定柔化。如果最终产品只做高程统计不做视觉分析回到bilinear也可以速度更快。重投影后注意检查NoData值有没有被保留有些GDAL版本在投影变换后会把NoData设置为0需要重新gdal_translate -a_nodata -32768修正。4.3 DSM生成DEM的差值法厦门高密度城区的特殊处理厦门地形分析的另一个常见需求是「dsm生成dem」——把包含建筑和植被的表面模型DSM滤出真实地面DEM。ALOS AW3D30这类光学雷达数据本身已经是DSM能看出城市建筑的轮廓是要付出代价的在厦门岛内老城区和软件园三期这种高密度区域DSM的地表高程比真实地面高几米到十几米直接拿它做坡度分析或洪水淹没会得到错误结论。最常用的DSM到DEM思路是「最小值滤波 局部差值」import numpy as np from scipy.ndimage import minimum_filter import rasterio with rasterio.open(/workspace/xiamen_dem/final/xiamen_dem_cgcs2000.tif) as src: dem src.read(1).astype(float32) profile src.profile # 用半径约50米的窗口取最小值抑制建筑和树冠 radius_px max(1, int(50 / 30)) # 30米分辨率下约2像素 ground minimum_filter(dem, sizeradius_px * 2 1, modenearest) # 只在地表高出地面超过3米的像素上做替换 mask (dem - ground) 3 dem_corrected np.where(mask, ground, dem) # 写回文件 profile.update(dtypefloat32, nodata-32768.0) with rasterio.open(/workspace/xiamen_dem/final/xiamen_dem_ground.tif, w, **profile) as dst: dst.write(dem_corrected, 1)这个方法的物理逻辑是在一个窗口内建筑和树冠的顶部不是地形取最小值能找出裸地但厦门大量农田和裸地区域本身是平坦的直接替换会影响真实地形所以只有差值超过3米才修正。厦门岛内建筑普遍在3层以上3米阈值刚好把低矮平房和雨棚排除在外。这个步骤在30米分辨率上效果一般因为30米网格里建筑和道路混在一起最小值滤波容易把道路挖成沟而在12.5米数据上效果会好得多建议有条件的时候优先对12.5米做这个处理。5. 厦门地形分析落地坡度、坡向与等高线的参数设置5.1 在厦门这种起伏地形下算坡度的参数选择完成DEM数据准备后最常见的分析任务是坡度、坡向、等高线。在厦门这些分析结果直接用在环评、选址、水土保持上比如同安的茶山地块坡度大于25度通常会被划为限制建设区。用GDAL直接生成坡度的命令是通用的但参数设置关乎结果正确性# 生成坡度输出单位为度 gdaldem slope /workspace/xiamen_dem/final/xiamen_dem_cgcs2000.tif \ /workspace/xiamen_dem/final/xiamen_slope.tif \ -p \ -s 1.0 # 生成坡向北为0度顺时针 gdaldem aspect /workspace/xiamen_dem/final/xiamen_dem_cgcs2000.tif \ /workspace/xiamen_dem/final/xiamen_aspect.tif-p参数让坡度以度为单位输出如果不加则默认是百分数tan值乘以100在厦门这种有山有平地的场景下百分数不利于直接设定阈值。-s是垂直高程单位与水平单位的比值因为数据已经重投影到米制坐标水平单位是米垂直单位也是米所以设1.0。一个常见错误是拿WGS84经纬度的DEM直接算坡度且忘记加-s这时GDAL默认水平单位是度算出来的坡度完全失真厦门山体动辄出现80度以上的假值。注意如果你跳过重投影直接在WGS84坐标系下运行gdaldem必须设置-s 111320.01度对应的米数否则坡度结果没有参考价值。5.2 等高线生成的关键参数间距与平滑等高线在厦门常用于精细的场地设计常见做法是用gdal_contour从DEM里抽线gdal_contour -a ELEV \ -i 5.0 \ /workspace/xiamen_dem/final/xiamen_dem_cgcs2000.tif \ /workspace/xiamen_dem/final/xiamen_contour_5m.shp-a ELEV指定属性字段名后面接SHP时这个属性就代表每根等高线的高程值。-i 5.0是等高距厦门岛内地形平缓5米比较合适同安山区高差大如果做的是全区域检索图建议用10米否则SHP文件过大且线条密集难读。厦门填海区东渡、高崎、翔安南部的高程通常在2-5米之间5米等高距会在这类区域产生非常少的线但实际上这些区域几乎完全平坦线少是正常的。gdal_contour生成的线默认是折线不做平滑。如果交付给设计院通常还需要在QGIS里对线做一次平滑或者用v.generalizeGRASS做Simplification。但平滑操作只影响视觉效果不影响高程属性所以精度敏感的工程验收不建议平滑。5.3 常见误用用DSM代替DEM做坡度分析的后果在厦门这类高密度城市直接用原始DSM做坡度分析会得到一个「非常难看」的结果岛内建筑密集区坡度普遍在10度以上看起来整个厦门老城区都建在山坡上。这显然是错的。用第4.3节修正后的地面DEM重新算一遍坡度同样的区域坡度会降到3度以下才符合真实地形感受。这个问题的根源在于DSM本身是「地表覆盖物顶部」的模型建筑屋顶、树木冠层都是它的组成部分。厦门绿化覆盖率高仙岳路沿线、狐尾山周边的树冠在高程模型上是连续起伏的「假山丘」直接分析会把真正的地形沟谷和排水路径全部掩盖。因此任何以「工程建设」为目的的分析都必须先做DSM到DEM的差值修正再出坡度和等高线。可视化或飞行漫游则不必做这个修正保留建筑轮廓反而信息量更大。6. 最后收口高程基准校验与DXF叠加的3个实用技巧6.1 用一个值得信赖的反例校验厦门附近GNSS水准点核对拿到处理完的厦门DEM我一般先做一次「单点校验」。厦门的实测高程控制点分布在岛内外不同区域从这些点里挑三到四个写入一个CSV包含经纬度和实测海拔然后用gdal_query.py或简单的rasterio采样脚本把DEM对应位置的高程取出来对比import csv import rasterio with rasterio.open(/workspace/xiamen_dem/final/xiamen_dem_ground.tif) as src: for row in csv.DictReader(open(control_points.csv)): lon, lat float(row[lon]), float(row[lat]) # 反算像素坐标 py, px src.index(lon, lat) dem_val src.read(1)[py, px] print(f{row[name]}: 实测 {row[h]}mDEM {dem_val}m差 {float(row[h])-dem_val:.2f}m)如果差值平均在±2米以内说明数据的垂直质量合格如果系统性偏向一侧比如所有点都比实测高2米基本可以判断是基准面差异而不是数据处理错误。厦门同安老城区的控制点数据常有基准不统一的情况碰到系统性误差时不要急着改数据先确认对照点使用的高程基准是否和DEM的EGM96一致。6.2 将DEM与DXF叠加的坐标系匹配技巧「dem dxf 叠加 python」是很多市政项目的刚需把CAD设计图通常带本地坐标叠加到DEM渲染的底图上做方案比对。这一步最常见的问题是两者相差几米到几十米不等对齐不上。此时用Python比在QGIS里手动拖拽来得可靠import geopandas as gpd import rasterio from rasterio.plot import show import matplotlib.pyplot as plt # 读取DXF指定EPSG:4549厦门常用CRS gdf gpd.read_file(/workspace/xiamen_dem/dwg/design.dxf) gdf gdf.to_crs(EPSG:4549) # 读取已重投影的DEM with rasterio.open(/workspace/xiamen_dem/final/xiamen_dem_cgcs2000.tif) as src: fig, ax plt.subplots(figsize(10, 10)) show(src, axax, cmapterrain) gdf.plot(axax, colorred, linewidth0.8) plt.savefig(/workspace/xiamen_dem/final/overlay_check.png, dpi150)如果dxf里的坐标不是米制而是毫米或者英尺转换时要在to_crs前用gdf.scale(xfact0.001, yfact0.001, origin(0,0))缩放。厦门本地院出图习惯不同这步检查必不可少。叠加后除了看套合度还要看DEM渲染的透明度和色带是否符合判别需求这些由参数补足。6.3 批量导出特定区域小范围地形切片最后送一个常用技巧从整市DEM里切出一小块矩形区域转成DXF或XYZ格式给无人机或其他专业软件用。gdal_translate配合-projwin参数可以做按坐标范围切片-of XYZ输出高度精简的三列文本同时这也是最容易读出来的格式gdal_translate -projwin 13425070 2751350 13428000 2748000 \ -of XYZ /workspace/xiamen_dem/final/xiamen_dem_ground.tif \ /workspace/xiamen_dem/final/slice_xiamen_island.xyz-projwin后跟的是左上角x、左上角y、右下角x、右下角y注意是「上」在前「下」在后很多人写成左右顺序结果切出反区域。切片完成后用head或wc -l快速检查点数确保没有把海洋区域的大面积NoData带进来。这样一整套流程从下载、预处理、修正到分析、校验就闭环了厦门的DEM数据在任何需要高程信息的GIS工程里都能直接成为可信的底图而不是一堆难以收拾的原始tif。本文还有配套的精品资源点击获取 厦门市DEM处理全流程：从数据源选型到地形分析实战
简介厦门市数字高程模型DEM是一份面向GIS从业者、规划人员及地理信息学习者的地形栅格数据包。该数据以ArcGIS Grid格式存储涵盖厦门市全域地表海拔信息可用于坡度坡向分析、流域划分、洪水风险评估及城市基建选址等场景使用者需具备一定ArcGIS或QGIS操作基础。压缩包共64个文件主要包括27个nit网格索引文件、27个dat属性数据文件以及adf、xml、log等辅助文件完整呈现了GRID栅格数据集的目录结构便于在ArcGIS中直接加载与二次处理。整包仅1.16MB轻量易下载。目前已有579人学习参考适合用于课堂实验、区域地形研究或规划项目初判。借助该DEM数据用户可快速提取厦门市任意区域的高程、坡度与地形起伏信息为环境分析、工程规划或灾害模拟提供基础数据支撑节省自行采集与预处理的时间。1. 「厦门市数字高程模型DEM」到底在解决什么问题把「厦门市数字高程模型DEM」这句话翻译成工程需求大多数人要的不是一份文件而是一条能跑通的数据链路从哪里拿到覆盖厦门全域的高程数据、用哪种分辨率、怎么裁剪到市界、怎么把坑填掉、最后落成坡度坡向或等高线产品。厦门的地形很特殊——本岛加岛外四区城区贴着海岸线北部和西部有低山丘陵高程从0米直接跳到200米以上这种破碎地形对DEM的处理要求比平原城市高得多。反直觉的一点是下载DEM从来不是瓶颈真正的瓶颈在投影、基准面和无效值。同一个厦门Xiamen大学附近的海平面0米、同安山区的负值、以及ALOS和SRTM之间一米多的系统性偏差会把未经预处理的RAW数据变成一张无法使用的「花屏」。这篇文章按「数据源选型 → 下载与预处理 → 裁剪拼接与DSM转DEM → 地形分析与校验」的顺序给出一套可以直接复现的厦门市DEM处理方案兼顾还在找入口的初学者和已经在调参数的老手。2. 厦门DEM数据源怎么选ALOS、SRTM还是ASTER2.1 三种全球DEM的分辨率与适用边界做厦门这种市级范围的DEM常见做法是先从覆盖全球或亚洲区域的开源数据源里选基座。真正在用的主要有三套ALOS World 3D简称AW3D30对外常叫ALOS 12.5米数据、SRTM航天飞机雷达地形测绘、ASTER GDEM。三套数据的分辨率、覆盖年份和高程基准都有差异不能简单说「越细越好」。数据源对外常见分辨率实际网格高程基准厦门适用性ALOS AW3D3030米/12.5米部分付费版本1角秒约30mEGM96大地水准面城区细节最均衡推荐默认SRTM v430米1角秒1角秒EGM96山体阴影质感好但厦门填海区数据旧ASTER GDEM v330米1角秒EGM96多云地区有异常坑点厦门夏季影像易出伪影Copernicus GLO-3030米1角秒大地水准面欧洲中心产的质量稳定可作交叉验证选择排序一般是这样优先ALOS AW3D30因为它在亚洲区域经过大量校正厦门周边的海岛、岸线轮廓比SRTM锐利。SRTM虽然经典但它的采集年代较早厦门港、翔安新城等填海造地路段存在明显的高程「断层」。ASTER的问题则是光学影像立体像对生成的DEM在云层覆盖区域会留下坑洞厦门春季多雾容易出现局部凹陷。2.2 为什么我建议直接下12.5米DEM而不是30米热词里频繁出现「12.5米dem下载」「alos 12.5米dem数据」这背后有一个实际需求厦门岛内的建筑、道路、岸堤在30米网格上会被抹平比如环岛路高架和地面道路在SRTM里几乎分不清而在12.5米数据里能看出结构轮廓。这里要澄清一个常见误解——12.5米数据通常来自ALOS PALSAR雷达数据重采样而成的DEM它的垂直精度RMSE在城区大约4-5米和30米版本相比提升有限但水平方向对地物边界的刻画好一档。提示如果你要算的是土方量、库容这类「体积」指标30米和12.5米的差异在厦门小范围地块上会放大到10%以上建议直接上12.5米。如果只是做全域坡度图或俯瞰可视化30米足够处理速度反而快一个量级。另一个不常被提及的点是文件体积。厦门市域面积约1700平方公里12.5米分辨率下栅格尺寸约1.1万个像素×1.4万个像素单张GeoTIFF约400-600MB而30米数据只有约100-150MB。内存不足的老机器处理这个量级会出现明显的卡顿但现代机器问题不大。考虑到厦门城中村和岛内老城区的地块尺度普遍很小12.5米能保住的边界细节是值得的。2.3 数据源官方入口和本地文件组织常见承载ALOS数据的入口是日本JAXA的AW3D平台以及OpenTopography。OpenTopography最大的好处是支持按多边形区域直接提交任务下载结果自带投影信息省去自己拼接的麻烦。下载之前建议先建立一套文件组织约定/workspace/xiamen_dem/ ├── source/ # 原始下载文件 │ ├── N24E117.tif │ ├── N24E118.tif │ └── N25E117.tif ├── boundary/ # 厦门市行政边界 │ └── xiamen_boundary.shp ├── cut/ # 裁剪后的中间产物 ├── final/ # 最终交付的DEM └── scripts/ # 处理脚本之所以必须按分幅文件管理是因为厦门的经纬度跨度大致在E117.5°-118.6°、N24.3°-24.8°之间不管用哪种数据源都不是单张影像能覆盖的。至少会横跨两到三个分幅而每幅的无效值范围和边缘重叠区域都不一致后期需要统一处理。3. 实操OpenTopography下载厦门12.5米DEM的最小步骤3.1 在OpenTopography上按边界选数据OpenTopography的常规流程是进入网站「Data」页签选择数据集类别里的「Global DEM」先在地图上框出厦门区域的AOIArea of Interest然后在数据源列表里勾选ALOS World 3D - 12.5m。提交后系统会用AWS的Lambda自动处理一般几分钟到十几分钟不等完成后会提供一个S3下载链接。按我的经验手工在网页上拖拽框选厦门全域有个隐蔽问题——很容易漏掉大小金门附近海域的边缘分幅以及厦门最东侧的翔安部分。更可靠的做法是先下载一份厦门市行政边界的GeoJSON或SHP文件然后在OpenTopography的AOI上传功能里直接提交这个边界文件。这样系统会针对多边形范围做瓦片拼接和裁剪下载结果就是「厦门市」形状的DEM而不是一个矩形范围。下载回来的文件通常叫merged_dem.tif之类但要注意它用的是WGS84经纬度坐标高程值单位是米。在厦门本地坐标系CGCS2000 / 3度分带中央经线118.5°E没有建立之前所有处理都先在这个通用坐标里做最后再重投影。3.2 命令行方式批量验证下载文件完整性如果下载的是分幅文件而不是合并结果拿到手后第一件事不是立刻处理而是先看范围、坐标系、无效值。用GDAL的命令行工具几个命令就能检查清楚# 查看栅格基本信息 gdalinfo N24E117.tif # 快速统计高程范围检测负值 gdalinfo -stats N24E117.tif # 输出为ASCII格式方便快速抽样 gdal_translate -of XYZ N24E117.tif sample.xyzgdalinfo输出的关键字段包括Size is 3601, 3601表示像素尺寸、Origin左上角坐标、Pixel Size网格间距ALOS是0.00027778度约30米12.5米版本则不同、TypeInt16高程值按整型存单位米。特别注意NoData Value字段——如果它是默认的-32768后面所有计算里必须显式指定它否则裁剪和叠加时会把海洋区域的高程拉进统计里。检查步骤里最容易忽略的一步是确认无效值的空间分布。海洋区域在ALOS原始数据里是有效的高程0米但OpenTopography输出的产品通常已经把海域设置为-32768或NaN因此厦门本岛周围的海面不会变成一马平川而是显示为无数据。这个特性在做淹没分析或海岸线提取时有用但在计算全区域高程均值时会造成严重偏移。3.3 用Python批量拼接厦门周边分幅瓦片如果不想依赖OpenTopography的在线拼接也需要在本地掌握批量拼接的脚本原因之一是万一最终下载结果只有单幅后续自行处理时用得到。厦门周边可能涉及N24E117、N24E118、N25E117、N25E118等分幅用rasterio.merge是最好的方式import rasterio from rasterio.merge import merge from rasterio.plot import show import glob # 列出所有ALOS分幅文件 tif_files glob.glob(/workspace/xiamen_dem/source/N*.tif) print(f共找到 {len(tif_files)} 个分幅文件) # 执行拼接注意输出路径和输出类型 mosaic, out_transform merge(tif_files, methodfirst) mosaic_meta rasterio.open(tif_files[0]).meta.copy() mosaic_meta.update({ driver: GTiff, height: mosaic.shape[1], width: mosaic.shape[2], transform: out_transform, dtype: int16 }) with rasterio.open(/workspace/xiamen_dem/source/xiamen_mosaic.tif, w, **mosaic_meta) as dest: dest.write(mosaic)这段代码的关键逻辑是merge默认把重叠区域的像素取第一个文件的数值methodfirst对于DEM来说重叠带两边是同一数据源的相邻瓦片取值差异通常在0.1米以内用first没问题。而mosaic_meta要从原始文件继承NoData和坐标信息不能直接新建一个空meta那样会丢坐标系。最后写入时dtype保持int16不要在拼接阶段就转Float32否则文件体积翻倍。4. 厦门DEM的裁剪、重投影与DSM到DEM的修正4.1 用行政边界裁剪DEM并处理边界锯齿拿到合并后的矩形DEM后要用厦门市行政边界把它裁剪出来。这里有个操作细节不能直接用gdalwarp -cutline默认参数因为默认情况下裁剪边界上的像素会做最邻近重采样导致边界线上出现「锯齿状」台阶。厦门海岸线和行政区界有大量曲线比如翔安海岸线需要指定-crop_to_cutline和合适的重采样算法gdalwarp -cutline /workspace/xiamen_dem/boundary/xiamen_boundary.shp \ -crop_to_cutline \ -tr 0.0002777777777777778 0.0002777777777777778 \ -r bilinear \ -of GTiff \ /workspace/xiamen_dem/source/xiamen_mosaic.tif \ /workspace/xiamen_dem/cut/xiamen_cut_30m.tif-tr的数值是目标分辨率的经纬度跨度0.00027777度约等于30米如果你用的是12.5米数据应该改成约0.00011574度1/86400度。用bilinear双线性插值配合crop_to_cutline边界上的高程过渡会比较平滑。gdalwarp后面可以加上--config GDALWARP_IGNORE_BAD_CUTLINE YES这样边界SHP即使有小瑕疵也不会中断整条命令但对厦门这种边界数据质量较高的场景不是必须的。裁剪完成后要验证两个东西一是裁剪后的像素行列数是否需要在后续处理中记录二是裁剪边缘的高程值有没有出现「一圈0米」的异常。厦门市界不是完全贴合海岸线边界外的海域应该保留NoData而不是0这在后续做表面积计算时才不会把海域算进去。4.2 重投影到厦门本地坐标系目标CRS怎么定DEM处理可以全程在WGS84经纬度里做但一旦涉及和已有CAD数据叠加比如热词里提到的「dem dxf 叠加 python」坐标统一就是绕不开的。厦门常见的本地坐标系有厦门94地方坐标系以及基于CGCS2000的3度带高斯投影中央经线118.5°E。在没有明确标注的情况下后台GIS数据通常采用CGCS2000 / 3-degree Gauss-Kruger zone 35或UTM zone 50N。用GDAL做重投影# 从WGS84重投影到CGCS2000 / 3-degree Gauss-Kruger CM 118.5E gdalwarp -t_srs EPSG:4549 \ -r cubic \ -ot Float32 \ /workspace/xiamen_dem/cut/xiamen_cut_30m.tif \ /workspace/xiamen_dem/final/xiamen_dem_cgcs2000.tif这里用EPSG:4549CGCS2000 / 3-degree Gauss-Kruger CM 118.5E是因为厦门正好落在118.5度中央经线上变形最小。重采样选cubic三次卷积因为经过投影转换后网格方向变了最邻近会产生轻微的马赛克感而双线性会对地形细节有一定柔化。如果最终产品只做高程统计不做视觉分析回到bilinear也可以速度更快。重投影后注意检查NoData值有没有被保留有些GDAL版本在投影变换后会把NoData设置为0需要重新gdal_translate -a_nodata -32768修正。4.3 DSM生成DEM的差值法厦门高密度城区的特殊处理厦门地形分析的另一个常见需求是「dsm生成dem」——把包含建筑和植被的表面模型DSM滤出真实地面DEM。ALOS AW3D30这类光学雷达数据本身已经是DSM能看出城市建筑的轮廓是要付出代价的在厦门岛内老城区和软件园三期这种高密度区域DSM的地表高程比真实地面高几米到十几米直接拿它做坡度分析或洪水淹没会得到错误结论。最常用的DSM到DEM思路是「最小值滤波 局部差值」import numpy as np from scipy.ndimage import minimum_filter import rasterio with rasterio.open(/workspace/xiamen_dem/final/xiamen_dem_cgcs2000.tif) as src: dem src.read(1).astype(float32) profile src.profile # 用半径约50米的窗口取最小值抑制建筑和树冠 radius_px max(1, int(50 / 30)) # 30米分辨率下约2像素 ground minimum_filter(dem, sizeradius_px * 2 1, modenearest) # 只在地表高出地面超过3米的像素上做替换 mask (dem - ground) 3 dem_corrected np.where(mask, ground, dem) # 写回文件 profile.update(dtypefloat32, nodata-32768.0) with rasterio.open(/workspace/xiamen_dem/final/xiamen_dem_ground.tif, w, **profile) as dst: dst.write(dem_corrected, 1)这个方法的物理逻辑是在一个窗口内建筑和树冠的顶部不是地形取最小值能找出裸地但厦门大量农田和裸地区域本身是平坦的直接替换会影响真实地形所以只有差值超过3米才修正。厦门岛内建筑普遍在3层以上3米阈值刚好把低矮平房和雨棚排除在外。这个步骤在30米分辨率上效果一般因为30米网格里建筑和道路混在一起最小值滤波容易把道路挖成沟而在12.5米数据上效果会好得多建议有条件的时候优先对12.5米做这个处理。5. 厦门地形分析落地坡度、坡向与等高线的参数设置5.1 在厦门这种起伏地形下算坡度的参数选择完成DEM数据准备后最常见的分析任务是坡度、坡向、等高线。在厦门这些分析结果直接用在环评、选址、水土保持上比如同安的茶山地块坡度大于25度通常会被划为限制建设区。用GDAL直接生成坡度的命令是通用的但参数设置关乎结果正确性# 生成坡度输出单位为度 gdaldem slope /workspace/xiamen_dem/final/xiamen_dem_cgcs2000.tif \ /workspace/xiamen_dem/final/xiamen_slope.tif \ -p \ -s 1.0 # 生成坡向北为0度顺时针 gdaldem aspect /workspace/xiamen_dem/final/xiamen_dem_cgcs2000.tif \ /workspace/xiamen_dem/final/xiamen_aspect.tif-p参数让坡度以度为单位输出如果不加则默认是百分数tan值乘以100在厦门这种有山有平地的场景下百分数不利于直接设定阈值。-s是垂直高程单位与水平单位的比值因为数据已经重投影到米制坐标水平单位是米垂直单位也是米所以设1.0。一个常见错误是拿WGS84经纬度的DEM直接算坡度且忘记加-s这时GDAL默认水平单位是度算出来的坡度完全失真厦门山体动辄出现80度以上的假值。注意如果你跳过重投影直接在WGS84坐标系下运行gdaldem必须设置-s 111320.01度对应的米数否则坡度结果没有参考价值。5.2 等高线生成的关键参数间距与平滑等高线在厦门常用于精细的场地设计常见做法是用gdal_contour从DEM里抽线gdal_contour -a ELEV \ -i 5.0 \ /workspace/xiamen_dem/final/xiamen_dem_cgcs2000.tif \ /workspace/xiamen_dem/final/xiamen_contour_5m.shp-a ELEV指定属性字段名后面接SHP时这个属性就代表每根等高线的高程值。-i 5.0是等高距厦门岛内地形平缓5米比较合适同安山区高差大如果做的是全区域检索图建议用10米否则SHP文件过大且线条密集难读。厦门填海区东渡、高崎、翔安南部的高程通常在2-5米之间5米等高距会在这类区域产生非常少的线但实际上这些区域几乎完全平坦线少是正常的。gdal_contour生成的线默认是折线不做平滑。如果交付给设计院通常还需要在QGIS里对线做一次平滑或者用v.generalizeGRASS做Simplification。但平滑操作只影响视觉效果不影响高程属性所以精度敏感的工程验收不建议平滑。5.3 常见误用用DSM代替DEM做坡度分析的后果在厦门这类高密度城市直接用原始DSM做坡度分析会得到一个「非常难看」的结果岛内建筑密集区坡度普遍在10度以上看起来整个厦门老城区都建在山坡上。这显然是错的。用第4.3节修正后的地面DEM重新算一遍坡度同样的区域坡度会降到3度以下才符合真实地形感受。这个问题的根源在于DSM本身是「地表覆盖物顶部」的模型建筑屋顶、树木冠层都是它的组成部分。厦门绿化覆盖率高仙岳路沿线、狐尾山周边的树冠在高程模型上是连续起伏的「假山丘」直接分析会把真正的地形沟谷和排水路径全部掩盖。因此任何以「工程建设」为目的的分析都必须先做DSM到DEM的差值修正再出坡度和等高线。可视化或飞行漫游则不必做这个修正保留建筑轮廓反而信息量更大。6. 最后收口高程基准校验与DXF叠加的3个实用技巧6.1 用一个值得信赖的反例校验厦门附近GNSS水准点核对拿到处理完的厦门DEM我一般先做一次「单点校验」。厦门的实测高程控制点分布在岛内外不同区域从这些点里挑三到四个写入一个CSV包含经纬度和实测海拔然后用gdal_query.py或简单的rasterio采样脚本把DEM对应位置的高程取出来对比import csv import rasterio with rasterio.open(/workspace/xiamen_dem/final/xiamen_dem_ground.tif) as src: for row in csv.DictReader(open(control_points.csv)): lon, lat float(row[lon]), float(row[lat]) # 反算像素坐标 py, px src.index(lon, lat) dem_val src.read(1)[py, px] print(f{row[name]}: 实测 {row[h]}mDEM {dem_val}m差 {float(row[h])-dem_val:.2f}m)如果差值平均在±2米以内说明数据的垂直质量合格如果系统性偏向一侧比如所有点都比实测高2米基本可以判断是基准面差异而不是数据处理错误。厦门同安老城区的控制点数据常有基准不统一的情况碰到系统性误差时不要急着改数据先确认对照点使用的高程基准是否和DEM的EGM96一致。6.2 将DEM与DXF叠加的坐标系匹配技巧「dem dxf 叠加 python」是很多市政项目的刚需把CAD设计图通常带本地坐标叠加到DEM渲染的底图上做方案比对。这一步最常见的问题是两者相差几米到几十米不等对齐不上。此时用Python比在QGIS里手动拖拽来得可靠import geopandas as gpd import rasterio from rasterio.plot import show import matplotlib.pyplot as plt # 读取DXF指定EPSG:4549厦门常用CRS gdf gpd.read_file(/workspace/xiamen_dem/dwg/design.dxf) gdf gdf.to_crs(EPSG:4549) # 读取已重投影的DEM with rasterio.open(/workspace/xiamen_dem/final/xiamen_dem_cgcs2000.tif) as src: fig, ax plt.subplots(figsize(10, 10)) show(src, axax, cmapterrain) gdf.plot(axax, colorred, linewidth0.8) plt.savefig(/workspace/xiamen_dem/final/overlay_check.png, dpi150)如果dxf里的坐标不是米制而是毫米或者英尺转换时要在to_crs前用gdf.scale(xfact0.001, yfact0.001, origin(0,0))缩放。厦门本地院出图习惯不同这步检查必不可少。叠加后除了看套合度还要看DEM渲染的透明度和色带是否符合判别需求这些由参数补足。6.3 批量导出特定区域小范围地形切片最后送一个常用技巧从整市DEM里切出一小块矩形区域转成DXF或XYZ格式给无人机或其他专业软件用。gdal_translate配合-projwin参数可以做按坐标范围切片-of XYZ输出高度精简的三列文本同时这也是最容易读出来的格式gdal_translate -projwin 13425070 2751350 13428000 2748000 \ -of XYZ /workspace/xiamen_dem/final/xiamen_dem_ground.tif \ /workspace/xiamen_dem/final/slice_xiamen_island.xyz-projwin后跟的是左上角x、左上角y、右下角x、右下角y注意是「上」在前「下」在后很多人写成左右顺序结果切出反区域。切片完成后用head或wc -l快速检查点数确保没有把海洋区域的大面积NoData带进来。这样一整套流程从下载、预处理、修正到分析、校验就闭环了厦门的DEM数据在任何需要高程信息的GIS工程里都能直接成为可信的底图而不是一堆难以收拾的原始tif。本文还有配套的精品资源点击获取