ARTICLE DETAIL

资讯详情

深耕网站建设与运营推广的一线实战洞察。

30m DEM与shp叠加分析:从坐标体检到坡度、等高线提取的完整指南

30m DEM与shp叠加分析:从坐标体检到坡度、等高线提取的完整指南 简介江西省赣州市三十米分辨率数字高程模型DEM及市级范围矢量边界为GIS使用者、城乡规划师、环境科研人员及应急管理人员提供可直接分析的地形基础数据。栅格文件记录逐点海拔可提取高程、坡度、坡向、地形起伏度等指标配套的SHP市界范围便于快速裁剪、统计与制图有效规避跨边界数据干扰。压缩包共十二个文件约一百四十MB包含TIF主数据、SHP矢量、DBF属性、PRJ投影、TFW配准、空间索引及XML元数据体系完整。该数据集已有1108人学习下载在ArcGIS、QGIS中加载即可用于汇水区划分、可视域分析、山体阴影模拟等场景。结合人口、土地利用等图层还可支撑水利设施选址、交通选线、生态敏感区识别与洪涝风险制图为区域规划与决策提供科学依据。1. 江西省赣州市30m DEM这个zip包到底能拿来干什么拿到“江西省赣州市DEM数字高程数据30m含本市级范围shp文件.zip”先别急着解压拖进ArcGIS。这个包的核心是一份覆盖赣州市范围的30米分辨率数字高程模型栅格再配一份赣州市市级行政范围shp矢量边界。30米分辨率意味着每个像元对应地面约30米×30米的方格能支撑市级尺度的坡度坡向、洪水淹没、选址、土方量估算和流域分析但别指望它看清一条几米宽的冲沟。适合谁国土空间规划、水利、林业、交通选线和GIS教学都能用新手想要快速出图熟手需要的是用shp把DEM裁出干净范围、算清参数。2. 解压后的第一件事坐标系体检与数据验收拿到zip先解压。这类数据包的常见结构是一个tif或img格式的DEM栅格一个以“市级范围”或“市界”命名的shp面文件还可能附带元数据说明。shp不是数据本身它是限定分析范围的“框”。很多项目翻车翻在没做坐标系体检就开始裁切、算坡度最后结果错得离谱。2.1 zip里通常装了哪些文件先看文件清单。DEM栅格如果是tif格式通常是单波段整型或浮点型栅格有的包会附tfw坐标文件shp则是一组同名不同扩展名的文件。打开shp之前先确认你要用的那个面图层到底是“市级边界”还是“市级范围”两者在语义上有细微差别直接影响裁剪范围。文件后缀作用常见问题.shp面矢量几何缺.prj时坐标系“裸奔”.shx几何索引缺失时图层打不开.dbf属性表字段名乱码常见于旧成果.prj坐标系定义与dem坐标系不一致.tif/.imgDEM栅格无数据值需单独确认属性表里的字段也要扫一眼。市级范围shp通常有行政区名称、面积、周长之类字段如果出现“省界线文件”场景注意省1、省2这类字段是省界类型编码不是市界别拿错字段做筛选。假如这份shp是某次外业从dwg转出来的属性表里很可能没有面积字段需要后期补充计算。2.2 坐标系体检shp和DEM不在一个坐标系会怎样常见的坑是dem文件是地理坐标系比如WGS84经纬度像元尺寸单位是度shp却是CGCS2000高斯投影单位是米。直接叠加时矢量会整体“找不着北”偏移几十到几百米。即便两者椭球基准一致经纬度与投影坐标混用后续坡度、坡向、面积、长度都会算错。我习惯用GDAL做快速体检。下面这段Python能同时读栅格和矢量信息适合刚到手的数据包from osgeo import gdal, ogr def inspect_raster(path): ds gdal.Open(path) if ds is None: print(raster open failed:, path) return gt ds.GetGeoTransform() srs ds.GetProjection() print(size:, ds.RasterXSize, x, ds.RasterYSize) print(pixel size:, gt[1], gt[5]) print(origin:, gt[0], gt[3]) print(srs:, srs[:200]) for i in range(1, ds.RasterCount 1): band ds.GetRasterBand(i) nd band.GetNoDataValue() stats band.ComputeStatistics(0) print(fband {i}: nodata{nd}, min{stats[0]}, max{stats[1]}) def inspect_vector(path): lyr ogr.Open(path).GetLayer(0) print(feature count:, lyr.GetFeatureCount()) print(extent:, lyr.GetExtent()) srs lyr.GetSpatialRef() print(srs:, srs.ExportToWkt()[:200] if srs else none) inspect_raster(赣州市dem.tif) inspect_vector(赣州市界.shp)gdal.Open读的是栅格头信息。GetGeoTransform返回六个数字依次是左上角X坐标、像元宽、旋转项、左上角Y坐标、旋转项、像元高。像元宽为正、像元高为负这正常表示数据从上往下排列。GetNoDataValue一定要看很多DEM把高程空值写成-32768ComputeStatistics(0)会算全图最小值和最大值如果数据里混入几万个-32768最小值一下就会露馅。参数上重点关注三点一是像元尺寸单位是度还是米二是SRS字符串里出现的是GCS_WGS_1984、CGCS2000还是别的关键词三是波段NoData值是否处于一个反常的数值区间。这三点确认正常再继续往下做。2.3 在ArcGIS里完成第一次加载与验收坐标系体检过后再把tif拖进ArcMap或ArcGIS Pro。加载顺序建议是先加shp再加dem先看shp属性表字段是否有名称字段再用Identify工具在dem上点几个点对照看高程值是否合理。赣州的地势有起伏几十米到一千多米都正常如果点出来出现负几万多半是NoData没处理干净。还有一个容易误会的点ArcGIS里“把tif转成dem文件”即另存为.dem后缀只是换了文件外壳像元值没有变化。真正影响分析的是坐标系、分辨率、NoData三个参数不是后缀。如果shp缺.prj文件图层会以未知坐标系加载和dem叠加时提示坐标系不一致。常见的补救做法是先看随包说明文档确认后给shp补一个.prj没有依据时不要乱选坐标系否则后面裁出来的范围会整体偏移。提示验收标准其实就三条——范围有交集、坐标系能对齐、NoData值明确。三条都满足这份数据才谈得上可用。3. 用市级shp裁剪出赣州市域的DEM两种可行路径及参数设置数据体检完就要“削足适履”了。原始DEM可能覆盖好几个图幅或超出赣州市界不少边角。用市级范围shp把DEM裁出来既是出图需要也是让后续统计不掺水分的前提。按掩膜提取与gdalwarp两条路径我都在用前者交互直观后者批处理省事。3.1 路径AArcGIS按掩膜提取的参数设置ArcGIS里对应的工具叫“按掩膜提取”英文Extract by Mask归属于Spatial Analyst工具。在ArcToolbox里搜索“按掩膜提取”输入栅格选dem输入掩膜数据选市级边界shp输出路径写清楚。点开环境设置把“像元大小”设为与输入栅格一致“捕捉栅格”也选dem这一步是避免输出栅格在对齐上悄悄跑偏。这里有个高频误用有人用“提取分析-裁剪”工具裁剪shp时只给外接矩形范围没有真正贴合多边形边界结果边缘出现大量NoData或包含市界以外的区域。按掩膜提取会逐像元判断面内面外多边形的凹凸都能贴合。输出前记得把NoData值保持默认不要随手填0否则后续坡度计算会把空白区当作平地。另外注意数据类型的差异dem是浮点型高程时输出栅格默认也是浮点dem是整型时输出会变成整型坡度分析时出现台阶状渲染这不是算法错了是高程精度被截断了。遇到这种情况先在栅格计算器里把dem转成浮点再分析。3.2 路径BGDAL命令行一条命令搞定不用ArcGIS时我更常用gdalwarp。它直接读shp做精确裁剪一条命令就能复现处理几百个文件时特别省事。# 用市级边界shp精确裁剪DEM输出tif gdalwarp -cutline 赣州市界.shp -crop_to_cutline -overwrite \ -tr 30 30 -r nearest -dstnodata -32768 \ 赣州市dem.tif 赣州裁剪.tif逻辑说明-cutline指定矢量边界shp面内保留、面外抛弃-crop_to_cutline表示输出范围贴合边界而不是外接矩形-tr 30 30是输出像元尺寸这里要注意如果源数据是经纬度坐标30代表的是30度而不是30米此时要先把度换算成约0.0002777度不能直接写30-dstnodata -32768把裁剪后生成的NoData统一写成-32768方便后期识别。推荐-r nearest。坡度、汇水分析对高程原始值敏感双线性或三次卷积会在边界和坡面产生平滑伪高程看起来好看实际每改一次值就引入一次误差。如果只做显示用的彩色影像双线性可以做定量分析一律最邻近。这条规则我很少破例。如果原始dem是分幅的先做镶嵌再裁剪。常见做法是gdal_merge或ArcGIS镶嵌工具输出一个完整覆盖赣州的tif再用上面的命令裁。直接对每一幅分别裁剪再拼接边界处容易出现重叠或缝隙后续统计各算各的很麻烦。3.3 裁剪结果的三查方法裁剪完成执行下面这条命令快速验收gdalinfo 赣州裁剪.tif | grep -E Size|Pixel Size|Origin|NoDatagrep能直接带出关键行。看到Size和Pixel Size时对照输入dem的尺寸如果像元尺寸变了说明重采样参数没设对面积比例异常多半是坐标系选错或shp范围有问题。此时回到shp的extent和dem的extent检查两者是否有交集如果根本没交集大概率是两套坐标系基准不一致。另一个验证办法是在ArcGIS里用“按属性选择”选中shp要素再做“缩放至所选要素”看dem能否正好覆盖它。如果dem边缘露出说明裁剪范围与shp之间有缝隙需要回到坐标系问题排查。若dem与shp只有部分相交先相交分析看重叠范围是否足够重叠太小的话这次裁剪结果不可靠别硬用。注意裁剪后边缘一圈黑边是NoData区不是数据损坏。后续所有工具只要在处理前声明NoData就不会把黑边算进高程统计。4. 从30m DEM产出等高线、坡度与流域成果参数选择与分析边界DEM不是拿来看的“彩图”它的价值在于派生数据。等高线、坡度、坡向、流向、汇水区是规划和水利最常用的几个产品。下面按我习惯的顺序展开参数一并说明。4.1 等高线生成等高距该取多少等高线生成在ArcGIS里叫“等值线”工具输入dem输出矢量线。关键参数是等高线间距。30m分辨率数据平缓丘陵取5m会出现大量细碎弯折实际项目里更多用10m或20m。我的习惯是先看dem的高程范围比如200到1400米按“范围除以20”估算初始间距再根据出图比例尺微调总图用20m局部大比例尺用10m特殊地段再加密。命令行方式更可控gdal_contour一次就能输出为GeoPackage省去shp多文件的麻烦# 每10米一条等高线高程值存进ELEV字段 gdal_contour -a ELEV -i 10.0 -f GPKG 赣州裁剪.tif 赣州等高线.gpkg参数含义-a ELEV指高程值写入的字段名-i 10.0是等高距单位与dem高程单位一致-f GPKG指定输出GeoPackage避免shp字段名截断。生成后到ArcGIS里检查河谷区域有没有等高线交叉、平地处是否过密。交叉线通常出现在陡崖区域出现时先回填洼再试一次。与这个方向相反的“用等高线生成DEM”也常用来补空白区ArcGIS里对应Topo to Raster工具。注意等高线转DEM时必须把河流、边界作为约束条件否则山谷里会出现平行条带伪影。如果30m DEM局部有空洞用等高线内插补洞可行但补出来的精度不要当作原厂同档看待。4.2 坡度坡向与日照分析30m分辨率的适用边界坡度工具输入dem输出单位可选度和百分比默认是度。它背后的计算窗口是3×3像元也就是约90米×90米的邻域。这意味着30m DEM根本表达不了十几米尺度的微地形坡度结果天然偏平滑。做地灾或边坡分析时直接拿这个坡度值做阈值判定要慎重最好先和现场实测对比再定门槛。另一个高频错误是Z因子。当dem是投影坐标、单位是米时Z因子设1当dem还是经纬度坐标时X和Y单位是度高程单位是米量纲不一致坡度工具会要求填Z因子。常见做法是用111319这个近似值把度换算成米但这是全球平均在赣州这种中纬度地区仍有误差。最稳妥的做法是先把dem投影到高斯或UTM坐标系再做坡度Z因子恒为1。坡向工具输出的0到360度方向角北向是0东向是90。做日照分析时还要结合当地纬度与太阳高度角不可只凭坡向断言阴坡阳坡。30m DEM用于大范围日照分区够用用于单栋建筑或小区级设计是不够的。如果手里有机载LiDAR的DSM常见做法是先做去地物处理把植被、建筑物滤掉再生成DEM这种DEM的局部精度远优于30m全球产品适合精细场景。4.3 填洼与流向做水文分析前必做的预处理水文分析的第一步不是流向而是填洼。ArcGIS的“填洼”工具输入dem输出填洼后的dem。参数只有一个Z限值默认填平所有洼地。河谷上游的小型洼地可能是真实地形全部填平后流域边界会失真。我的习惯是先给Z限值比如5米只填浅坑保留明显洼地。填洼后紧接着计算流向Flow Direction再做流量Flow Accumulation然后才能提取河网和流域。流向算法默认是D8即水流只流入8邻域中落差最大的下坡方向。D8在30m网格下对宽浅河谷的模拟比较粗对山区沟谷反而表现不错。提取河网时通常在流量累积栅格上设一个阈值比如1000个像元以上作为河网阈值越小河网越密没有绝对标准要对着影像反复试。填洼是有损操作它改变了原始高程。开工前务必保留一份原始dem后续做淹没分析或与实测高程对比时都用原始数据不要用填洼后的结果。5. DEM与shp配合避坑指南偏移、NoData、锯齿、分辨率与单位以下是拆过多份DEM数据包、翻过几次车之后沉淀下来的问题清单。每一条都按现象、原因、解决展开可以对照处理。5.1 坑一shp与dem坐标系基准不一致叠加时整体偏移现象shapefile边界线与dem上的地形错开一两百米看起来像“套歪了”。原因dem是WGS84经纬度shp却是CGCS2000高斯投影或者shp用了旧坐标系成果没有做基准转换就直接叠加。解决先用第2章的脚本分别打印两者的SRS确认椭球和投影。统一的做法是把dem和shp都转到同一个坐标系。以赣州为例常见选用CGCS2000高斯投影或WGS84 UTM投影。ArcGIS里分别用“投影栅格”和“投影”工具处理投影栅格的重采样方法选最邻近。不要依赖软件弹窗警告后点忽略继续弹窗点掉之后的结果基本是错的。5.2 坑二NoData像元被当成0坡度图出现假悬崖现象坡度分析后河谷底部和边界边缘出现一圈夸张的高坡度带地形像刀削过一样。原因裁剪或镶嵌时NoData值被写成了0或者没有识别原始NoData坡度计算时把空白区当作海拔0米参与差分产生十几米甚至几十米的伪高差。解决在任何坡度、等高线、填洼工具执行前先确认dem的NoData。如果已经混成0用栅格计算器把这些0值先设为NoData再重算。执行gdalwarp裁剪时必须用-dstnodata明确指定一个负值比如-32768让它不会默认变成0。埋在地形里的零点可以保留但边界与洞内的0要清理。5.3 坑三裁剪后的边界锯齿严重甚至出现透明缝隙现象用shp裁剪后dem边界一圈锯齿局部有透明空像元叠加到底图时露出背景色。原因多边形边界与栅格网格天然不同步。按掩膜提取时边界像元只能整体保留或丢弃像元中心落在面外的就被丢掉如果shp边界本身折点不干净缝隙更明显。解决先给shp做一个轻微缓冲比如一个像元即30米再裁剪让mask多覆盖一圈余量之后用gdal_fillnodata把空洞补上。若shp文件打不开或几何损坏先用shapechk这类工具检查和修复shp修好再拿来做掩膜否则缝隙问题会反复出现。不要为了边界顺滑而把dem重采样到更粗分辨率那是用精度换美观。5.4 坑四30m分辨率在山谷与城区不够用等高线在陡坡上交叉现象把等高线放大到1:5000看山谷细沟完全没表达陡坡处等高线挤成一团甚至交叉。原因30m像元对地形做了空间平均宽度小于30米的沟壑已经被平滑掉陡坡上3×3窗口只有9个高程点插值出来的等高线自然不稳定。解决城区、河谷、矿权范围等精细区域换用12.5m分辨率的ALOS数据或机载LiDAR生成的DEM不要硬啃30m旧数据。大范围趋势分析用30m局部精细设计用更高分辨率二者配合使用而不是只挑一个。做洪水淹没这类对地形细节敏感的模拟时尤其要意识到30m网格对窄河谷过流能力偏高这个事实。5.5 坑五地理坐标系的dem直接量距算面积数值离谱现象dem还是经纬度坐标时量一条河长得到“0.05”不知道是什么单位计算面积又得到一堆奇怪小数。原因栅格没有投影时X方向单位是度高程单位是米从像元数乘像元大小得到的只是“度平方”不能当作实际面积。量距时也一样测出来是度。解决先把dem和shp统一投影到高斯或UTM再做任何长度面积统计。投影后的像元尺寸不再是0.000几度而是确定的30米等高线间距、Z因子、面积量算都会变正常。顺手把投影后的dem另存一份作为所有派生分析的统一输入避免每个工具各投影一次造成新的错位。6. 让这套DEM数据发挥更大价值批量裁切、精度与自洽验证6.1 按县级shp批量裁切赣州下辖章贡区、赣县区、瑞金市等县级区域。用市级shp裁完整体后可以用一段bash循环按县级shp挨个裁几分钟得到一套分县dem文件。for f in 章贡区.shp 赣县区.shp 瑞金市.shp; do name$(basename $f .shp) gdalwarp -cutline $f -crop_to_cutline -tr 30 30 -r nearest \ -dstnodata -32768 赣州裁剪.tif ${name}_dem.tif done这段循环的可复制性很高cutline切换不同县级边界面输出文件名自动带上县名。如果某个县级shp坐标与整体不一致先统一投影再进循环。6.2 用实测高程点验证DEM精度精度验证别靠肉眼。常见做法是用现场实测的高程点在ArcGIS里用“多值提取至点”把dem上的高程提取出来与实测值比较计算均方根误差。30m DEM与实测点的RMSE在10米以内市级尺度分析基本可用超出20米就要警惕数据源是旧版或存在局部空洞。RMSE可以作为是否换高分辨率数据的判据而不是一句“感觉差不多”。6.3 用山体阴影和等高线做自洽验证最后一个技巧生成山体阴影并叠加等高线。山体阴影里有明显条带或方形接缝说明原始数据是分幅拼接的存在高程台阶。等高线与山体阴影的明暗分界完全背离说明dem有插值伪影这时可以用Topo to Raster按等高线重建一版做交叉验证。外业核对时把重点区域转换成kml导入手机地图现场点几个高程点即可。我现在拿到任何DEM加shp包都会先走完“坐标体检、NoData清点、投影统一”这三步再做分析。这个习惯帮我把返工率压得很低。希望帮到你。本文还有配套的精品资源点击获取
返回列表