多彩编程 多彩编程MZPH · CODE BLOG
ARTICLE DETAIL

文章详情

深耕前端与后端开发技术的一线实战笔记与踩坑复盘。

上海DEM 12.5m数字高程数据从体检到裁剪的完整处理指南

上海DEM 12.5m数字高程数据从体检到裁剪的完整处理指南 简介上海市12.5米分辨率数字高程模型DEM数据包附带完整市域行政边界矢量文件适用于GIS专业学生、城市规划与地学研究人员进行地形分析、水文模拟、灾害风险评估及三维地形可视化等场景。基于WGS84坐标系兼容性好可直接接入主流GIS平台使用。压缩包共15个文件、约48.05MB以TIFF/ADF栅格数据、SHP/DBF矢量边界及TFW/PRJ/XML辅助文件为主可满足加载、配准与投影转换等需求。已有806人学习下载。12.5米分辨率可清晰呈现城市微地形起伏如黄浦江沿岸、松江丘陵结合行政边界SHP能便捷地裁剪局部区域、按行政区统计高程也可提取等高线、计算坡度坡向、模拟汇水路径为科研与工程项目提供扎实的基础地理数据支撑。1. 上海市DEM数字高程12.5m包里有什么够不够你做项目用别人扔给你一个“上海市DEM数字高程12.5m含区上海市范围shp文件”压缩包里面大概率就是一块全市覆盖的DEM栅格外加一份用来裁剪范围的shp边界文件。12.5m说的是栅格像元大小不是高程精度更不是“连马路牙子都看得清”的精度。拿它做上海市尺度的坡度、坡向、水文分析、填挖方概算和通视分析完全够用拿它报施工图阶段的土方量那就是给自己埋雷。这个数据包真正值钱的地方是DEM和shp已经配套好了省去你去网上找边界、再手动对齐坐标系的时间。适合规划、水利、交通选线、GIS开发的人但动手前一定要先把坐标系、NoData和裁剪方式理清楚否则后面每一步都在补前面的坑。2. 先体检再干活拿到上海DEM的第一件事是查文件、投影和NoData很多人的习惯是解压完直接拖进ArcMap看一眼有地形起伏就开干。我建议你改掉这个习惯因为这类数据包里的坐标系和NoData设置比数据本身更容易翻车。先花10分钟把文件清单、投影信息和无效值查清楚后面能省下半天返工时间。2.1 文件清单tif、tfw、prj、shp各管什么压缩包解开后通常会看到下面几类文件我先说清它们的分工。文件作用DEM_Shanghai_12p5m.tif主数据12.5m分辨率的DEM栅格真正参与分析的高程数据.tfw世界文件和tif同名放一起时提供栅格位置信息.prj坐标系统文本描述了经纬度或投影坐标系shp/shx/dbf边界矢量三件套shp存几何shx做索引dbf存属性这里有个新手最容易犯的错只复制了.shp漏掉.shx和.dbf结果在ArcMap里打不开。三个文件必须放同一个目录。而.tfw文件不是必须的因为GeoTIFF本身就能把空间位置写在tif文件头里但如果压缩包里给了.tfw说明原始栅格可能是早期版本一定要检查坐标系和范围。如果shp文件名是“上海市范围”或“区县边界”打开属性表看一下。有的包把上海外轮廓和区县级边界分两个shp有的只给一个整体边界裁剪时按你的需求选。属性表里如果带“区县名”字段后面就能按区批量裁剪。2.2 在ArcMap/QGIS里加载并检查投影与范围先在ArcMap里把tif和shp都拖进去。右键图层属性切到“源”选项卡重点看两个地方一个是“像元大小”一个是你现在看不到但左上角会显示的“投影坐标系”。像元大小如果显示的是“0.000138888”说明这是经度纬度的度数单位不是米后面不能直接算坡度。如果shp显示的是“CGCS2000_3_Degree_GK_CM_117E”这类投影坐标栅格却是经纬度两者叠加会错位几百米甚至几公里表现为shp跑到海里去。QGIS用户直接在图层右键“属性”信息里看“CRS”和“范围”。另一个技巧是加载后把图层透明度调到50%叠加一个高分辨率影像或路网确认边界没有明显偏移。不要用肉眼在“缩放到图层”模式下判断普通小范围错位肉眼根本看不出来。2.3 用GDAL做无头体检gdalinfo与ogrinfo的判读如果你在服务器上处理或者想留一份可复现的体检记录用GDAL命令行最靠谱。打开终端执行gdalinfo DEM_Shanghai_12p5m.tif重点看这几行Size is 6800, 5200行列数决定数据覆盖范围结合像元大小能大致算出面积。Coordinate System is如果显示的是一长串经纬度要留意后续操作是否需要投影。Pixel Size (12.500000, -12.500000)这才是真正的12.5m分辨率。NoData Value这里写了什么值特别重要常见的是-9999、-32768或空。没有写NoData的DEM不代表没有无效值可能只是没声明。Corner Coordinates四个角点坐标能帮你判断覆盖范围是否真的包含上海全境。然后检查shpogrinfo -so -al 上海市范围.shp这里的-so表示只输出概要信息-al表示列出所有图层。重点看“Extent”和“Layer SRS”“Extent”是shp的外包矩形范围“Layer SRS”是它的坐标系统。我们把两者和刚才的gdalinfo结果放在一起判断shp和栅格范围是否对得上。如果嫌命令行输出不够直观用Python一次性读出来更方便。下面是最小检查脚本from osgeo import gdal, ogr dem_path rD:\Shanghai_DEM\DEM_Shanghai_12p5m.tif shp_path rD:\Shanghai_DEM\上海市范围.shp raster gdal.Open(dem_path, 0) gt raster.GetGeoTransform() print(像元大小X:, gt[1]) print(像元大小Y:, abs(gt[5])) print(栅格宽高:, raster.RasterXSize, raster.RasterYSize) print(投影:, raster.GetProjection()[:200]) band raster.GetRasterBand(1) print(NoData:, band.GetNoDataValue()) vec ogr.Open(shp_path, 0) lyr vec.GetLayer(0) print(shp范围:, lyr.GetExtent())代码里的gt[1]是X方向的像元大小gt[5]是Y方向像元大小通常为负值表示影像从左上角开始向下排列所以取绝对值。投影信息比较长先截前200个字符判断坐标系类型。最后打印的shp范围要和栅格范围做对比如果shp范围明显大于栅格范围说明DEM不一定完全覆盖上海所有区县后面裁剪时要先确认交集。这套体检做完你对这套数据的坐标系、范围和NoData心里就有底了。如果发现shp是投影坐标、栅格是经纬度那不要继续往下走先统一坐标系。统一的方式我建议优先把栅格投影到shp所在的投影坐标系因为DEM重投影到投影坐标后算坡度、算距离都更方便。3. 用上海范围shp裁剪DEM栅格三种做法与参数差异数据体检没问题后下一步就是裁剪。这一步是所有后续分析的地基裁得不对后面坡度、水文、土方全错。常见做法是依靠面图层裁剪DEM栅格但同一个需求在ArcMap里至少有两条路一个是“裁剪”工具一个是“按掩膜提取”工具命令行还有gdalwarp。很多人在这一步就开始纠结我先说清楚它们到底有什么区别。3.1 三种裁剪路径的选型掩膜提取、Clip、gdalwarp怎么选“按掩膜提取”和“裁剪”最大的区别在于NoData处理方式。按掩膜提取会把shp面之外的像元全部置为NoData相当于给DEM戴上一个面具而ArcToolbox里的“裁剪”工具更像用剪刀剪它更关注输出范围是否把外部置为NoData取决于参数设置。实际工作中如果只跑一次且后续要用ArcGIS空间分析模块用“按掩膜提取”顺手如果要批量处理多个区县或者要写进自动化脚本用gdalwarp。我一般按照这个表快速决策方式适用场景关键参数坑ArcMap 按掩膜提取单次交互式裁剪掩膜图层、输出路径边缘像元易出现锯齿ArcMap 裁剪带背景值栅格输入范围、NoData值容易把0值当有效高程gdalwarp批量、自动化-cutline、-crop_to_cutline坐标系不一致时自动重投影可能引入偏移记住一个原则凡是后面还要做坡度、填挖方这些栅格计算裁出来的文件必须明确NoData。否则0值背景会被当成真实海拔0米等高线在边界上密集到吓人。3.2 ArcMap里用面图层掩膜提取DEM参数怎么设ArcMap中操作路径是ArcToolbox选“Spatial Analyst工具”再选“提取分析”双击“按掩膜提取”。打开对话框后填四项输入栅格原始DEM不要填已经做过拉伸显示的图层名。输入栅格或要素掩膜数据上海市范围shp。输出栅格给一个带路径的GeoTIFF名称。可选的环境设置点右侧“环境”把“处理范围”设为“与要素上海市范围相同”把“栅格分析”的“像元大小”设为“与原DEM相同”把“坐标系”定为“与原DEM相同”。这四步里最容易漏的是环境设置。不做坐标系和像元大小设置ArcMap会以掩膜shp的坐标系为基准重采样输出栅格可能从12.5m变成别的分辨率这是一条非常隐蔽的翻车路径。如果你用的是ArcGIS Pro工具箱里同样能找到“按掩膜提取”界面和参数基本一致。操作完成后用识别工具点一下就明白成功与否边界内部应该是正常高程边界外应该是无数据。3.3 命令行批量裁剪gdalwarp -cutline 与 -crop_to_cutline当shp里既有全市范围又有区县边界时命令行就是最好的批量方案。单个要素裁剪的最简命令是gdalwarp -cutline 上海市范围.shp -crop_to_cutline \ -dstnodata -9999 \ -of GTiff \ DEM_Shanghai_12p5m.tif 上海DEM_crop.tif这里几个参数我拆开讲-cutline指定裁剪边界必须是矢量文件。-crop_to_cutline让输出范围跟随要素几何本身而不是shp的外包矩形。不加这个参数裁出来的结果仍然是完整矩形只是边界外的值是空。-dstnodata把输出背景值设为-9999目的是让无效区域在后续计算中自动跳过。-of GTiff输出格式写成GeoTIFF通用性最好。如果shp属性表里有字段能区县比如“区县名”可以在命令中加上-cwhere裁剪一个特定区域gdalwarp -cutline 上海区县.shp -cwhere 区县名浦东新区 \ -crop_to_cutline -dstnodata -9999 \ DEM_Shanghai_12p5m.tif DEM_浦东新区.tif需要注意的是-cwhere后面的SQL语法和字段名大小写由shp原始dbf决定中文环境下经常踩编码坑。我更推荐先把字段问题解决再做批量。如果shp已经包含多个区县面可以用Python循环批量裁剪import os from osgeo import gdal dem_path rD:\Shanghai_DEM\DEM_Shanghai_12p5m.tif shp_path rD:\Shanghai_DEM\上海区县.shp out_dir rD:\Shanghai_DEM\裁剪结果 districts [黄浦区, 徐汇区, 静安区, 浦东新区] os.makedirs(out_dir, exist_okTrue) for name in districts: out_tif os.path.join(out_dir, fDEM_{name}.tif) options gdal.WarpOptions( cutlineDSNameshp_path, cutlineWherefNAME{name}, cropToCutlineTrue, dstNodata-9999, ) gdal.Warp(out_tif, dem_path, optionsoptions) print(完成:, out_tif)代码中的cutlineWhere是按字段过滤字段名一定要改成你shp里实际存在的字段。cropToCutlineTrue对应命令行的-crop_to_cutline。如果内存紧张在WarpOptions里加warpMemoryLimitMB512再配合options[NUM_THREADSALL_CPUS]可以明显提升速度。我通常还会传给gdal.Warp的第三个参数注意这个写法在不同版本GDAL里兼容性略有不同建议把GDAL升级到较新版本。批量裁剪最怕的是某一个区县范围特别小裁出来的栅格只有几十行几十列但12.5m分辨率下这种问题不多。裁剪完成后打开属性表确认每个区县的像元数没有异常少再继续下一步。4. 裁剪完只是第一步把12.5m DEM变成真正能用的地形产品得到一份干净的上海DEM裁剪结果后就可以开始做真正有价值的地形分析了。这一步的目标是把DEM变成坡度、坡向、山体阴影、水文流向这些工程上能直接用的产品。很多人以为裁完就收工其实DEM最大的价值在派生分析。4.1 生成坡度、坡向、山体阴影坐标系和Z因子是两个命门坡度坡向是几乎所有后续分析的基础。ArcGIS里在Spatial Analyst工具中找到“表面分析”点“坡度”输入上一步裁好的DEM输出角度还是百分比看你用途。水利上习惯看百分比交通规划习惯看角度我一般两个都输出一份。这里有个硬前提输入DEM的坐标系必须是米制投影坐标系不能是WGS84经纬度。如果直接拿经纬度DEM算坡度每个像元在东西方向和南北方向对应的地表距离不一样算出来的坡度就没有物理意义。所以如果你在数据体检阶段发现原始DEM是经纬度坐标系请先转成CGCS2000或UTM这类投影坐标。ArcGIS里用“投影栅格”工具或者在“环境”里统一坐标系。如果你偏向GDAL命令行坡度命令是这样的gdaldem slope 上海DEM_crop.tif 坡度.tif -s 1 -p这里的-s 1是Z因子只有高程单位与XY单位一致时才取1。-p表示输出百分比不加的话输出度数。坡向命令把关键字换成aspectgdaldem aspect 上海DEM_crop.tif 坡向.tif山体阴影是用来做效果图和视觉检查的gdaldem hillshade 上海DEM_crop.tif 山体阴影.tif -z 1.0 -alt 45 -az 315-z是垂直拉伸倍数做汇报效果图时可以调到2~3让平地也看得出光影但分析类场景保持1。-alt是太阳高度角-az是方位角。上海这种平缓地形我习惯把太阳高度角调到30~45阴影层次更明显。4.2 填洼与流向提取水文分析的预处理顺序水文分析是DEM最典型的使用场景思路是先填洼再算流向再算汇流累积量。如果直接拿原始DEM做流向那些因数据采集误差产生的伪洼地会截断汇水路径导致水流到一半断掉。上海市区虽然地势平坦但高架桥、基坑、高楼地基都会在DEM上留下一个一个小坑不填洼汇水分析基本没法看。ArcGIS流程水文工具集里先用“填洼”再“流向”最后“汇流累积量”。“填洼”工具一般默认全部填入零限制实际上它默认没有Z限制会把每个洼地填平如果你关心真实洼地蓄水可以在环境里设置Z限制但普通汇水区划分不需要管那么多。流向算法用默认的D8就行12.5m分辨率配合D8是主流组合。QGIS用户可以用SAGA的“填洼”算法它的效果和ArcGIS差不多。注意填洼和流向这两步必须按顺序先填洼后流向不要颠倒也不要跳步。填洼文件生成后再执行“流向”然后把“流向”作为“汇流累积量”的输入三个栅格一定要保持12.5m像元大小任何一步被重采样都会让流向变得极其奇怪。4.3 估算填挖方与剖面用DEM算土方量的公式与误差DEM最容易被高估的能力是土方量计算。两个不同时相的DEM相减就可以用差值栅格估算填挖方。ArcGIS里直接用“填挖方”工具但更可控的是先做差值再自己统计。命令如下gdal_calc.py -A 现状DEM.tif -B 设计DEM.tif \ --outfile差值.tif \ --calcA-B \ --NoDataValue-9999这里-A和-B是输入的两个DEM--calcA-B表示用现状减去设计结果是正的区域就是挖方负的区域就是填方。--NoDataValue-9999让无效区域不参与计算。然后可以用Python快速统计总土方量from osgeo import gdal import numpy as np ds gdal.Open(rD:\Shanghai_DEM\差值.tif) band ds.GetRasterBand(1) arr band.ReadAsArray() nodata band.GetNoDataValue() valid arr[arr ! nodata] cell_area 12.5 * 12.5 cut_volume cell_area * valid[valid 0].sum() fill_volume cell_area * np.abs(valid[valid 0].sum()) print(f挖方约: {cut_volume:.0f} m3) print(f填方约: {fill_volume:.0f} m3)这个脚本的原理是每个像元代表12.5m×12.5m面积用像元面积乘以高差再累加所有挖方像元和填方像元得到总体积。作用是把上千个像元的计算压缩成两行numpy统计。注意这里的结果是一个概算量它的误差主要来自DEM本身12.5m分辨率的像元会把小水沟、小土堆抹平概算可以看趋势最终设计阶段请用实测断面控制。5. 上海12.5m DEM避坑最常见的5个翻车点与排查方法做上海这种高度城市化的地形数据最容易踩的坑不是DEM格式而是那些看起来不显眼、一旦踩中就让你反复返工的问题。我把这几年遇见的几类问题列出来每条都按“现象、原因、解决”的顺序梳理成台账希望你能拿着直接对号入座。5.1 裁剪出来一片“0值海”NoData被当成有效高程现象裁完的DEM在ArcMap里显示整片陆地高程为0边界内起伏的等高线反而全没了最小值统计是0。原因原始DEM的背景NoData在裁剪前就是0或者ArcMap按掩膜提取时把输出背景值设成了0导致0被当成有效高程参与显示与统计。解决回到裁剪环节明确给输出栅格设置一个不影响分析的NoData值比如-9999。已经裁坏了的数据不要盲目用栅格计算器去改因为上海沿海地区真实海拔确实有接近0米的地方把0全部替换掉会误伤真实地形。最好的做法是重新执行gdalwarp -dstnodata -9999再查看波段属性确认NoData值已经写进文件头。5.2 高处是楼顶DEM叠加影像后屋顶比地面还高现象把DEM叠加到卫星影像上发现市中心到处是密集的凸起道路断断续续高架桥下面的地面被抬高或压低。原因这份12.5m DEM很可能来自光学立体像对或雷达数据属于表面模型DSM而非纯地形模型DTM。它记录的是楼顶、树冠、高架桥顶面不是地表。解决城市水文分析、地下管网规划别直接拿DSM当DEM用。除非你能确认原始数据做了严格的地形滤波否则要么换一份经过滤波的公开DEM做交叉对比要么用LiDAR数据在关键区域修正。对于上海这种大型城市很多高阶项目都会选择采购更精确的高程数据来替代公开DEM但这套12.5m数据用于坡度分区、淹没分析、通视分析仍然够看。5.3 坐标系错位shp和tif叠加差了500米现象在ArcMap里同时拖入DEM和上海shp边界对不上一个在陆地上另一个飘在海上。原因这是最常见的坑shp和tif不是同一个坐标系。常见情况是shp是CGCS2000投影坐标DEM是WGS84经纬度或者shp用的是上海地方坐标DEM用了UTM坐标两边对不上。解决先确认谁是谁再做一个统一投影。把DEM投影到shp所在坐标系是更好的选择因为DEM后续算坡度坡向都需要投影坐标。在GDAL里可以直接在gdalwarp时加-t_srs但更稳妥的做法是先查询shp的.prj文件再用ArcGIS“投影栅格”或gdalwarp -t_srs统一坐标。做完后重新叠加确认边界与海岸线吻合。5.4 shp属性表中文乱码没法按区县批量裁剪现象打开shp属性表区县名字段全是“”或者“锟斤拷”在裁剪对话框里写SQL都写不进去。原因shp的dbf属性文件是GBK编码而ArcGIS、QGIS在某些版本里默认按UTF-8读取导致中文乱码。解决QGIS打开shp时在“数据源管理器”选择“GBK”编码ArcMap中乱码不那么常见如果遇到了用“属性表”里的“添加字段”后重新赋值或者用ArcToolbox“复制要素”到新的shapefile让软件重新写dbf编码。命令行用户可以用ogr2ogr指定编码ogr2ogr 上海区县_utf8.shp 上海区县.shp -lco ENCODINGUTF-8转完再重新打开中文字段就正常了。这里要注意乱码问题不解决后面所有按区县名的SQL都抓不到记录。5.5 分辨率不等于精度12.5m能不能用要看项目阶段现象用12.5m DEM做完土方量报告拿去和外业GPS实测断面一比悬殊大到说不过去。原因很多人把“分辨率”理解成“精度”其实12.5m分辨率只说明每隔12.5m采集一个高程值垂直方向上的误差可能是几十厘米甚至米级。特别是在上海这种平坦地区高程变化小垂直误差的影响会被放大。解决做全市尺度的淹没分析、坡度分级、方案比选12.5m数据完全够用做具体地块的场地设计、道路纵断面和工程量结算至少要用激光雷达或1:2000实测地形。用DEM土方量前先找几个明显的已知点做偏差统计如果平均偏差超过0.5m最后结果就只适合当估算参考。6. 把上海DEM用得更顺验证、压缩与批量处理的三个技巧最后一个部分我想给你留三个能直接改善工作流的技巧。第一个技巧是“验证裁剪成果”。很多人在裁剪后只看一眼范围不核对边界结果把边界外的一条狭长NoData带算进了汇水面积。我习惯的做法是在ArcGIS里把DEM和shp同时打开把DEM设为半透明再用“栅格计算器”统计NoData像元比例。如果NoData面积超过边界外残留区域说明裁剪出了缝隙。第二个技巧是输出优化。12.5m分辨率的上海全境DEM如果保存为普通GeoTIFF单文件很可能超过1GB。转成压缩GeoTIFF能明显减小体积gdal_translate 上海DEM_crop.tif 上海DEM_压缩.tif \ -co COMPRESSDEFLATE -co TILEDYES -co BIGTIFFYES这里COMPRESSDEFLATE是无损压缩TILEDYES适合分块读取BIGTIFFYES避免单个文件超过4GB。再加一句添加金字塔gdaladdo 上海DEM_压缩.tif 2 4 8 16这样ArcGIS/QGIS加载时按需读金字塔不再卡顿。第三个技巧是区县批量处理。如果你的shp已经分好了区县面直接用一层循环就能把全部区县裁出来。上面第3章已经写过Python版本这里推荐对大批量文件更好用的bash循环for shp in 区县/*.shp; do name$(basename $shp .shp) gdalwarp -cutline $shp -crop_to_cutline -dstnodata -9999 \ 上海DEM_crop.tif 输出/DEM_${name}.tif done循环原理很简单每轮取一个shp文件名basename去掉后缀得到区县名再把它透传给gdalwarp的-cutline参数。跑之前先加一个空目录测试能不能正常生成第一个文件别直接把整个目录扔进去否则某个shp边界的拓扑错误会让循环惨痛中断。我自己第一次处理这类上海数据时就是没检查坐标系直接裁剪结果边界整整偏了一条街。后来养成了“先gdalinfo、再投影、再裁剪、最后验证”的习惯踩坑几率大幅下降。希望这些经验能帮你在拿到这套DEM时少走弯路把精力放在真正决定分析质量的地方。希望帮到你。本文还有配套的精品资源点击获取
返回列表