ARTICLE DETAIL

资讯详情

深耕郑州网站建设与运营推广的一线实战洞察。

邯郸市30米DEM数据处理全流程:从解压到地形分析

邯郸市30米DEM数据处理全流程:从解压到地形分析 简介河北省邯郸市DEM数字高程数据由三十米分辨率的栅格高程模型和行政边界矢量文件构成覆盖邯郸市及周边区域适合GIS从业者、城乡规划人员、测绘人员及高校相关专业学生使用每个像元代表三十米乘三十米区域的平均海拔高度可支撑城市规划、环境研究、灾害预防、农业区划等地理分析任务。压缩包内共含十二个文件核心为邯郸市高程栅格文件与邯郸市范围矢量边界文件附带投影定义、坐标参考、金字塔、元数据、属性表、空间索引等配套内容整体约二十九点三兆结构清晰便于直接导入主流GIS平台。使用者可完成等高线生成、坡度坡向分析、剖面量测、汇水区提取、淹没模拟、选址评价等操作边界矢量文件用于裁剪、叠加统计和专题制图配套元数据与投影文件则有助于正确配准和追溯数据来源。既可快速开展地形可视化教学也能支撑流域分析、工程选址等进阶建模已有八百九十一人学习或下载适合项目申报、课题研究和制图实验等场景。1. 邯郸市 DEM 数据包里到底装了什么拿到「河北省邯郸市DEM数字高程数据含区域范围shp文件.zip」后如果你只把里面的「邯郸市DEM.tif」拖进 QGIS 就开工等于把整套数据的价值用掉一半。这个压缩包里真正值钱的除了主栅格还有配套的 tfw 世界文件、ovr 金字塔、aux.xml 元数据以及一组描述邯郸市行政边界的 shp 文件组合。30 米分辨率的意思是每个像元代表地面 30 米乘 30 米区域的平均海拔这个精度适合做邯郸市全境的坡度分析、洪水淹没范围估算、选址评估和三维地形可视化也能和城市规划、测绘、地质调查里的常见矢量数据直接叠加。对 GIS 开发、测绘工程师和做空间分析的从业者来说这是一份可以直接复用的基础地形数据而不是只能看一眼的静态图。2. 拆解 ziptif、tfw、ovr 和 shp 侧车文件的作用2.1 栅格数据家族的“标配文件”各有分工解压 zip 后你会看到「邯郸市DEM.tif」「邯郸市DEM.tfw」「邯郸市DEM.tif.ovr」「邯郸市DEM.tif.aux.xml」「邯郸市DEM.tif.xml」这些文件。很多人以为只有 tif 是数据其余都是“垃圾文件”其实每一个都承担了明确职责。先看一张文件角色对照表文件定位核心作用缺失影响邯郸市DEM.tif主数据存储完整的高程像素矩阵是进行地形分析的基础数据丢失无法恢复邯郸市DEM.tfw世界文件记录左上角坐标、像元大小和旋转参数用于地理配准如果 tif 内嵌坐标缺失影响不大若缺少内嵌地理参考则无法定位邯郸市DEM.tif.ovr金字塔多级降采样加快大范围地图缩放和浏览速度缺失时软件会临时构建首轮打开明显变慢邯郸市DEM.tif.aux.xml辅助元数据GDAL 写入的统计信息、色彩解释和坐标系信息可自动重建不影响像素值邯郸市DEM.tif.xml元数据记录数据来源、生成时间和坐标系统说明属于文档属性不影响运算tif 是一个很灵活的栅格容器地理参考信息可以直接写进文件头。正因如此tfw 这个文件在现代 GIS 软件中经常显得“多余”。但我不建议你随手删掉它——当你以后要把 DEM 转成不带空间信息的普通图片格式或者把 tif 单独拷贝给别人时tfw 会成为唯一能恢复空间位置的依据。ovr 是金字塔文件它本质上是原始栅格的降采样副本软件在缩小时直接读取这些低分辨率版本避免每次重采样都去读原始大文件。如果这个文件丢失QGIS 和 ArcGIS 都会在用户打开数据时临时重建但第一次缩放会卡上几秒到几十秒尤其在邯郸市全境这种范围不小的数据集上体验很明显。2.2 shp 不是单个文件而是一组必须成套出现的侧车文件常见的「邯郸市范围.shp」并不是一个文件而是一族文件shp 存几何形状shx 存几何索引dbf 存属性字段prj 存坐标系定义sbn 和 sbx 是 ArcGIS 的二进制空间索引xml 则描述矢量数据的元信息。真正打开 shapefile 的最低条件是 shp、shx、dbf 三个文件同时存在如果缺失任何一个ArcGIS 会直接拒绝打开。prj 文件的重要性常被忽略它用 WKT 文本描述坐标系统一旦丢失QGIS 可能会把数据默认为 WGS84导致后续所有与 DEM 的叠加都出现偏移。dbf 里存的是属性内容比如区县名称或行政区代码这类字段如果出现乱码往往是 dbf 的编码声明与实际编码不一致最常见的是 GBK 和 UTF-8 之间的错位处理时可以在 QGIS 的「数据源管理」里指定编码或者用iconv工具转换后重新生成 dbf。2.3 用命令行快速检视压缩包和解压结果unzip -l 邯郸市DEM.zip unzip 邯郸市DEM.zip -d ./handan_dem gdalinfo handan_dem/邯郸市DEM.tif | head -n 20第一行只列出 zip 内的文件名不实际解压适合先确认压缩包是否完整第二行解压到指定目录避免当前目录被文件填满第三行调用 GDAL 的 gdalinfo 读取主栅格头信息查看 tif 是否已经内置地理参考。头 20 行通常会包含行列数、像元尺寸、坐标系统和 NoData 值。如果你机器上没有安装 GDAL 命令行工具可以直接用 QGIS 的「图层属性 - 元数据」面板读取同样的信息。我个人的习惯是拿到任何 DEM 数据都先跑一遍 gdalinfo因为后面所有裁剪、投影和可视化都依赖这些头信息的正确性。3. 用 GDAL/QGIS 读取、裁剪和重投影 30 米 DEM3.1 先读懂 gdalinfo 输出里的关键字段拿到数据不要急着拖进软件先用命令确认底层信息gdalinfo -stats 邯郸市DEM.tif重点看四类输出。第一是Size is 宽度, 高度它决定后续计算范围第二是Pixel Size如果坐标系是投影坐标单位是米直接能看到30, -30这样的值如果原始数据是经纬度则会看到接近0.00027, -0.00027这就是 30 米在赤道附近的弧度长度第三是Coordinate System决定是否能和邯郸市范围.shp 正确叠加第四是Band 1下的NoData Value常见值是-9999或-3.40282e38后续统计和高程处理时一定要把这个值排除。-stats参数会让 gdalinfo 顺便计算每个波段的 min、max、mean、stddev如果文件此前没有统计信息这个操作会临时扫描全图数据量大时稍等几秒属正常。3.2 用范围 shp 裁剪 DEM让数据贴合工作区邯郸市 DEM 数据如果覆盖了整个行政区划甚至周边区域直接分析会带来大量无效计算。常见做法是用压缩包里的「邯郸市范围.shp」作为掩膜边界进行精确裁剪gdalwarp -cutline 邯郸市范围.shp -crop_to_cutline -dstnodata -9999 \ 邯郸市DEM.tif 邯郸市_裁剪.tif-cutline指定裁剪用的矢量边界-crop_to_cutline让输出栅格的范围严格贴合边界的外接矩形并栅格化边界外的区域-dstnodata -9999把边界外设成统一的无值标记。如果不加-crop_to_cutlinegdalwarp 会保留输入栅格的完整范围只在裁剪区域外写 NoData结果是文件尺寸没有缩小意义不大。在 QGIS 里对应的工具是「栅格 - 提取 - 按掩膜图层裁剪栅格」对话框中的「目标范围」选择“裁剪到切割线”并在“为输出栅格分配指定的无效值”里填-9999效果与命令一致。这里要特别提醒被漏斗边界裁掉的往往是山区边缘如果你的研究区正好在邯郸市边界线上裁剪后一定要检查边界处的像元是否有连续的高程跳变避免切到一半的山体会在后续坡度计算里产生虚假陡坡。3.3 重投影到适合面积计算的投影坐标系如果 gdalinfo 显示原始坐标系是 WGS84 经纬度那么直接计算坡度、面积、坡长都会因为纬度方向的变形产生误差。我一般会把 DEM 重投影到 CGCS2000 3 度带投影邯郸市大致位于东经 114 度中央经线附近所以选择 EPSG:4547CGCS2000 / 3-degree Gauss-Kruger CM 114Egdalwarp -t_srs EPSG:4547 -r bilinear -tr 30 30 \ 邯郸市DEM.tif 邯郸市_4547.tif-t_srs指定目标坐标系-tr 30 30强制输出像元尺寸为 30 米乘 30 米避免重投影后像元变成非规则尺寸-r bilinear是重采样方法DEM 是高程连续表面双线性插值比最近邻更平滑不容易在陡坡处出现台阶感。这里有一个常见误区很多人看到原始 DEM 已经是 30 米就认为不需要-tr但经纬度数据重投影到投影坐标后每个像元的实际地面尺寸会随纬度变化必须显式指定目标像元尺寸才能保持分辨率一致。重投影一次就够了反复进行 GDAL 重采样会损失原始高程精度尤其是峰值高程。3.4 用 Python 批量读取高程值做质量检查在裁剪和重投影之后我习惯用 rasterio 快速跑一遍质量检查确认 NoData 设置、数值范围和像元数量是否符合预期import rasterio import numpy as np with rasterio.open(邯郸市_4547.tif) as src: dem src.read(1) nodata src.nodata valid dem[dem ! nodata] print(有效像元数:, valid.size) print(高程范围: {:.1f} ~ {:.1f} m.format(valid.min(), valid.max())) print(平均高程: {:.1f} m.format(valid.mean()))读取后先排除 NoData 像元再做统计否则所有统计值都会被-9999污染。有效像元数可以用于估算邯郸市行政区域内被有效覆盖的面积将像元数乘以 900 平方米即可得到大约面积。如果输出的有效像元数明显低于预期通常是 shp 坐标系和 DEM 坐标系不匹配导致掩膜裁剪出来大部分是空值这时回头检查两个数据源的CRS是否一致并优先将 shpto_crs到 DEM 的坐标系。4. 结合范围 shp 的坡度、坡向、山体阴影与等高线分析4.1 30 米 DEM 的适用分析尺度以及和 12.5 米数据的取舍30 米分辨率意味着一个像元覆盖 900 平方米足够支撑邯郸市全域尺度的地表过程模拟。比如主城区扩张的地形适宜性评价、大范围洪水淹没范围估算、流域汇水区划分30 米数据都够用。如果项目要细化到某个县城内部的小沟谷、景观单元或单体建筑周边排水30 米就会显得粗糙这时候很多人会去换成 ALOS 12.5 米数据分辨率确实更高地形细节更丰富但同时数据量增加约 5 倍处理时间随之上升而且 12.5 米包含更多地表噪声不经过滤波直接做坡度分析产生的结果在小范围上反而更碎。我的取舍标准是研究区面积超过几十平方公里时用 30 米 DEM只盯住一个镇或一条沟时换 12.5 米。30 米数据在宏观分析上已经能很好描绘邯郸西部山区和东部平原的过渡特征。4.2 gdaldem 一键生成坡度、坡向和山体阴影接下来是用 DEM 做地形分析最频繁的三个操作gdaldem slope 邯郸市_4547.tif 邯郸市_slope.tif -p gdaldem aspect 邯郸市_4547.tif 邯郸市_aspect.tif gdaldem hillshade -azimuth 315 -altitude 45 \ 邯郸市_4547.tif 邯郸市_hillshade.tifslope计算每个像元的坡度-p表示输出百分比坡度如果不加则输出单位为度的角度值两种单位在色带拉伸时差别很大建议根据后续用途统一aspect输出坡向值范围 0 到 360 度0 表示北方向90 表示东方向负值或 0 表示平地hillshade生成山体阴影栅格-azimuth 315表示光源方位角为西北方向-altitude 45表示太阳高度角 45 度。这三个参数不会改变 DEM 像素值只影响渲染立体感。如果生成的坡向图在平地上出现大量随机噪声说明 DEM 原始数据包含微小起伏噪声常见做法是先对 DEM 做一次低通滤波或者直接接受这个现象因为高地起伏区域的坡向信号远强于噪声。4.3 用范围 shp 叠加制图并统计边界内的高程坡度、坡向只是基础因子真正要落到邯郸市行政边界内做统计时shp 文件就派上用场了。比如计算邯郸市行政区域内的平均坡度、高程分位数可以通过 rasterio.mask 完成import geopandas as gpd import rasterio from rasterio.mask import mask shp gpd.read_file(邯郸市范围.shp) with rasterio.open(邯郸市_4547.tif) as src: if shp.crs ! src.crs: shp shp.to_crs(src.crs) out_image, out_transform mask(src, shp.geometry, cropTrue, nodata-9999) vals out_image[out_image ! -9999] print(区域内高程 min/mean/max: {:.1f} {:.1f} {:.1f}.format( vals.min(), vals.mean(), vals.max())) print(耕地陡坡占比模拟: {:.2f}%.format( ((vals 300) (vals 200)).sum() / vals.size * 100))这里先判断 shp 和 DEM 的坐标系是否一致不一致则转换否则 mask 会直接报错或返回全空。mask的cropTrue裁剪边界范围nodata-9999将边界外设为无值。统计示例里顺带演示了如何用高程区间粗略估算某个海拔带的占比在实际项目里可以替换成坡度栅格做类似计算。注意out_image是一个三维数组读取时用out_image[0]或out_image[out_image ! -9999]都能得到有效值。4.4 从 DEM 生成等高线并导出 shp 和文本等高线是 DEM 最直观的衍生成果。用 GDAL 自带工具可以快速生成gdal_contour -a elev -i 50 邯郸市_4547.tif 邯郸市_contour.shp-a elev指定生成的属性字段名为elev这个字段会存放每条等高线对应的高程值-i 50表示每隔 50 米生成一条。如果遇到陡峭山区50 米间隔会导致等高线过密可以改成-i 100或根据地形起伏动态调整。生成的是 shapefile可以直接叠加到山顶阴影底图上制图。如果需要把等高线属性导出成 txt 给其他程序做进一步分析常见做法是用ogr2ogr转成 CSV 再处理ogr2ogr -f CSV 邯郸市_contour.csv 邯郸市_contour.shp这种方式会保留每条等高线的 elev 属性后续在文本编辑器或 Python 里读取都非常方便。注意 CSV 导出会丢掉几何坐标信息如果还要坐标改用-f GeoJSON或者用 Python 的 geopandas 读取后按需输出。5. 数据可用性验证tfw 一致性、ovr 重建与 shp 完整性5.1 对比 tfw 和 tif 内嵌地理参考找出“漂移”源头拿到数据先做一次坐标一致性检查否则后续所有叠加分析都不可信。把 tfw 内容打印出来再对比 gdalinfo 里的 Origin 和 Pixel Sizecat 邯郸市DEM.tfw gdalinfo 邯郸市DEM.tif | grep -E Origin|Pixel Sizetfw 六行数据分别代表 x 方向像元尺寸、y 旋转项、x 旋转项、y 方向像元尺寸、左上角 x 坐标、左上角 y 坐标。如果 gdalinfo 显示的 Origin 与 tfw 最后两行不一致说明 tif 文件在拷贝或解压过程中被重新写过头信息。这种情况下以 tfw 为准通常更可靠我一般会复制一份 tfw命名为邯郸市DEM.tfw然后重新用gdal_translate -co TFWYES 邯郸市DEM.tif 邯郸市_修正.tif生成一个内嵌坐标修正后的新文件。5.2 重建 ovr解决大文件缩放卡顿如果打开数据时发现缩放明显迟缓或者 QGIS 提示金字塔文件无效可以直接重建 ovrgdaladdo -r average 邯郸市DEM.tif 2 4 8 16-r average指定重采样方法为平均值适合高程这种连续型栅格后面跟着的金字塔层级表示在原始分辨率基础上分别降采样 2、4、8、16 倍。执行后会在 tif 旁边生成一个.ovr文件再次在 QGIS 或 ArcGIS 中加载时缩放预览速度会显著提升。如果处理的是已经重投影过的邯郸市_4547.tif记得对它也执行一次金字塔构建否则后续每次做地形分析都会卡在文件读取上。5.3 检查 shp 文件组是否完整修复属性编码问题最后检查「邯郸市范围.shp」这组文件是否完整并查看其概要信息ogrinfo 邯郸市范围.shp -so -al | head -n 30-so表示只输出概要-al表示列出全部图层。输出里应包含几何类型、要素数量、属性字段列表和坐标系。如果报错提示缺少文件第一反应是检查目录下是否存在 shx 和 dbf。如果缺失可以用ogr2ogr -overwrite -f ESRI Shapefile 邯郸市范围_修复.shp 邯郸市范围.shp强制生成一套完整的新 shapefile。属性字段如果出现中文乱码常见做法是用 QGIS 重新指定源编码为 GBK 或 UTF-8 后另存一份而不是手工去改 dbf 二进制结构后者极易破坏文件。处理 DEM 和 shp 时任何一步先跑一遍 gdalinfo 和 ogrinfo能省下大半排查时间。本文还有配套的精品资源点击获取
返回列表