汕头30米DEM数据全流程处理：从Shapefile裁剪到Python分析-北京尧图网络

简介汕头市30米分辨率DEM数字高程数据包面向GIS分析、城市规划、灾害评估等场景提供完整的地形栅格与行政边界矢量数据可直接用于坡度计算、可视域分析、径流模拟、环境评估等专业应用。压缩包内共12个文件涵盖TIFF栅格、Shapefile矢量、DBF属性表、投影定义文件等核心类型并附空间索引、坐标参考和金字塔概览文件保障数据在ArcGIS、QGIS等平台中的准确加载与快速浏览整体大小约4.94MB。目前已有456人学习下载。数据以汕头市全域为范围附带.shp边界矢量且因以市界裁切也可能包含周边过渡区域用户无需再自行裁剪和配准配合30米地面分辨率足以支撑中尺度地形分析与工程前期调研尤其适合地理信息相关专业学生、科研人员和城乡规划从业者使用。1. 拿到这份30m汕头DEM先别急着解压第一次看到“广东省汕头市DEM数字高程30m含区域范围shp文件.zip”时大多数人只在意里面的高程栅格直到在 ArcGIS 里打开才发现矩形范围比汕头市辖区大一圈裁剪、统计、出图都得多做一步。这个包好就好在把汕头市行政边界拆成了标准 Shapefile和 DEM 一起交付省掉了从全国矢量库里抠边界的步骤。30 米分辨率意味着每个栅格覆盖约 900 平方米在城市规划和中小流域分析里已经能看出地形起伏趋势又比 12.5 米数据体量小加载和计算都快得多。适合刚接触 DEM 的 GIS 从业者也适合需要快速复现完整地形分析流程的人。2. 拆包之前先认全文件Shapefile 与 DEM 的完整拼图2.1 压缩包里的十二个文件各管什么解压后看到一堆后缀新手容易只盯着汕头市DEM.tif其余文件似乎无关紧要。实际上没有汕头市范围.shp那组文件你连准确的市域边界都要重新找没有.tfw栅格在高斯平面里的位置可能错位一整格。下面这张表把每个文件的角色列清楚。文件角色说明汕头市DEM.tif主栅格30m 分辨率数字高程模型单波段浮点或整型数据汕头市DEM.tfw世界文件与 .tif 同名的纯文本配准文件记录栅格像元在地图坐标下的位置和分辨率汕头市DEM.tif.aux.xml栅格辅助元数据保存 NoData 值、统计信息、坐标系等QGIS 和 GDAL 都会读取汕头市DEM.tif.ovr金字塔文件预先建立的降采样层加速大图缩放显示汕头市DEM.tif.xml栅格元数据ArcGIS 生成内容和 aux.xml 部分重叠汕头市范围.shp边界几何汕头市行政区域面要素汕头市范围.shx空间位置索引配合 shp 快速定位要素汕头市范围.dbf属性表存名称、面积、编码等属性字段汕头市范围.prj投影定义WKT 文本告诉 GIS 这个 shp 用的什么坐标系汕头市范围.sbn / .sbx空间索引ArcGIS 生成用于加速空间查询删除不影响数据本身汕头市范围.shp.xml元数据可选一般不必手工编辑shp实际上不是单个文件而是shp shx dbf的最小组合加上prj才能正确投影。这里还多了sbn/sbx说明这个边界很可能是从 ArcGIS 环境里导出的。给这些文件做备份时最好保留原压缩包因为.sbn丢了以后 ArcGIS 会重新生成但.shp一旦缺了dbf属性表就保不住了。2.2 .tfw、.aux.xml 和 .ovr地理配准的底层约定.tfw是 ESRI 世界文件只有六行文本分别代表 X 方向像元宽度、旋转项、旋转项、Y 方向像元高度、左上角 X、左上角 Y。GDAL 读取 TIF 时会优先使用内嵌的 GeoTIFF 标签如果没有内嵌标签再读.tfw。这个包同时保留.tif和.tfw说明数据可能在某个环节做过一次格式归一化。.tif.aux.xml里的STATISTICS_MINIMUM和STATISTICS_MAXIMUM能提前告诉你高程范围可以先用文本编辑器打开看一眼。.ovr金字塔只影响显示性能不影响原始分析结果如果觉得文件太大删掉后 QGIS 会重新生成。常见做法是保留它因为 30m 的 DEM 覆盖汕头全市后缩放到县级视角时没有金字塔会明显卡顿。2.3 用 gdalinfo 和 ogrinfo 检查数据质量解压后建议做的第一件事不是直接拖进 ArcGIS而是用 GDAL 命令行验证数据完整性。我一般会先进入解压目录执行下面两条命令。cd /path/to/shantou_dem gdalinfo 汕头市DEM.tif ogrinfo -so -al 汕头市范围.shpgdalinfo输出里重点看四段内容Size表示栅格行列数Pixel Size是像元分辨率Coordinate System是投影描述NoData Value用于区分真实高程和无效区。ogrinfo -so -al中的-so是 summary only只汇总不打印全要素-al表示读取全部图层。输出中的Extent就是边界的左上右下坐标和 DEM 的范围比对后能确认边界是否完全落在 DEM 覆盖范围里。如果发现中文文件名在终端里乱码不是文件坏了是字符编码问题。Windows 终端建议先chcp 65001切换 UTF-8或者用英文文件夹名重新解压。GDAL 的 Python 绑定对中文路径的支持比命令行更友好后面第 4 章会用到。3. 用边界 shp 裁剪 30m DEM并生成坡度、坡向与等高线3.1 裁剪前先确认两个数据集的空间参考一致直接用gdalwarp裁剪前最好把上一步拿到的两个Extent摆在一起看。DEM 如果是经纬度坐标Extent 单位是度像元大小大约是 0.000277778shp 如果是高斯投影Extent 单位是米。这两者坐标系不同不能直接叠加。用ogrinfo返回的Geometry: Polygon和GEOGCS判断或者统一交给gdalwarp的-t_srs参数处理。如果 shp 的坐标系和 DEM 不一致常见做法是先给 shp 做投影转换再裁剪。不过gdalwarp本身支持-cutline和-crop_to_cutline它会在裁剪的同时完成重投影所以直接执行下面的命令也能得到正确结果。gdalwarp -cutline 汕头市范围.shp -crop_to_cutline -dstnodata -9999 -overwrite 汕头市DEM.tif shantou_clip.tif-cutline指定边界矢量-crop_to_cutline让输出范围严格按面边界生成而不是矩形包络-dstnodata -9999把边界外的区域统一写成 -9999后续做坡度、三维渲染时不会把空值当 0 米处理。如果不加这个参数GDAL 默认继承原文件的 NoData 设置但一旦原文件没有明确的 NoData边界外很可能出现黑色斑块。-overwrite用于覆盖已有输出文件避免重复运行报错。3.2 坡度与坡向单位选择比算法更关键坡度计算是地形分析最常用的步骤gdaldem提供了现成实现不需要自己写滑动窗口。执行下面的命令gdaldem slope shantou_clip.tif shantou_slope.tif -p -s 111120 -of GTiff gdaldem aspect shantou_clip.tif shantou_aspect.tif -of GTiff-p表示坡度输出为百分比而不是角度如果你后续要按水利规范分级通常用百分比阈值不写-p时输出角度。-s 111120是水平和垂直方向的尺度比只有 DEM 坐标系为度时才需要。它的含义是纬度方向 1 度约等于 111120 米gdaldem slope需要让水平和垂直方向统一到同一单位才能算出真实坡度角。如果 DEM 已经是投影坐标系就不要加-s否则会把真实坡度缩小一半以上。gdaldem aspect输出的坡向以正北为 0 度按顺时针到 360 度。坡向本身没有单位争议但要注意 NoData 区也会参与计算导致边界一圈出现异常坡向。解决办法是在计算前把 -9999 设置为 NoDataGDAL 在-dstnodata -9999时已经写入元数据这里一般不会出错。3.3 等高线生成与 shp 转 KML、TXT等值线是拿 DEM 做图件时的“快速交付物”。gdaldem contour可以直接从裁剪后的栅格生成矢量等高线。gdaldem contour -i 10 shantou_clip.tif shantou_contour.shp-i 10表示等高线间隔 10 米适合汕头这种低海拔地区如果是丘陵区间隔取 20 或 25 米更合适。生成出来的 shp 包含ELEV属性记录每条线的高程值。注意gdaldem contour输出的 shp 会默认带上该字段不需要手工拼接。如果要把等高线分享给外业同事或放入其他 GIS 平台常见做法是转成 KML 或者带表头的 TXT。ogr2ogr可以一条命令完成ogr2ogr -f KML shantou_contour.kml shantou_contour.shp ogr2ogr -f CSV shantou_contour.csv shantou_contour.shpCSV 实际上就是“shp 转 txt”的一种变体属性表被扁平化成行列格式几何字段会丢。如果你需要每个顶点坐标建议用ogr2ogr的-segmentize参数把曲线打散成点后再导出。这里给一个常用对比表。命令适合场景注意事项gdalwarp -cutline按行政区裁剪 DEM先检查两个数据集的坐标系gdaldem slope生成坡度百分比或角度经纬度坐标系需要加 -sgdaldem aspect生成坡向NoData 区会影响边界gdaldem contour抽等高线间隔需结合地形高差ogr2ogr -f KML移动端浏览建议限制字段数量ogr2ogr -f CSV数据交换几何信息会丢失这一套流程走完后你会得到一个按汕头市边界裁剪的 DEM、配套的坡度、坡向和等高线。下一步可以做更有业务味道的分析比如视域、汇水但在此之前我通常先做一次高程统计判断数据里是否存在不合理负值或空洞。4. 用 Python 批量处理 30m DEM高程统计、边界叠加与点位采样4.1 为什么选择 rasterio geopandasGDAL 命令行适合单张处理但当你需要跑几十个站点的节点高程或者把 DEM 和道路 shp、行政区 shp 叠加做批量出图时Python 更合适。rasterio负责读写栅格geopandas负责矢量 io 和投影转换。这两者底层分别是 GDAL 和 Fiona所以和前面命令行看到的数据行为完全一致。汕头市的 DEM 大约几千乘几千像元直接read(1)读入内存没问题。如果是 12.5m 或更高分辨率建议用window参数按块读取避免内存溢出。下面的代码先读取 DEM 基本信息并统计高程范围。import numpy as np import rasterio with rasterio.open(汕头市DEM.tif) as src: elev src.read(1) nodata src.nodata res src.res bounds src.bounds if nodata is not None: elev np.where(elev nodata, np.nan, elev) print(分辨率: {}x{}, 像素数: {}.format(res[0], res[1], elev.shape)) print(高程范围: {:.2f} 到 {:.2f} 米.format(np.nanmin(elev), np.nanmax(elev)))这里src.res返回 (dx, dy)dy 通常是负值但方向不影响分析src.bounds是栅格矩形范围。把nodata替换为np.nan很关键否则最小值会被 -9999 这类占位符污染。如果输出出现 -9999 作为最小值说明 NoData 没有被正确识别。4.2 裁剪到 shp 边界并叠加行政边界图有了 DEM 数组接下来用rasterio.mask做出和gdalwarp一样的结果但保留在 Python 内存里方便继续画图和采样。import geopandas as gpd from rasterio.mask import mask with rasterio.open(汕头市DEM.tif) as src: shp gpd.read_file(汕头市范围.shp, encodingutf-8) if shp.crs.to_string() ! src.crs.to_string(): shp shp.to_crs(src.crs) clip_img, clip_transform mask(src, shp.geometry, cropTrue, nodata-9999) clip_img np.squeeze(clip_img) clip_img np.where(clip_img -9999, np.nan, clip_img)shp.crs.to_string()和src.crs.to_string()比较的是 PROJ 字符串如果一个是 EPSG:4490一个是 EPSG:4618字符不同但都是 CGCS2000 下的经纬度坐标严格判断会触发投影转换。这里多一步转换并不会引入明显误差。mask(src, shp.geometry, cropTrue, nodata-9999)返回裁剪后的数组和新的仿射变换矩阵裁剪结果是从原图行列中抽出的窗口并不丢失空间参考。把边界直接叠加在栅格图上用matplotlib的extent参数把数组行列坐标映射到地理坐标即可。import matplotlib.pyplot as plt from matplotlib import cm fig, ax plt.subplots(figsize(12, 10)) im ax.imshow(clip_img, cmapcm.terrain, extent[clip_transform[2], clip_transform[2] clip_transform[0] * clip_img.shape[1], clip_transform[5] clip_transform[4] * clip_img.shape[0], clip_transform[5]]) gdf shp.boundary.plot(axax, colorblack, linewidth1.2) plt.colorbar(im, labelElevation (m)) plt.savefig(shantou_dem_boundary.png, dpi300)extent四元组顺序是左、右、下、上。注意clip_transform[5]是左上角 Yclip_transform[4]是行方向像元高度二者相乘后再加上首值得到的是右下角 Y。如果直接拿bounds做 extent在旋转栅格里会出错但这里 30m DEM 的旋转项一般为 0所以两种写法结果一致。4.3 从 DEM 提取站点高程shp 转 txt 的常见场景地形分析常遇到“有一批站点坐标要查这些位置的高程”的需求。用rasterio.sample可以在不逐像元读取的情况下完成采样。假设你有一份stations.csv包含 lon 和 lat 两列。import pandas as pd stations pd.read_csv(stations.csv) points list(zip(stations[lon], stations[lat])) with rasterio.open(汕头市DEM.tif) as src: vals [x[0] for x in src.sample(points)] if nodata is not None: vals [np.nan if v nodata else v for v in vals] stations[elev] vals stations.to_csv(stations_elev.txt, sep\t, indexFalse)src.sample(points)需要点坐标和栅格坐标一致如果你的站点是经纬度而 DEM 是高斯投影必须先投影。最稳妥的方式是用geopandas.GeoSeries创建点图层然后to_crs(src.crs)。这段代码导出的是制表符分隔的stations_elev.txt既方便 Excel 打开也和 GIS 里的“shp 转 txt”数据交换习惯一致。4.4 批量处理多幅 DEM 的策略如果后续从多个来源收集了相邻区域的 DEM 瓦片比如从 OpenTopography 下载的 12.5m 数据需要先做图幅拼接再裁剪。GDAL 提供一个vrt中间格式不实际合并文件只建立懒加载的索引视图。gdalbuildvrt mosaic.vrt tile1.tif tile2.tif tile3.tif gdalwarp -cutline 汕头市范围.shp -crop_to_cutline mosaic.vrt shantou_merge_clip.tifmosaic.vrt是 XML 文件记录各瓦片的范围和分辨率。把多张 DEM 拼成 VRT 后再用前面第 3 章的gdalwarp命令裁剪这样不用生成中间大文件节省磁盘空间。需要注意的是不同来源的 DEM 高程基准可能差几米拼接前用gdalinfo -stats检查重叠带的均值差若差值大于 1 米说明两个产品的高程参考不一致不能直接混用。5. 解压与打开 DEM 时的排错zip 损坏、EOCD 缺失和高程异常的判断5.1 zip 完整性测试与 invalid zip archive 修复从网盘或旧硬盘拷出来的 zip 包最常见的问题是下载截断导致解压时报could not find EOCD。这个错误的本质是解压器在文件末尾找不到中央目录记录文件不完整或结构被破坏。先不要急着用各种工具反复点“修复”第一步是做完整性测试。unzip -t 广东省汕头市DEM数字高程30m.zip-t参数逐文件测试 CRC 校验如果输出里有bad CRC或mismatch说明某个文件已经损坏。此时用zip -FF把能恢复的部分抢救出来会损坏内部数据不建议作为常规手段。最常见的可靠做法是重新下载原包下载后先比对文件大小和源页面标注的字节数一致再解压。如果 zip 能解压但 GDAL 打开 TIF 报error read zip archive先确认.tif是否被安全软件拦截或占用了别名为folded的随机盘符。移动文件到纯英文路径下再试一次通常能解决 Windows 长路径或中文路径导致的文件访问异常。5.2 DEM 显示全黑或范围不对时检查 tfw 和 aux.xml裁完的 DEM 在 QGIS 里显示全黑先不要怀疑数据坏了打开图层样式把渲染类型改成“单波段伪彩色”并重新计算最值。如果显示范围偏离真值大概率是.tfw和.tif内嵌的 GeoTIFF 标签不一致。可以用gdalinfo查看输出的 Corner Coordinates如果左上角坐标和 shp 的范围相距很远就手动读取.tfw的前两行和后两行确认像元尺寸和原点。以下命令可以帮助重建丢失的.tfw但前提是你知道栅格左上角坐标和像元大小。gdal_translate -a_srs EPSG:4490 -a_ullr 116.2 23.6 117.2 23.0 汕头市DEM.tif corrected.tif-a_ullr后面依次是左上经度、左上纬度、右下经度、右下纬度这是应急方案不能凭记忆填写必须从辅助元数据或原始来源中获取否则坐标系信息会被错误覆盖。5.3 高程异常检测用gdalinfo -stats或 Python 快速定位高程最小值如果是 -9999、0 或负几万说明 NoData 没有被正确屏蔽。先用gdalinfo -stats 汕头市DEM.tif查看统计信息和STATISTICS_VALID_PERCENT。如果有效像元占比低于 95%说明大量区域是空值边界裁剪、坡度计算都会产生异常。Python 侧的判断更快import numpy as np import rasterio with rasterio.open(汕头市DEM.tif) as src: data src.read(1) valid data[data -1000] print(有效像元: {}%.format(100 * valid.size / data.size)) print(5%-95%分位: {} {}.format(np.percentile(valid, 5), np.percentile(valid, 95)))超过 5% 像元落在 -1000 以下基本可以判定 NoData 掩膜丢失5%-95%分位比最大值最小值更稳健适合发现极端孤立异常点。对异常点位可以配合站点高程 CSV 输出做矩阵定位再用 shp 边界叠加看是否落在汕头市辖区外。这一步做完数据就能放心进入后续的填洼、视域或洪水淹没分析了。本文还有配套的精品资源点击获取

