ARTICLE DETAIL

资讯详情

深耕编程入门与网站建设的一线实战洞察。

普洱30m DEM数据处理全流程:从坐标对齐到坡度提取与裁剪避坑

普洱30m DEM数据处理全流程:从坐标对齐到坡度提取与裁剪避坑 简介这份资源面向地理信息、城乡规划、环境研究及遥感分析方向的学习者与从业者提供云南省普洱市30米分辨率的DEM数字高程数据并附带区域行政边界矢量文件可用于地形分析、制图渲染、洪水模拟与空间规划等场景。压缩包共12个文件约211.48MB以tif栅格高程数据为核心配合shp、dbf、prj、shx等Shapefile组件描述区域边界与坐标系统另有ovr、tfw、xml等辅助文件保障影像显示与元数据完整。目前已有527人学习下载。借助这套数据读者可直接在QGIS或ArcGIS中完成地形起伏提取、坡度坡向计算与边界裁剪省去自行搜集与配准的环节是开展普洱市地理分析与制图练习的实用底图素材。1. 普洱30m DEM到手之后先搞清这套数据能干什么、不能干什么拿到“云南省普洱市DEM数字高程数据30m含区域范围shp文件.zip”这类数据包很多人第一反应是解压、拖进ArcGIS然后发现——咦怎么只有一片灰不溜秋的栅格说好的普洱市呢其实这个包的核心就两样东西一份30m分辨率的数字高程栅格通常是tif或img格式一份普洱市行政边界的shp面文件。前者告诉你每一格地面有多高后者告诉你哪些格子属于普洱。30m这个精度意味着每个像素代表地面30米×30米的区域对市级尺度的地形分析、坡度坡向提取、流域划分、选址踏勘来说完全够用但你要是想拿它做茶园梯田的施工放样那就属于用错了工具。这篇笔记就按一线干活的顺序把从数据检查、坐标对齐、按shp裁剪DEM到坡度提取和常见翻车点一步步讲清楚。适合手里已经拿到这套数据、准备做普洱本地地形分析的人。2. 拆开压缩包先别急着拖进软件数据自检与坐标系统一2.1 包内文件到底哪个是主角一个典型的DEM数据包解压后通常包含以下文件。不要看到一堆文件就懵按用途分就三类高程栅格、边界矢量、附带的元数据。文件类型常见扩展名作用是否必须高程栅格.tif / .img / .dem存储每个格网的高程值必须区域范围.shp .shx .dbf .prj普洱市行政边界用于裁剪和统计必须元数据.xml / .txt记录坐标系、精度、来源建议保留金字塔.rrd / .ovr加速显示可重建可删重点看两个东西栅格的坐标系和shp的坐标系。如果两者不一致后面裁剪一定出问题。用ArcGIS的话在Catalog里右键属性看Spatial Reference用QGIS的话右键图层属性看CRS。普洱市常见的坐标系有CGCS2000地理坐标系经纬度单位度和CGCS2000 3度带投影坐标系单位米。30m分辨率的数据如果原始是经纬度实际地面分辨率会随纬度变化在普洱一带大约相当于经度方向30m、纬度方向也接近30m但做面积和距离计算时最好转成投影坐标。2.2 用Python快速读取栅格和矢量基本信息不想开笨重的桌面软件可以用rasterio和geopandas几行代码把底细摸清。下面这段脚本我经常用来做数据入场检查。import rasterio import geopandas as gpd # 读取DEM栅格 dem_path puer_dem_30m.tif with rasterio.open(dem_path) as src: print(栅格尺寸:, src.width, x, src.height) print(波段数:, src.count) print(坐标系:, src.crs) print(地理范围:, src.bounds) print(像元大小:, src.res) # 输出如 (0.000277, 0.000277) 表示度 print(无效值:, src.nodata) # 读取区域范围shp shp_path puer_boundary.shp gdf gpd.read_file(shp_path) print(要素数量:, len(gdf)) print(shp坐标系:, gdf.crs) print(边界范围:, gdf.total_bounds)逻辑说明rasterio.open读取栅格元数据重点看crs和res。如果res是0.000277这种小数说明是地理坐标系单位是度如果是30左右说明是投影坐标系单位是米。geopandas读取shp后看crs是否与栅格一致。参数方面nodata很关键常见是-9999或-32768裁剪和统计时要把它排除否则最小值会变成-9999坡度计算直接崩。2.3 坐标系不一致时的统一策略如果发现栅格是地理坐标系、shp是投影坐标系或者反过来不要硬裁。正确做法是统一到投影坐标系因为30m分辨率做地形分析用米为单位更直观。用QGIS的话右键图层导出时选择目标CRS用ArcGIS的话用Project Raster工具转栅格用Project工具转矢量。转完再确认一次两者的范围是否重叠。普洱市大致在东经99°到102°、北纬22°到24°之间如果转出来范围跑到国外去了说明投影参数选错了。提示转坐标系之前先备份原始文件。投影变换是有损的尤其是栅格重采样双线性插值和最近邻结果不同高程数据建议用双线性或三次卷积。3. 按shp裁剪DEMArcGIS、QGIS和Python三条路都走一遍3.1 ArcGIS里用面图层裁剪栅格的标准操作这是热搜里反复出现的问题在ArcMap中依靠面图层裁剪DEM栅格tif文件。步骤不复杂但参数选错的人不少。第一步打开ArcMap或ArcGIS Pro加载DEM栅格和普洱市边界shp。第二步打开ArcToolbox找到Spatial Analyst Tools → Extraction → Extract by Mask。第三步Input raster选DEMInput raster or feature mask data选shpOutput raster指定输出路径。第四步点环境设置确认Output Coordinates与输入一致Processing Extent选与shp相同或自动。第五步运行。关键参数在Extract by Mask里其实就一个是否勾选“Use Input Features for Clipping Geometry”。勾上输出范围严格按shp边界不勾输出范围是shp的外接矩形。做普洱市地形分析通常要勾上否则你会得到一块矩形区域里面包含相邻州市的高程。3.2 QGIS里用Clip raster by mask layer更顺手QGIS的裁剪工具在Raster → Extraction → Clip Raster by Mask Layer。参数更直白Mask layer选shpSource CRS和Target CRS保持一致勾选“Match the extent of the clipped raster to the extent of the mask layer”。输出格式选GeoTIFF。如果shp有多个面要素比如普洱市下辖各区县可以勾选“Keep resolution of input raster”这样输出分辨率不变。一个容易忽略的点QGIS默认会用掩膜图层的范围裁剪但如果shp的坐标系和栅格不一致工具会报错或输出空白。所以第2章的自检步骤不能省。3.3 Python批量裁剪适合多个区县一次跑完如果普洱市下辖的每个区县都要单独裁一份DEM用arcpy或rasterio写循环最省事。下面用rasterio.mask实现。import rasterio from rasterio.mask import mask import geopandas as gpd import os dem_path puer_dem_30m.tif shp_path puer_counties.shp out_dir output_counties os.makedirs(out_dir, exist_okTrue) gdf gpd.read_file(shp_path) with rasterio.open(dem_path) as src: for idx, row in gdf.iterrows(): geom [row.geometry.__geo_interface__] try: out_image, out_transform mask(src, geom, cropTrue, nodatasrc.nodata) out_meta src.meta.copy() out_meta.update({ height: out_image.shape[1], width: out_image.shape[2], transform: out_transform, nodata: src.nodata }) name row.get(NAME, fcounty_{idx}) out_path os.path.join(out_dir, f{name}_dem.tif) with rasterio.open(out_path, w, **out_meta) as dest: dest.write(out_image) print(f已输出: {out_path}) except Exception as e: print(f第{idx}个要素裁剪失败: {e})逻辑说明mask函数接收栅格数据集和几何对象列表cropTrue表示按几何边界裁剪而不是外接矩形nodata继承源数据。out_meta更新宽高和仿射变换保证输出栅格的地理位置正确。参数方面如果shp字段名不是NAME改成实际字段如果某个区县与DEM范围无交集会抛异常用try包住避免中断。3.4 裁剪结果验证别只看图要看数值裁剪完不要只肉眼看形状对不对。用rasterio读一下统计值最小值、最大值、均值。普洱市海拔范围大致在300米到3300米之间如果裁剪结果最小值是-9999说明nodata没处理好如果最大值超过4000可能混入了其他区域的数据。再检查像元数量普洱市面积约4.5万平方公里30m分辨率下大约5000万个像元裁剪后数量应该在这个量级附近。4. 从DEM到坡度坡向提取地形因子的参数怎么设4.1 坡度提取ArcGIS和Python结果为什么不一样坡度是DEM最常用的衍生因子。ArcGIS的Slope工具在Spatial Analyst → Surface → Slope输入DEM输出单位可选Degree或Percent。Python里用richdem或gdaldem也能算。但很多人发现两个软件算出来的坡度有差异原因通常是算法不同ArcGIS默认用Horn算法三阶反距离平方权差分GDAL的gdaldem slope默认用Zevenbergen Thorne算法。两者在平坦地区差异小在陡坡地区可能差1到2度。如果只是做宏观地形分析这点差异可以接受如果要做工程边坡分级建议统一用同一种算法并在报告里注明。参数上Z因子Z factor很关键。地理坐标系下Z因子要设成纬度对应的换算系数普洱一带大约取1.0左右因为30m分辨率下经纬度与米的换算接近1:1投影坐标系下Z因子保持1。设错了坡度会整体偏大或偏小。4.2 坡向和山体阴影参数少但容易出玄学坡向Aspect输出的是0到360度的方向0代表正北90代表正东。ArcGIS的Aspect工具没有太多参数但要注意平坦区域坡度接近0的坡向会被赋值为-1统计时要排除。山体阴影Hillshade需要设方位角和高度角默认315度和45度适合大多数场景。如果做普洱这种多山地区可以调成270度和30度让山脊和沟谷的立体感更强。4.3 用gdaldem命令行批量生成坡度坡向如果不想开图形界面GDAL的命令行工具非常高效。下面三条命令分别生成坡度、坡向和山体阴影。# 生成坡度单位度Z因子1 gdaldem slope puer_dem_30m.tif puer_slope.tif -p -z 1 -of GTiff # 生成坡向 gdaldem aspect puer_dem_30m.tif puer_aspect.tif -of GTiff # 生成山体阴影方位角315高度角45 gdaldem hillshade puer_dem_30m.tif puer_hillshade.tif -az 315 -alt 45 -of GTiff逻辑说明-p表示输出百分比坡度不加则输出度-z是Z因子-az和-alt控制光照方向。这些命令可以直接写进批处理脚本对多个区县DEM循环执行。注意gdaldem要求输入栅格有正确的nodata值否则边缘会出现异常坡度。5. 避坑与排查普洱DEM处理中最容易翻车的5件事5.1 裁剪后栅格全黑或全白现象Extract by Mask跑完加载输出栅格显示全黑或全白拉伸后也看不到地形。原因最常见的是nodata值设置冲突。源DEM的nodata是-9999裁剪后输出栅格的nodata被自动改成0或其他值导致显示时把0当成有效高程整个栅格看起来就是一片平地。另一个原因是输出栅格的统计值没有重新计算ArcGIS默认用旧统计值渲染。解决在ArcGIS里右键输出栅格 → Data → Export Raster时勾选“Use Renderer”或手动计算统计值Catalog里右键 → Calculate Statistics。用Python的话在out_meta里显式写入nodatasrc.nodata并确保写入后重新打开检查。5.2 shp边界和DEM对不上裁剪出来是空现象运行裁剪工具后报错“ERROR 999999”或输出栅格没有任何像元。原因shp和DEM的坐标系不一致或者shp的范围与DEM完全不重叠。有时候shp看起来在普洱但坐标系被错误定义为WGS84 UTM 47N而DEM是48N差一个带号位置就偏了几百公里。解决先用第2章的脚本打印两者的bounds对比经纬度范围。如果shp的bounds明显偏离普洱用Define Projection重新定义正确的坐标系再用Project做转换。不要直接用Define Projection改投影那只是改标签不改坐标值。5.3 坡度计算结果出现大量0或异常值现象坡度栅格中大片区域值为0或者出现超过90度的值。原因DEM中存在nodata区域但坡度工具没有正确识别nodata把-9999当成高程参与计算导致相邻像元高差巨大坡度异常。另一个原因是Z因子设错地理坐标系下没设Z因子坡度值整体偏小。解决计算坡度前用Set Null或Con工具把nodata区域掩膜掉。在ArcGIS的Slope工具环境设置里把Mask设为DEM的nodata范围。Python里用numpy把nodata替换成np.nan再计算。5.4 裁剪边界出现锯齿或黑边现象按shp裁剪后边界不是平滑的行政界线而是锯齿状或者边界外有一圈黑边。原因栅格是格网结构shp边界是矢量裁剪时按像元中心点是否在面内判断边界像元会被取舍产生锯齿。黑边通常是输出时背景值没设为nodata。解决锯齿是正常现象30m分辨率下无法避免。如果要求边界平滑可以先用shp生成掩膜栅格再做栅格乘法。黑边问题在输出时设置nodata为-9999或0并在显示时设为透明。5.5 文件太大处理速度慢现象普洱市全域30m DEM大约几十MB到上百MB但加上坡度、坡向、山体阴影后文件夹迅速膨胀ArcGIS操作卡顿。原因GeoTIFF默认不压缩每个像元占4字节或8字节。5000万像元的浮点栅格就是200MB到400MB。解决输出时启用压缩。ArcGIS的Export Raster里选LZW或DEFLATE压缩GDAL用-co COMPRESSLZW。如果只做显示可以建金字塔.ovr文件。另外把不需要的中间文件及时清理别让它们堆在工程目录里。6. 进阶技巧用DEM做普洱市地形起伏度分级和流域快速划分6.1 地形起伏度一个窗口统计的实用参数地形起伏度是指单位面积内最高海拔与最低海拔之差常用窗口有3×3、5×5、10×10像元。30m分辨率下5×5窗口相当于150m×150m范围适合分析普洱这种山地丘陵区。在ArcGIS里用Focal Statistics邻域选Rectangle 5×5统计类型选Range输出就是起伏度栅格。然后按自然间断点分级一般分成微起伏30m、小起伏30-70m、中起伏70-150m、大起伏150m。普洱市大部分区域属于小起伏和中起伏局部澜沧江沿岸有大起伏。6.2 用DEM和shp快速划分小流域流域划分通常需要填洼、计算流向、流量累积、提取河网、生成流域边界。ArcGIS的Hydrology工具集可以一键完成但参数多。一个简化做法先用Fill填洼再用Flow Direction计算流向然后用Flow Accumulation计算汇流累积量设定阈值比如1000个像元提取河网最后用Watershed工具以河网节点为出口生成流域。普洱市属于澜沧江水系划分结果可以和实际水系图对比验证。6.3 我踩过的坑和固定习惯早期做普洱DEM分析时我图省事直接用地理坐标系算坡度结果坡度值整体偏小后来才发现Z因子没设。还有一次裁剪完没检查nodata出图时整个区域一片灰被同事笑了半天。现在我的固定习惯是拿到数据先跑一遍第2章的自检脚本裁剪后必看统计值坡度计算前必确认Z因子和nodata。这套流程虽然多花十分钟但省掉了后面反复返工的后悔药。希望帮到你。本文还有配套的精品资源点击获取
返回列表