肇庆30米DEM与shp边界数据:从裁剪到地形因子提取全流程
简介《广东省肇庆市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不是只有SRTMALOS与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 48NEPSG: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。这样得到的tifshp之外没有像元文件体积明显减小。如果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是否等于-9999Size是否明显小于原图。再打开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不然阴影会糊成一片。# 生成山体阴影z2.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, cropTrue, nodata-9999, filledTrue ) dem out_img[0] valid dem[dem ! -9999] if valid.size 0: print(no valid pixel, check crs or nodata) else: print(fvalid_pixels{valid.size}) print(fmean_height_m{valid.mean():.1f}) print(fmin_height_m{valid.min():.1f}) print(fmax_height_m{valid.max():.1f})mask()的作用是用shp里的多边形从DEM中裁剪出一块ndarraycropTrue让输出矩阵范围贴近shpnodata-9999把掩膜外区域填成-9999filledTrue把输入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 GEOMETRYAS_WKTKML文件在Google Earth和部分地图App里能直接打开。若对方只需要坐标或属性文本CSV本身就是最通用的文本表格改扩展名不影响读取。-lco GEOMETRYAS_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同时叠加确认边界区域的等高线没有明显错层数据包就算真正对齐、可以进入业务分析了。三维导出的视角不要用截屏用三维视图自带的导出场景图片功能分辨率按需调高输出结果能直接进汇报材料。本文还有配套的精品资源点击获取