简介：《广东省肇庆市DEM数字高程30m》是一份面向地理信息学习与研究者的实用数据集，包含肇庆市行政边界范围文件，适合用于地形分析、地表水资源模拟、城市规划辅助、环境研究与灾害风险评估等场景。压缩包内共12个文件，核心是30米分辨率的DEM高程栅格（tif格式），并配套坐标配准文件、投影定义文件、影像金字塔文件，以及行政边界的矢量图形和属性数据库，可同时满足栅格与矢量两类地理数据的教学演示需要；数据在主流地理信息系统软件中可直接加载，能够还原肇庆市及周边区域的三维地表形态。整个资料包约50.22MB，下载便捷，目前已有446人学习浏览；数据覆盖肇庆市全域并向周边适度延伸，可支撑区域尺度的地形对比与专项制图。除高程数据外，还提供带投影信息和元数据的行政边界，初学者可借此理解数字高程模型与Shapefile空间数据的组织方式；研究者也能利用该数据开展坡度坡向计算、径流模拟、环境评估等综合实践，是一份兼顾教学与科研价值的高质量地形数据。
1. 这就是那个“30米高程+边界矢量”一包打天下的地形底图
卫星影像告诉你地表长什么样，DEM告诉你的则是地表每个点的海拔是多少。广东省肇庆市DEM数字高程30m（含区域范围shp文件）.zip 这个包，题面已经把两件关键事说清楚：一个格网间距30米的数字高程模型，外加一份肇庆市范围的shp矢量边界。把栅格和矢量压进同一个压缩包，是地信数据分发中最常见的打包逻辑，好处是拿到手不用再到处找边界文件，解压就能按区域处理。后面无论做坡度坡向、生成山体阴影、按行政区统计海拔，还是出三维地形，起点都在这个DEM上；这份shp则负责把所有计算限制在真正的肇庆范围内，不把周边地区的高程一起算进来。接下来我按拿到这个包后的正常处理顺序展开：先验数据、再对齐坐标系、用shp裁剪，然后批量提取地形因子，最后把结果导成能直接交给别人使用的格式。
2. 先看清包里的两种数据：30米DEM的来头与shp的坐标底细
拿到zip先别急着拖进GIS，解压后先清点文件、确认格式和坐标系。这一步能避免后续至少一半的错位和空值问题，也顺手确认数据包是不是完整可用的。
2.1 30米DEM不是只有SRTM，ALOS与NASADEM也常以30m形式出现
“30米”指的是栅格像元在地面上的边长，一个像元代表30m乘30m的区域。标题只写了分辨率，没写数据源，实际常见的30m高程产品有SRTM、ALOS PALSAR、NASADEM和ASTER GDEM，不同来源在地形细节和空洞处理上差别很大。如果你以前下载过12.5米DEM数据，会明显感觉到30m格网在山谷和坡面转折处更“钝”，对区域尺度分析足够，但在单条冲沟、道路堑坡上会丢失细节。所以拿到包后第一件事，是用工具读出文件头信息：
# 读取DEM文件头信息，重点看坐标系、像元尺寸与NoData gdalinfo zq_dem_30m.tif如果文件名不同，先用ls -R列出解压目录。gdalinfo输出中的Driver告诉你实际格式，Size是行列数，Coordinate System是坐标系，Pixel Size的数值就是像元尺寸。比如Pixel Size = (30, -30)表示x方向30米、y方向向下30米，负号只是栅格坐标系的y轴方向约定，不是反转。顺便记一下NoData Value，后面裁剪统计都要靠它排除无值区。gdalinfo只做只读解析，不会修改文件，可以放心重复执行。
常见来源的特征差异可以用下面这张表快速对照：
| 数据源 | 常见特征 |
|---|---|
| SRTM 1弧秒 | 全球覆盖，空洞多分布在高山积雪区 |
| ALOS PALSAR 12.5m/30m | 细节更丰富，文件体积偏大 |
| NASADEM | 修复SRTM空洞，可作为替代底图 |
| ASTER GDEM | 全球30m，水域云雾区可能出现异常值 |
2.2 shp文件不是单个文件，是“4件套”或更多
shapefile在zip包里通常不是一个孤立文件。能正常打开的shp必须带.shx和.dbf，最好还带.prj。.shp只存几何坐标，.shx提供几何索引，.dbf存属性记录，.prj说明坐标系。缺少.dbf时GIS会报属性表无法打开，缺少.prj时图层能显示，但坐标系未知，后面叠加结果不可信。收到数据包后先检查同一前缀下的文件是否齐全，再顺手用ogrinfo读一下图层概要：
# 只输出图层概要，验证shp能否被正常解析 ogrinfo -so -al zhaoqing_boundary.shp-so是summary only，只输出概要不展开要素；-al表示处理所有图层。输出里能看到图层类型是Polygon、属性字段列表和坐标系统信息。这里的属性字段往往比你想像的要多，比如“市”“区县”“面积”等，后续做按区域统计时要直接引用字段名。道路shp、行政区划shp这类矢量专题也遵循同样结构，换一张图纸，套路不变。
2.3 坐标系的坑：先统一再裁剪
DEM和shp最常见的错位原因是坐标系不一致。shp边界为了通用，常用WGS84或CGCS2000经纬度，而DEM为了保持地面分辨率，常被重投影成UTM或高斯平面坐标。判断方法是把两个图层在GIS里叠加，预览中位置差几公里但属性都正常，那基本就是坐标系不同。以肇庆为例，它位于北纬23度、东经111度到112.5度附近，适合用WGS84 / UTM zone 48N（EPSG:32648）作为统一平面坐标。
# 将DEM重投影到UTM 48N，双线性插值避免坡度台阶 gdalwarp -t_srs EPSG:32648 -r bilinear zq_dem_30m.tif zq_dem_32648.tif-t_srs指定目标坐标系，-r bilinear表示用双线性重采样。高程数据重采样不适合用最近邻，否则山坡上会出现明显的台阶状伪影。DEM转完后再用同一EPSG转shp，两个图层就真正叠上了：
# 把shp也转到UTM 48N，属性字段会原样保留 ogr2ogr -t_srs EPSG:32648 zhaoqing_bnd_32648.shp zhaoqing_boundary.shpogr2ogr默认会保留属性字段和要素几何。转换成UTM后，后续裁剪、面积统计都能直接用米作单位，不必再把经纬度换算一次。如果shp本身是CGCS2000基准，与WGS84在30米网格级别上的差异可忽略，但涉及国土成果交付时仍要以项目要求的基准为准。
3. 把DEM和shp叠在一起：解压、配准、按shp裁剪
进入操作阶段。这里按实际处理数据包的顺序写，兼顾QGIS、ArcGIS和命令行三种场景。三者的逻辑完全一致，只是入口不同。
3.1 zip解压的细节：目录结构与中文文件名
这个zip文件名是带括号和中文的完整命名，属于典型下载包。在Windows双击解压通常没问题，但在Linux或macOS终端里，中文引号和括号会被shell解析成特殊字符，需要给文件名加双引号。解压路径最好也改成纯英文：
# 创建英文目录，避免中文路径带来的工具链问题 mkdir -p /data/gis/zhaoqing # 解压zip到指定目录 unzip "/data/gis/download/广东省肇庆市DEM数字高程30m（含区域范围shp文件）.zip" -d /data/gis/zhaoqing # 列出解压结果，检查shp四件套是否齐全 ls -lh /data/gis/zhaoqing-d指定解压目标目录。解压时常见问题是.shx、.dbf等文件被邮件系统或安全软件判定为“未知附件”而没放进zip，结果只剩一个.shp，这种情况下GIS打不开，报错还比较隐晦。如果解压后看到tif和shp的几个组成部分都在，就可以进行下一步。若文件名在Linux下乱码，使用支持编码转换的7-Zip或unzip的-O参数重新解压，别在读完文件后再批量改名，容易连带丢投影元数据。
3.2 用gdalwarp对齐坐标系，用ogr2ogr转shp
数据包里的DEM和shp如果原坐标系是经纬度，先把DEM重投影到shp的坐标系，或者反过来，但要保持统一。以EPSG:32648为例，前面已经给出gdalwarp命令。需要特别说明的是三个高频参数：
| 参数 | 作用 | 说明 |
|---|---|---|
-t_srs | 目标坐标系 | 可写EPSG代码或Proj4字符串 |
-r | 重采样方法 | DEM推荐bilinear或cubic，避免nearest |
-srcnodata | 原始数据空值 | 不指定可能把nodata重采样成边缘灰值 |
另外，重投影不是越多越好。DEM每重投影一次，像元值会经过一次插值，坡度、坡向这些衍生数据都会引入人为噪声。如果shp没有特别复杂的投影需求，尽量让DEM向shp的坐标靠拢，而不是反过来为了“保持原始DEM”把shp做成地理坐标，再用经纬度做面积统计。
3.3 用shp把DEM裁剪出来：cutline与crop_to_cutline
范围shp最典型的用途就是做不规则裁剪。直接在GIS里用常规“裁剪”往往把DEM裁成矩形，边界外的像素只是被设成0或nodata，文件不但没变小，后续统计还要多做一次排除。正确做法是让输出栅格边界严格贴合shp几何。命令行做法是：
# 用shp边界做不规则裁剪，输出范围与shp严格一致 gdalwarp \ -cutline zhaoqing_bnd_32648.shp \ -crop_to_cutline \ -dstnodata -9999 \ zq_dem_32648.tif zq_dem_clipped.tif参数作用：-cutline指向shp文件；-crop_to_cutline让输出行列数按shp几何范围计算；-dstnodata -9999指定输出无值区为-9999。这样得到的tif，shp之外没有像元，文件体积明显减小。如果shp里有多个辖区，比如镇街边界，可以用-cl指定图层，或配合-cwhere加SQL条件只保留某条记录，生成对应的单独文件。QGIS里对应的入口是“栅格→提取→按掩膜图层裁剪”，ArcGIS里是Extract by Mask或Clip工具并勾选“使用输入要素裁剪几何”，底层思路和上面的cutline一致。
没用shp时，矩形裁剪可以用gdal_translate -projwin，四个数字是左上角x、y和右下角x、y。要注意projwin按像元边界计算，与shp范围存在不到一个像元的误差，要求高精度对齐时不要混用。
3.4 裁剪后的检查：nodata、空洞和直方图
裁剪不是一锤子买卖。用gdalinfo再看一遍裁剪输出：
# 检查输出tif的NoData与尺寸是否正常 gdalinfo zq_dem_clipped.tif重点看NoData Value是否等于-9999，Size是否明显小于原图。再打开QGIS的直方图面板，如果-9999处有极高柱状，说明裁剪有效但图像仍含少量无值像素；如果直方图在0处出现尖峰，多半是原始DEM用0填充海平面以下，或裁剪时把背景0值留了下来。此时用gdal_translate补一道：
# 强制把输出NoData标记改为-9999 gdal_translate -a_nodata -9999 zq_dem_clipped.tif zq_dem_clean.tif-a_nodata只改元数据，不重算像素值，执行很快。若要把原值为0的像素统一改成nodata，得用gdal_calc.py或QGIS栅格计算器。看起来多一步，但总比追着一堆异常高点检查数小时强。
4. 从30米DEM里生成坡度、坡向与山体阴影
DEM原始高程只是半成品，大部分项目要的是地形因子。GDAL自带gdaldem工具，30m DEM配上shp边界，可以快速派生山体阴影、坡度和坡向。
4.1 山体阴影：方位角与垂直拉伸是两个核心参数
先做最能直观检查DEM质量的产物——山体阴影，它用一个假想光源模拟光照，输出0到255灰度图。显示时的立体感受两个参数控制，一个是太阳方位角-az，一个是太阳高度角-alt。传统地图常用315度方向光、45度高度角，立体感均衡；地形南北走向明显时，把方位角调到225度或135度可以强化侧向沟壑。垂直拉伸-z控制高差被放大的倍数，丘陵地带用2.0到2.5，山区则用1.0到1.5，不然阴影会糊成一片。
# 生成山体阴影，z=2.2增强起伏感 gdaldem hillshade zq_dem_clipped.tif zq_hillshade.tif -z 2.2 -az 315 -alt 45-z 2.2意味着高程被乘以2.2再计算阴影，这对海拔起伏不大的区域尤其有用。如果输出的山体阴影在shp边缘出现明显黑色三角，优先检查裁剪后的DEM是否残留nodata，而不是改光源参数。
4.2 坡度计算：经纬度DEM必须加scale参数
坡度是相邻像元高差除以水平距离。当DEM是经纬度坐标系时，x、y方向单位是度，z方向单位是米，直接算出来的坡度会严重失真。GDAL提供-s参数，把水平距离从度换算成米，最常用的经验值是111120，即一纬度约111.12公里：
# 经纬度DEM算坡度，必须加scale，-s把度转成米 gdaldem slope zq_dem_wgs84.tif zq_slope_deg.tif -s 111120如果已经用gdalwarp把DEM转到UTM上，水平单位就是米，直接省略-s；加了反而会把30米当成30度来算。想要百分比坡度，再加上-p，不写时输出的是度数坡度。两者换算关系是 degrees = atan(percent/100)。30m DEM在平缓的珠三角边缘地带容易出现大量0到1度的“湖面”，这是格网分辨率决定的，不代表地面绝对平整。坡度分级阈值在不同行业不统一，下表是常用可视化参考：
| 坡度（度） | 常见描述 |
|---|---|
| 0-2 | 平地 |
| 2-6 | 缓坡 |
| 6-15 | 斜坡 |
| 15-25 | 陡坡 |
| >25 | 急坡 |
做水土保持等专题时，要按项目规范替换阈值，不要直接拿这张表当分析标准。
4.3 用shp批量统计区域内的高程与坡度：Python最小脚本
实际业务中，光有坡度和高程分布还不够，通常要汇报肇庆全市的“平均海拔”“最大坡度”。如果手头有Python环境，用rasterio和geopandas就能完成。脚本核心是让shp与DEM保持统一坐标系，然后让shp的几何作为掩膜从DEM中提取像元：
import geopandas as gpd import rasterio from rasterio.mask import mask DEM_TIF = "zq_dem_clipped.tif" BOUNDARY_SHP = "zhaoqing_bnd_32648.shp" with rasterio.open(DEM_TIF) as src: gdf = gpd.read_file(BOUNDARY_SHP) geoms = [g for g in gdf.geometry if not g.is_empty] out_img, out_transform = mask( src, geoms, crop=True, nodata=-9999, filled=True ) dem = out_img[0] valid = dem[dem != -9999] if valid.size == 0: print("no valid pixel, check crs or nodata") else: print(f"valid_pixels={valid.size}") print(f"mean_height_m={valid.mean():.1f}") print(f"min_height_m={valid.min():.1f}") print(f"max_height_m={valid.max():.1f}")mask()的作用是用shp里的多边形从DEM中裁剪出一块ndarray，crop=True让输出矩阵范围贴近shp，nodata=-9999把掩膜外区域填成-9999，filled=True把输入nodata先填成-9999再参与后续判断。如果把同样流程换成slope.tif，就能统计平均坡度，脚本不用动结构。这里的几何字段读取依赖geopandas自动匹配.shp属性，注意不要在多进程循环里反复打开同一个rasterio文件句柄，容易触发GDAL文件锁问题。
4.4 分块批处理：镇街边界与渔网分割shp
需要按镇街统计时，遍历shp里的每一行要素，用同样的mask逻辑生成多个子tif，再逐个统计。如果shp里没有镇街，而是想按规则网格分块处理，可以先在QGIS里用“创建渔网”工具生成1km乘1km的格网shp，再把每个格网单独导出成shp。得到grid_0.shp、grid_1.shp等一组文件后，用循环逐个裁剪：
# 对每个格网shp执行一次cutline裁剪，输出单独tif for grid in grid_*.shp; do gdalwarp -cutline "$grid" -crop_to_cutline \ zq_dem_clipped.tif "${grid%.shp}_dem.tif" done这里有一个容易踩的坑：gdalwarp的-cutline一次命令只处理一个shp的所有要素，不会按要素自动拆分。如果只有一个shp包含很多网格要素，必须先用ogr2ogr -where或Python按FID拆成单要素shp，再进循环。脚本里shp文件名要保证没有空格，否则for循环会把词拆开，这也是批处理脚本报错的常见原因。
5. 收尾技巧：完整性自检、范围导出与三维展示
到了最后一步，数据已经能用，但交付和复盘时还有几个顺手的小习惯。
5.1 用unzip -t验证zip完整性
数据包很大，网盘下载容易在中途丢字节。解压前先跑一下完整性测试：
# 只校验压缩包，不释放文件 unzip -t "广东省肇庆市DEM数字高程30m（含区域范围shp文件）.zip"这个命令逐个读取压缩条目的CRC校验值，输出末尾出现No errors detected in compressed data才算通过。如果报CRC failed，或提示End-of-central-directory signature not found，不要强行解压继续用，DEM某个局部区域的像元可能损坏，在坡度图上表现为规则的细碎高值点，回看原始数据时反而看不出问题。重新下载后再检测，比直接解压更省时间。
5.2 把shp导出成txt或kml，交给没装GIS的人看范围
区域范围shp经常要发给外业或被下游脚本消费。GDAL提供现成的KML和CSV导出：
# 导出KML，直接用Google Earth打开 ogr2ogr -f KML zhaoqing_outline.kml zhaoqing_boundary.shp # 导出带WKT几何的CSV，扩展名改成txt即可 ogr2ogr -f CSV zhaoqing_boundary.csv zhaoqing_boundary.shp -lco GEOMETRY=AS_WKTKML文件在Google Earth和部分地图App里能直接打开。若对方只需要坐标或属性文本，CSV本身就是最通用的文本表格，改扩展名不影响读取。-lco GEOMETRY=AS_WKT让输出CSV包含一列WKT字符串，字段顺序受属性顺序影响，因此读取txt时不要用cut对列编号做硬编码。想批量转kml，就用同样命令循环遍历目录下的所有shp文件。
5.3 在QGIS里快速查看三维效果
30米DEM和shp已经叠加后，打开QGIS的3D Map View，在场景选项里把Elevation选为DEM图层，再把山体阴影作为覆盖层叠加到DEM上，垂直比例调成1到3，预览窗口里就能看到肇庆山体的大致起伏。检查山脊线走向是否与shp边界重合，或者把DEM的等高线与shp同时叠加，确认边界区域的等高线没有明显错层，数据包就算真正对齐、可以进入业务分析了。三维导出的视角不要用截屏，用三维视图自带的导出场景图片功能，分辨率按需调高，输出结果能直接进汇报材料。
本文还有配套的精品资源，点击获取