汕头30米DEM数据全流程处理：从Shapefile裁剪到Python分析

相关推荐

CubeSandbox 浏览器沙箱实战：在 MicroVM 中运行无头 Chromium 并用 Playwright CDP 远程驱动

STM32G431嵌入式开发：CubeMX驱动封装与BSP分层实战

ReClip多实例部署指南：如何快速跑多个ReClip隔离负载的完整教程

更多资讯

SpringBoot整合Nacos常见报错排查与解决方案

Comsol EBG能带计算与伪模式处理技术详解

游戏超分辨率替换：OptiScaler 让你在 DLSS、FSR、XeSS 之间任选，帧率与画质兼得

多阈值Otsu的MATLAB实现：从类间方差到递归分割实战

Presidio完整指南：免费实现PII检测与数据脱敏，三步跑通敏感信息匿名化

DataHub DataFlow 与 DataJob 实体管理实战：用 Python SDK 构建数据处理管线元数据

今日更新

基于MATLAB的CRI显色指数计算：从SPD光谱到Ra的完整流程

JSP+Servlet+MySQL博客系统源码部署与优化全攻略

蜂窝调度算法 MATLAB 实现：比例公平、轮询与最大 C/I 对比

本周精选

AI驱动的自动化专利分析平台技术解析

LeetCode-Go 题解：1252. Cells with Odd Values in a Matrix（奇数值单元格计数）双解法剖析

在线答疑系统Java毕设实战：Spring Boot前后端分离与状态流转

本月热门

基于YOLO的焊缝缺陷检测：从数据标注到工业部署全流程实践

分布式服务的决策记录

创业产品的首版范围

读完文章，想聊聊您的网站？