简介这份武汉市30米分辨率DEM数字高程数据包以栅格形式刻画了武汉市及周边区域的地表起伏同时附带完整的市界矢量文件便于在GIS中直接进行地形分析、规划选址与环境评估面向地理信息从业者、城市规划师及高校相关专业师生。压缩包共12个文件整体大小25.27MB除核心TIFF高程数据外还配齐了Shapefile边界、投影定义、属性表、空间索引与元数据文件栅格与矢量两类数据相互配套可满足快速裁剪、叠加分析及专题制图等常见需求。目前已有2069人学习下载。借助ArcGIS、QGIS等软件即可加载使用文件结构完整、坐标系信息明确省去自行配准和整理数据的步骤可直接用于坡度坡向提取、等高线生成、汇水区分析、视域分析等多种地形建模场景实用性强。1. 武汉30m DEM数据包从解压到用得上的完整路径一份「湖北省武汉市DEM数字高程数据30m含本市级范围shp文件.zip」能做的事比大多数人想象的多用它做洪水淹没分析、日照分析、建筑选址坡度评估、通信基站视距仿真甚至直接给三维城市模型当基底。zip里装的是一张武汉市范围的30米分辨率数字高程模型外加一个市级行政区划边界shp。30m这个数字意味着每个栅格像元对应地面30米乘30米武汉全市域大约能铺出千万级像元文件不大但精度足够做区域级分析。适合谁用GIS工程师、城市规划相关方向的研究生、做城市级数据产品的人。拿到手第一件事不该是双击打开而是先想清楚坐标系、裁剪边界和有效值范围这组问题。2. DEM与30m分辨率的实际含义先看懂再动手2.1 DEM的三种常见形态和30m数据的定位DEMDigital Elevation Model在学术和工程语境里有三种常被混用的形态DSMDigital Surface Model记录地表以上所有物体的顶面包括树冠、屋顶DTMDigital Terrain Model剥离地物只留裸地面DEM一词在实际数据包里通常指代DTM但严格说它只是「高程栅格」的统称。这份标题写明「数字高程数据30m」没有写DSM所以按裸地面高程理解是对的。分辨率30m是区域分析的黄金档位。ALOS 12.5m数据和NASA的SRTM 30m是开源高程数据里最常用的两组后者的全球覆盖度和精度校验都更成熟国内很多城市级项目直接把30m当作默认底图。12.5m在城区能看到更细的街道切割痕迹但数据量翻倍处理速度下降对山坡地这种高程变化剧烈的地方12.5m与30m的坡度计算结果差距会达到5度以上。如果只做市级宏观分析30m完全够用没必要追高分辨率。反直觉的一点是30m DEM做坡度分析时大于30度的坡面会被系统性低估因为栅格化天然抹平了细碎地形。真要做坡耕地识别或地质灾害危险性评价至少要用12.5m甚至1m的数据。2.2 本市级范围shp文件在数据包里的角色标题括号里的「含本市级范围shp文件」是这份数据最值钱的部分。DEM栅格是矩形的必然包含武汉市行政边界以外的区域分析时如果直接拿整个矩形范围算面积、算体积结果一定偏大。shp文件就是用来做掩膜裁剪的裁完之后的统计结果才和武汉市的实际辖区对得上。shp不是一个文件而是一组文件的集合最少包含四个扩展名作用缺失后果.shp几何信息点线面的坐标无法显示图形.shx几何索引部分软件打不开.dbf属性表行政区名称等字段属性读取失败.prj坐标系定义WKT格式坐标系错乱或无法叠加拷贝或传输时把这四个文件漏掉任何一个都会导致打不开或坐标漂移。zip打包正好规避了这个问题这也是这类数据多用zip分发的原因。解压后如果发现 .prj 缺失最常见的后果是QGIS和ArcGIS把它当未知坐标系叠加DEM时几何位置对不上。2.2.1 判断shp坐标系是否与DEM一致在QGIS里加载shp后右键图层打开属性在「信息」标签页里看「CRS」一栏。国内市级数据大概率是CGCS2000 / 3-degree Gauss-Kruger zone 38EPSG:4547或4548取决于经度带也可能是WGS84经纬度EPSG:4326。DEM的坐标系同样在图层属性里查看。两者不一致时先做重投影再裁剪顺序不能反先统一坐标系后执行裁剪。2.3 坐标系选错是90%的叠加错位根源30m DEM如果坐标系是经纬度一个像元的跨度就是0.00027778度左右如果投影坐标系是CGCS2000高斯克吕格像元对应的就是地板平整的30米。做面积、坡度、体积计算必须用投影坐标系因为经纬度不是等距的——在武汉这个纬度上1度经度的实地距离约为100公里1度纬度约111公里直接用经纬度算距离会偏大约10%。用投影坐标系的逻辑是把地球曲面展平成平面让每个像元对应真实的30米乘以30米。处理流程里我一般先跑一行GDAL命令确认两件事一个看坐标系一个看像元大小gdalinfo DEM_merge.tif | grep -E Pixel Size|Coordinate System|EPSG|Size is输出里Pixel Size (30.000000, -30.000000)表示正南北的30m像元负号是正常的代表行号向下时纬度递减。Coordinate System is: PROJCRS[CGCS2000 / 3-degree Gauss-Kruger zone 38...]这类描述说明栅格已经带投影。如果看到GEOGCRS[WGS 84]且 Pixel Size 是0.000几的小数说明还是经纬度数据后续算面积前必须先投影。3. 解压之后先做数据体检完整性、可视化与无效值3.1 用命令行检查zip完整性和内部清单拿到zip先不急着解压先验证压缩包的完整性。zip本身有CRC32校验解压工具会自动校验但命令行能看得更清楚# 列出压缩包内容不实际解压 unzip -l 武汉DEM30m.zip # 测试压缩包完整性 unzip -t 武汉DEM30m.zip-l参数会列出所有文件、原始大小和压缩后大小。重点看是否有.tfw世界文件或.aux.xml辅助文件——有.tfw说明栅格自带地理配准信息.aux.xml里存了统计信息和无效值定义对后续处理很有帮助。-t参数执行CRC校验看到No errors detected再往下走。这一步不是洁癖而是zip在传输过程中丢数据的概率远高于本地复制栅格数据丢几个字节在视觉上可能看不出来但做统计分析时会出现异常值。解压后的一等公民操作是统计有效值范围gdalinfo -stats 武汉DEM30m.tif | grep -E Minimum|Maximum|NoData-stats会让GDAL先扫描一遍全图把最小值、最大值、均值和标准差写进元数据。如果最小值是 -9999 或 0 附近的大片区域说明有NoData值混进来了后面做坡度分析前必须屏蔽。这一步和shp转txt之类的操作无关纯粹是数据质量的底线。3.2 在QGIS中把DEM和shp叠在一起看不需要写任何代码打开QGIS拖拽DEM栅格文件进入图层窗口再拖入shp边界文件。栅格默认用灰度渲染一片黑需要调整双击图层打开「符号体系」渲染类型选「单波段假彩色」色带选「Spectral」或「Turbo」色带方向按「值递减」反转让低海拔显示深色、高海拔显示亮色。这样武汉的南高北低地势结构立即显形——江夏、黄陂北部的低山丘陵区会明显凸起长江汉江冲积平原一片平坦。shp边界叠加后如果边界线正好贴合DEM里的地形突变处说明坐标匹配良好如果边界跑到影像外面或者歪斜优先怀疑坐标系不一致。QGIS里可以随时查看每个图层的CRS不需要额外插件。小技巧把shp边界设成半透明红色在DEM渲染之上看得最清楚。3.3 无效值和异常高程的处理顺序DEM数据里常见的异常点表现为「盐粒」——某个像元高程比周围高出几百米或低出几百米视觉上像一个亮斑。原因是原始数据生产中的回波异常或插值失败。处理顺序建议是先用gdalinfo确认NoData值再做一个简单的邻域滤波gdal_fillnodata.py 武汉DEM30m.tif 武汉DEM30m_fill.tif -md 5-md参数设定最大搜索距离单位是像元5表示以异常点为中心向外搜索5个像元半径之内的有效值来插值。这个工具只处理NoData区域不管「非NoData但明显异常」的盐粒。粗查盐粒的方法是加载到QGIS里拉大对比度或者用Python读数组后看分位数这一步在第4章代码里展开。4. 用Python实现武汉DEM的裁剪、重投影与地形分析4.1 用rasterio和geopandas按市级边界批量裁剪Python处理DEM的首选组合是rasterio负责栅格、geopandas负责矢量。完整流程是读shp边界、统一坐标系、按边界裁剪、输出结果。示例代码import rasterio from rasterio.mask import mask import geopandas as gpd # 读取shp边界武汉市级范围 boundary gpd.read_file(武汉市级边界.shp) # 读取DEM栅格 with rasterio.open(武汉DEM30m.tif) as src: # 边界必须和栅格坐标系一致不一致时先重投影 if boundary.crs ! src.crs: boundary boundary.to_crs(src.crs) # 用边界几何对象裁剪栅格 out_image, out_transform mask( src, boundary.geometry, cropTrue, # 裁剪到边界严格外接矩形 nodatasrc.nodata # 保留原来的NoData值 ) # 更新元数据并写新文件 out_meta src.meta.copy() out_meta.update({ driver: GTiff, height: out_image.shape[1], width: out_image.shape[2], transform: out_transform }) with rasterio.open(武汉DEM_裁剪.tif, w, **out_meta) as dst: dst.write(out_image)to_crs只做矢量重投影比整幅栅格重投影快得多所以先统一边界坐标再让rasterio在裁剪时自动处理像元对齐。mask函数的cropTrue是关键参数它会自动把输出范围收缩到边界的矩形外框上不设这个参数的话输出范围和原栅格一样大只是边界外变NoData。裁剪结果的实际几何形状仍然是不规则的因为边界多边形是任意形状矩形外框之外的区域会被NoData填充。4.2 坡度、坡向和山体阴影的参数选择裁完之后的DEM可以直接做地形因子分析。GDAL自带gdaldem工具不需要装额外库# 坡度单位是度 gdaldem slope 武汉DEM_裁剪.tif 武汉DEM_坡度.tif -p -s 1.0 # 坡向0-360度 gdaldem aspect 武汉DEM_裁剪.tif 武汉DEM_坡向.tif # 山体阴影用于可视化底图 gdaldem hillshade 武汉DEM_裁剪.tif 武汉DEM_阴影.tif -z 2.0 -az 315 -alt 45-p让坡度输出单位为度而不是百分数更直观-s是垂直比例因子平地数据设为1.0即可但如果数据是经纬度坐标系这里必须改成实际的尺度因子约为 111320 / 90000否则坡度会失真。gdaldem的算法用的是3x3窗口的中心差分法边缘像元会自动扩展不需要预填充。hillshade里-az是太阳方位角315度是西北方向-alt是太阳高度角45度模拟春秋季日照做展示时这两个参数调出来的阴影效果比默认值立体得多。4.3 用numpy做体积和面积统计裁完后用numpy直接算填挖方量是DEM最常用到的需求之一比如估算某个区域平整场地的土方量import numpy as np import rasterio with rasterio.open(武汉DEM_裁剪.tif) as src: dem src.read(1).astype(np.float32) nodata src.nodata pixel_width src.transform[0] # 像元宽度 pixel_height -src.transform[4] # 像元高度负号是因为Y轴向下 # 屏蔽NoData valid_mask dem ! nodata dem np.where(valid_mask, dem, np.nan) pixel_area pixel_width * pixel_height # 假设设计高程为25米计算填挖方 design_level 25.0 diff dem - design_level fill_volume np.nansum(np.where(diff 0, -diff, 0)) * pixel_area cut_volume np.nansum(np.where(diff 0, diff, 0)) * pixel_area print(f填方量{fill_volume / 1e6:.2f} 百万立方米) print(f挖方量{cut_volume / 1e6:.2f} 百万立方米)src.transform[0]是X方向像元分辨率transform[4]是Y方向分辨率因为GDAL的transform中第5个元素索引4是负的表示向下递减取绝对值得到正的高度值。整个逻辑是先计算每个像元与设计高程的差值正值是挖负值是填再乘以每个像元的实地面积后累加。这套代码可以直接替换成任意设计高程例如做水库淹没分析时把design_level换成水位高程输出结果就是淹没区域的库容量。5. 三个高频踩坑点与数据自检技巧5.1 shp属性表乱码问题国内shp的dbf属性表用GBK编码存储QGIS和ArcGIS默认按UTF-8读取时会出现「武汉市」显示为乱码的情况。解决方式不是改数据而是让软件正确识别编码。QGIS中双击图层在「数据源」选项卡里找到「数据源编码」手动改成GBK或CP936属性表立即正常。命令行下用ogrinfo -so -al 武汉市级边界.shp也可以查看到属性字段内容如果输出乱码同样加一句-oo ENCODINGGBK。这个坑高发在从国内数据商或政府平台下载的数据包里。5.2 裁剪后边缘出现一圈黑色NoData用QGIS的「裁剪栅格」工具或第4章的Python代码裁剪后输出图像四周可能出现一圈黑色边界。黑色区域是NoData栅格裁剪的本质是用矩形外框承装不规则几何边界。解决思路有两个层面。视觉层面在QGIS渲染设置里把NoData颜色设为透明右键图层属性「透明度」选项卡里勾选「NoData值」填入栅格的NoData数值。分析层面做统计计算时用np.nansum或 GDAL 的-srcnodata参数显式屏蔽NoData避免把NoData当成0值计入。0值在高程数据里可能代表海平面或真实平地混在一起会污染统计结果。5.3 像元对齐与重采样方法选择当有多幅DEM需要镶嵌比如武昌、汉阳、汉口三个子区域各自的DEM时两幅图之间的重叠区可能出现像元错位。原因是不同子图可能来自不同轨道或生产批次像元的网格原点有微小偏移。最稳的解法是设定统一的目标网格再做重采样gdalwarp -t_srs EPSG:4547 -tr 30 30 -r cubic -dstnodata -9999 \ 武昌.tif 汉阳.tif 汉口.tif 武汉镶嵌.tif-tr 30 30强制输出像元为30米网格-r cubic选择三次卷积重采样比默认的near算法在连续地形表面上少很多锯齿但计算量大约多三倍。-dstnodata把输出NoData统一为-9999避免三幅图各带各的NoData值。这里有一个重要取舍如果后续要做坡度分析建议用-r bilinear或-r cubic因为near会在斜坡上产生阶梯状伪影如果只是做分类比如按高程分带near反而更安全不会生成原数据里不存在的高程值。5.4 数据自检的半小时清单最终交付前用一组固定流程做数据质量自检。第一查直方图在QGIS里打开栅格直方图如果分布是平滑的单峰或双峰基本正常如果出现「远离主峰的孤立细柱」说明有坏像元。第二查边界对齐把裁好的DEM和原始shp边界叠加放大到1:5000检查边界线两侧是否有带状NoData正常裁剪边界贴合度应该在一个像元以内。第三查坐标系用gdalinfo查看输出文件的EPSG代码确认没有在中间步骤丢失.prj信息。第四查文件大小30m分辨率的武汉市域未压缩GeoTIFF大约在50-150MB之间取决于是否带LZW压缩如果只有几MB大概率是裁剪过头或压缩过度。这套流程跑完数据就可以放心进入后续分析管线了。本文还有配套的精品资源点击获取 武汉30m DEM数据处理全流程：从解压到地形分析
简介这份武汉市30米分辨率DEM数字高程数据包以栅格形式刻画了武汉市及周边区域的地表起伏同时附带完整的市界矢量文件便于在GIS中直接进行地形分析、规划选址与环境评估面向地理信息从业者、城市规划师及高校相关专业师生。压缩包共12个文件整体大小25.27MB除核心TIFF高程数据外还配齐了Shapefile边界、投影定义、属性表、空间索引与元数据文件栅格与矢量两类数据相互配套可满足快速裁剪、叠加分析及专题制图等常见需求。目前已有2069人学习下载。借助ArcGIS、QGIS等软件即可加载使用文件结构完整、坐标系信息明确省去自行配准和整理数据的步骤可直接用于坡度坡向提取、等高线生成、汇水区分析、视域分析等多种地形建模场景实用性强。1. 武汉30m DEM数据包从解压到用得上的完整路径一份「湖北省武汉市DEM数字高程数据30m含本市级范围shp文件.zip」能做的事比大多数人想象的多用它做洪水淹没分析、日照分析、建筑选址坡度评估、通信基站视距仿真甚至直接给三维城市模型当基底。zip里装的是一张武汉市范围的30米分辨率数字高程模型外加一个市级行政区划边界shp。30m这个数字意味着每个栅格像元对应地面30米乘30米武汉全市域大约能铺出千万级像元文件不大但精度足够做区域级分析。适合谁用GIS工程师、城市规划相关方向的研究生、做城市级数据产品的人。拿到手第一件事不该是双击打开而是先想清楚坐标系、裁剪边界和有效值范围这组问题。2. DEM与30m分辨率的实际含义先看懂再动手2.1 DEM的三种常见形态和30m数据的定位DEMDigital Elevation Model在学术和工程语境里有三种常被混用的形态DSMDigital Surface Model记录地表以上所有物体的顶面包括树冠、屋顶DTMDigital Terrain Model剥离地物只留裸地面DEM一词在实际数据包里通常指代DTM但严格说它只是「高程栅格」的统称。这份标题写明「数字高程数据30m」没有写DSM所以按裸地面高程理解是对的。分辨率30m是区域分析的黄金档位。ALOS 12.5m数据和NASA的SRTM 30m是开源高程数据里最常用的两组后者的全球覆盖度和精度校验都更成熟国内很多城市级项目直接把30m当作默认底图。12.5m在城区能看到更细的街道切割痕迹但数据量翻倍处理速度下降对山坡地这种高程变化剧烈的地方12.5m与30m的坡度计算结果差距会达到5度以上。如果只做市级宏观分析30m完全够用没必要追高分辨率。反直觉的一点是30m DEM做坡度分析时大于30度的坡面会被系统性低估因为栅格化天然抹平了细碎地形。真要做坡耕地识别或地质灾害危险性评价至少要用12.5m甚至1m的数据。2.2 本市级范围shp文件在数据包里的角色标题括号里的「含本市级范围shp文件」是这份数据最值钱的部分。DEM栅格是矩形的必然包含武汉市行政边界以外的区域分析时如果直接拿整个矩形范围算面积、算体积结果一定偏大。shp文件就是用来做掩膜裁剪的裁完之后的统计结果才和武汉市的实际辖区对得上。shp不是一个文件而是一组文件的集合最少包含四个扩展名作用缺失后果.shp几何信息点线面的坐标无法显示图形.shx几何索引部分软件打不开.dbf属性表行政区名称等字段属性读取失败.prj坐标系定义WKT格式坐标系错乱或无法叠加拷贝或传输时把这四个文件漏掉任何一个都会导致打不开或坐标漂移。zip打包正好规避了这个问题这也是这类数据多用zip分发的原因。解压后如果发现 .prj 缺失最常见的后果是QGIS和ArcGIS把它当未知坐标系叠加DEM时几何位置对不上。2.2.1 判断shp坐标系是否与DEM一致在QGIS里加载shp后右键图层打开属性在「信息」标签页里看「CRS」一栏。国内市级数据大概率是CGCS2000 / 3-degree Gauss-Kruger zone 38EPSG:4547或4548取决于经度带也可能是WGS84经纬度EPSG:4326。DEM的坐标系同样在图层属性里查看。两者不一致时先做重投影再裁剪顺序不能反先统一坐标系后执行裁剪。2.3 坐标系选错是90%的叠加错位根源30m DEM如果坐标系是经纬度一个像元的跨度就是0.00027778度左右如果投影坐标系是CGCS2000高斯克吕格像元对应的就是地板平整的30米。做面积、坡度、体积计算必须用投影坐标系因为经纬度不是等距的——在武汉这个纬度上1度经度的实地距离约为100公里1度纬度约111公里直接用经纬度算距离会偏大约10%。用投影坐标系的逻辑是把地球曲面展平成平面让每个像元对应真实的30米乘以30米。处理流程里我一般先跑一行GDAL命令确认两件事一个看坐标系一个看像元大小gdalinfo DEM_merge.tif | grep -E Pixel Size|Coordinate System|EPSG|Size is输出里Pixel Size (30.000000, -30.000000)表示正南北的30m像元负号是正常的代表行号向下时纬度递减。Coordinate System is: PROJCRS[CGCS2000 / 3-degree Gauss-Kruger zone 38...]这类描述说明栅格已经带投影。如果看到GEOGCRS[WGS 84]且 Pixel Size 是0.000几的小数说明还是经纬度数据后续算面积前必须先投影。3. 解压之后先做数据体检完整性、可视化与无效值3.1 用命令行检查zip完整性和内部清单拿到zip先不急着解压先验证压缩包的完整性。zip本身有CRC32校验解压工具会自动校验但命令行能看得更清楚# 列出压缩包内容不实际解压 unzip -l 武汉DEM30m.zip # 测试压缩包完整性 unzip -t 武汉DEM30m.zip-l参数会列出所有文件、原始大小和压缩后大小。重点看是否有.tfw世界文件或.aux.xml辅助文件——有.tfw说明栅格自带地理配准信息.aux.xml里存了统计信息和无效值定义对后续处理很有帮助。-t参数执行CRC校验看到No errors detected再往下走。这一步不是洁癖而是zip在传输过程中丢数据的概率远高于本地复制栅格数据丢几个字节在视觉上可能看不出来但做统计分析时会出现异常值。解压后的一等公民操作是统计有效值范围gdalinfo -stats 武汉DEM30m.tif | grep -E Minimum|Maximum|NoData-stats会让GDAL先扫描一遍全图把最小值、最大值、均值和标准差写进元数据。如果最小值是 -9999 或 0 附近的大片区域说明有NoData值混进来了后面做坡度分析前必须屏蔽。这一步和shp转txt之类的操作无关纯粹是数据质量的底线。3.2 在QGIS中把DEM和shp叠在一起看不需要写任何代码打开QGIS拖拽DEM栅格文件进入图层窗口再拖入shp边界文件。栅格默认用灰度渲染一片黑需要调整双击图层打开「符号体系」渲染类型选「单波段假彩色」色带选「Spectral」或「Turbo」色带方向按「值递减」反转让低海拔显示深色、高海拔显示亮色。这样武汉的南高北低地势结构立即显形——江夏、黄陂北部的低山丘陵区会明显凸起长江汉江冲积平原一片平坦。shp边界叠加后如果边界线正好贴合DEM里的地形突变处说明坐标匹配良好如果边界跑到影像外面或者歪斜优先怀疑坐标系不一致。QGIS里可以随时查看每个图层的CRS不需要额外插件。小技巧把shp边界设成半透明红色在DEM渲染之上看得最清楚。3.3 无效值和异常高程的处理顺序DEM数据里常见的异常点表现为「盐粒」——某个像元高程比周围高出几百米或低出几百米视觉上像一个亮斑。原因是原始数据生产中的回波异常或插值失败。处理顺序建议是先用gdalinfo确认NoData值再做一个简单的邻域滤波gdal_fillnodata.py 武汉DEM30m.tif 武汉DEM30m_fill.tif -md 5-md参数设定最大搜索距离单位是像元5表示以异常点为中心向外搜索5个像元半径之内的有效值来插值。这个工具只处理NoData区域不管「非NoData但明显异常」的盐粒。粗查盐粒的方法是加载到QGIS里拉大对比度或者用Python读数组后看分位数这一步在第4章代码里展开。4. 用Python实现武汉DEM的裁剪、重投影与地形分析4.1 用rasterio和geopandas按市级边界批量裁剪Python处理DEM的首选组合是rasterio负责栅格、geopandas负责矢量。完整流程是读shp边界、统一坐标系、按边界裁剪、输出结果。示例代码import rasterio from rasterio.mask import mask import geopandas as gpd # 读取shp边界武汉市级范围 boundary gpd.read_file(武汉市级边界.shp) # 读取DEM栅格 with rasterio.open(武汉DEM30m.tif) as src: # 边界必须和栅格坐标系一致不一致时先重投影 if boundary.crs ! src.crs: boundary boundary.to_crs(src.crs) # 用边界几何对象裁剪栅格 out_image, out_transform mask( src, boundary.geometry, cropTrue, # 裁剪到边界严格外接矩形 nodatasrc.nodata # 保留原来的NoData值 ) # 更新元数据并写新文件 out_meta src.meta.copy() out_meta.update({ driver: GTiff, height: out_image.shape[1], width: out_image.shape[2], transform: out_transform }) with rasterio.open(武汉DEM_裁剪.tif, w, **out_meta) as dst: dst.write(out_image)to_crs只做矢量重投影比整幅栅格重投影快得多所以先统一边界坐标再让rasterio在裁剪时自动处理像元对齐。mask函数的cropTrue是关键参数它会自动把输出范围收缩到边界的矩形外框上不设这个参数的话输出范围和原栅格一样大只是边界外变NoData。裁剪结果的实际几何形状仍然是不规则的因为边界多边形是任意形状矩形外框之外的区域会被NoData填充。4.2 坡度、坡向和山体阴影的参数选择裁完之后的DEM可以直接做地形因子分析。GDAL自带gdaldem工具不需要装额外库# 坡度单位是度 gdaldem slope 武汉DEM_裁剪.tif 武汉DEM_坡度.tif -p -s 1.0 # 坡向0-360度 gdaldem aspect 武汉DEM_裁剪.tif 武汉DEM_坡向.tif # 山体阴影用于可视化底图 gdaldem hillshade 武汉DEM_裁剪.tif 武汉DEM_阴影.tif -z 2.0 -az 315 -alt 45-p让坡度输出单位为度而不是百分数更直观-s是垂直比例因子平地数据设为1.0即可但如果数据是经纬度坐标系这里必须改成实际的尺度因子约为 111320 / 90000否则坡度会失真。gdaldem的算法用的是3x3窗口的中心差分法边缘像元会自动扩展不需要预填充。hillshade里-az是太阳方位角315度是西北方向-alt是太阳高度角45度模拟春秋季日照做展示时这两个参数调出来的阴影效果比默认值立体得多。4.3 用numpy做体积和面积统计裁完后用numpy直接算填挖方量是DEM最常用到的需求之一比如估算某个区域平整场地的土方量import numpy as np import rasterio with rasterio.open(武汉DEM_裁剪.tif) as src: dem src.read(1).astype(np.float32) nodata src.nodata pixel_width src.transform[0] # 像元宽度 pixel_height -src.transform[4] # 像元高度负号是因为Y轴向下 # 屏蔽NoData valid_mask dem ! nodata dem np.where(valid_mask, dem, np.nan) pixel_area pixel_width * pixel_height # 假设设计高程为25米计算填挖方 design_level 25.0 diff dem - design_level fill_volume np.nansum(np.where(diff 0, -diff, 0)) * pixel_area cut_volume np.nansum(np.where(diff 0, diff, 0)) * pixel_area print(f填方量{fill_volume / 1e6:.2f} 百万立方米) print(f挖方量{cut_volume / 1e6:.2f} 百万立方米)src.transform[0]是X方向像元分辨率transform[4]是Y方向分辨率因为GDAL的transform中第5个元素索引4是负的表示向下递减取绝对值得到正的高度值。整个逻辑是先计算每个像元与设计高程的差值正值是挖负值是填再乘以每个像元的实地面积后累加。这套代码可以直接替换成任意设计高程例如做水库淹没分析时把design_level换成水位高程输出结果就是淹没区域的库容量。5. 三个高频踩坑点与数据自检技巧5.1 shp属性表乱码问题国内shp的dbf属性表用GBK编码存储QGIS和ArcGIS默认按UTF-8读取时会出现「武汉市」显示为乱码的情况。解决方式不是改数据而是让软件正确识别编码。QGIS中双击图层在「数据源」选项卡里找到「数据源编码」手动改成GBK或CP936属性表立即正常。命令行下用ogrinfo -so -al 武汉市级边界.shp也可以查看到属性字段内容如果输出乱码同样加一句-oo ENCODINGGBK。这个坑高发在从国内数据商或政府平台下载的数据包里。5.2 裁剪后边缘出现一圈黑色NoData用QGIS的「裁剪栅格」工具或第4章的Python代码裁剪后输出图像四周可能出现一圈黑色边界。黑色区域是NoData栅格裁剪的本质是用矩形外框承装不规则几何边界。解决思路有两个层面。视觉层面在QGIS渲染设置里把NoData颜色设为透明右键图层属性「透明度」选项卡里勾选「NoData值」填入栅格的NoData数值。分析层面做统计计算时用np.nansum或 GDAL 的-srcnodata参数显式屏蔽NoData避免把NoData当成0值计入。0值在高程数据里可能代表海平面或真实平地混在一起会污染统计结果。5.3 像元对齐与重采样方法选择当有多幅DEM需要镶嵌比如武昌、汉阳、汉口三个子区域各自的DEM时两幅图之间的重叠区可能出现像元错位。原因是不同子图可能来自不同轨道或生产批次像元的网格原点有微小偏移。最稳的解法是设定统一的目标网格再做重采样gdalwarp -t_srs EPSG:4547 -tr 30 30 -r cubic -dstnodata -9999 \ 武昌.tif 汉阳.tif 汉口.tif 武汉镶嵌.tif-tr 30 30强制输出像元为30米网格-r cubic选择三次卷积重采样比默认的near算法在连续地形表面上少很多锯齿但计算量大约多三倍。-dstnodata把输出NoData统一为-9999避免三幅图各带各的NoData值。这里有一个重要取舍如果后续要做坡度分析建议用-r bilinear或-r cubic因为near会在斜坡上产生阶梯状伪影如果只是做分类比如按高程分带near反而更安全不会生成原数据里不存在的高程值。5.4 数据自检的半小时清单最终交付前用一组固定流程做数据质量自检。第一查直方图在QGIS里打开栅格直方图如果分布是平滑的单峰或双峰基本正常如果出现「远离主峰的孤立细柱」说明有坏像元。第二查边界对齐把裁好的DEM和原始shp边界叠加放大到1:5000检查边界线两侧是否有带状NoData正常裁剪边界贴合度应该在一个像元以内。第三查坐标系用gdalinfo查看输出文件的EPSG代码确认没有在中间步骤丢失.prj信息。第四查文件大小30m分辨率的武汉市域未压缩GeoTIFF大约在50-150MB之间取决于是否带LZW压缩如果只有几MB大概率是裁剪过头或压缩过度。这套流程跑完数据就可以放心进入后续分析管线了。本文还有配套的精品资源点击获取