简介这份资源是云南省临沧市及周边地区的30米分辨率DEM数字高程模型地理信息数据包适合GIS学习者、地理专业学生与城市规划相关人员使用。数据以GeoTIFF格式存储地形高程信息并附带临沧市行政范围Shapefile边界文件方便在ArcGIS、QGIS等软件中直接加载开展坡度分析、洪水模拟、可视域分析等练习。压缩包共12个文件包含tif主数据、shp/dbf/prj/sbn/sbx/shx等矢量要素及索引文件、tfw空间参考文本与多个xml元数据记录整体大小约87.69MB。已有330人学习下载包内文件结构完整既能支撑地形分析实践也有助于理解GIS栅格与矢量数据的组织逻辑。1. 拿下临沧市30m DEM压缩包先弄懂这份数据能干什么以前做滇西南项目最难的不是分析是找一份能直接用的高程数据。临沧那种山高谷深、动不动就是几十公里无人区的地方90m分辨率的DEM压根看不出支沟和陡坎自己用GPS实测又跑不过来。所以当看到云南省临沧市DEM数字高程数据30m含区域范围shp文件.zip这种压缩包第一反应是找数据这步终于省了。它本质上就是一张覆盖临沧全境的30m分辨率高程栅格TIF外加一个临沧市或指定研究区的shp矢量边界把DEM和行政范围捆在一起发省去了自己从全球分幅数据里裁剪的工序。适合谁常年在arcmap里做地形分析、坡度坡向提取、河网划分的GIS工程师以及水利、交通、地灾、林业这些需要快速出地形底图的从业者。2. 解压后先查家底投影、Value范围、边界三个必检项拿到DEM之后别急着拖进arcmap里就开始算坡度先花五分钟把三个基础项查清楚坐标系、高程值域、与shp的范围匹配。数据出问题多半不是玄学全是这几个参数没对齐引起的。下面按我处理DEM的固定顺序说一遍。2.1 打开tif先确认坐标系WGS84还是CGCS2000直接决定要不要投影在ArcCatalog里右键tif选属性或者直接在ArcGIS Pro的目录窗格里看栅格源信息重点看空间参考里的两个字段坐标系名称和投影坐标系名称。如果显示 GCS_WGS_1984 或 GCS_China_Geodetic_Coordinate_System_2000说明是地理坐标系单位是度后面做坡度坡向、算面积长度之前必须先投影如果显示 WGS_1984_UTM_Zone_47N 或 CGCS2000_3_Degree_GK_CM_99E 这种带投影的名称说明单位已经是米可以直接做距离和面积运算。沉余数据我在ArcCatalog里看批量处理时就用arcpy一行解决import arcpy raster_path rD:\lincang\data\dem30m.tif desc arcpy.Describe(raster_path) sr desc.spatialReference print(坐标系名称:, sr.name) print(投影坐标系:, sr.projectionCode if sr.projectionCode else 无地理坐标系) print(线性单位:, sr.linearUnitName if sr.linearUnitName else 度) print(WKID:, sr.factoryCode)这段代码的逻辑是获取栅格的空间参考对象再分别输出名称、投影代码和线性单位。投影代码有值时说明数据已经带了投影坐标系没有时说明仅地理坐标系。注意arcpy.Describe的linearUnitName属性只有投影坐标系下有值地理坐标系会打印度这正好帮我们判断后面要不要先投影。shp也用同样方法查一遍两个数据的坐标系完全一致最好不一致就用Project或ProjectRaster统一这个操作放到第3章裁剪前完成。2.2 看Value类型与高程范围临沧的海拔带和NoData长什么样第二步是打开栅格属性表或者图层属性里的源选项卡看像元类型和统计数据。临沧市整体处于横断山系南段余脉澜沧江、怒江穿境而过河谷地带海拔最低在几百米分水岭山脊普遍在2500米到3000米以上所以一份正常的临沧DEM高程值应该落在一个从低海拔江面到高海拔山脊连续分布的区间里。右键图层属性在符号系统选项卡里查看直方图或者用栅格计算统计之后看波段统计信息。正常情况直方图是一条连续的山丘形曲线值域大概从几百米延伸到三千多米不正常的情况是在0附近或者数据两端出现一根孤立的柱子这种极端值常见的是-32768、-9999、32767之类它们不是真实高程是NoData标记或数据生产时的填充值。如果带着这些值去填洼、算坡度ArcGIS会当成真实地形处理于是你会在河谷里看到一个断崖式的假坎或者凭空多出一个负海拔的坑。发现这类值后用SetNull把它置为NoData或者用CopyRaster工具重新写一遍并指定NoData值再做后续分析。2.3 对比DEM覆盖范围与shp边界边角是差一个像元还是差一块把shp直接拖到DEM图层上给shp设置一个醒目的红色空心符号叠加看一遍边界。这一步信息量很大shp完全落在DEM范围内且离DEM边缘还有几十个像元以上是理想状态shp有的地方伸出DEM边界说明裁剪后边缘会缺数据shp整体和DEM对不上差出几公里那基本就是坐标系不一致。把三个检查项汇总成一张自查表贴在工程里随时对照检查项操作方法合格标准常见异常坐标系右键属性里的源选项卡或arcpy.Describe与shp完全一致地理坐标与投影坐标混用高程值域图层属性的直方图500m至3000m连续分布出现-32768、-9999等极值范围匹配shp叠加DEM目视检查shp完全落在tif内整体偏移或边界伸出如果三项都正常这份数据可以直接进入裁剪流程如果第一项就有问题先去投影。很多人在arcmap里看到数据显示叠上了就觉得没问题忘了动态投影只影响显示不影响实际坐标计算后面做掩膜提取时错误会原形毕露。3. 用shp面图层裁剪dem栅格tif掩膜提取与Clip的取舍这是整个压缩包最核心的使用方式依靠区域范围shp面图层把30m的DEM栅格裁剪成研究区范围。其实两条路都通按掩膜提取Extract by Mask和裁剪栅格Clip Raster。选错的人不在少数成矩形、黑边、像元错位基本都是这两个工具用岔了。3.1 解压目录规划与栅格环境设置先把zip解压到没有中文、没有空格的路径下面。我习惯把工程目录固定成三个子目录data放原始DEM和shpoutput放裁好的结果scratch放中间缓存。这样做的直接原因是arcpy和地理处理工具对中文路径的支持时好时坏早期我吃过新建文件夹2这种路径的亏报错报得莫名其妙最后把路径全改成英文一次跑通。mkdir -p D:/lincang/data mkdir -p D:/lincang/output mkdir -p D:/lincang/scratch目录建好后把压缩包里的dem30m.tif和lincang_boundary.shp文件名以实际解压结果为准都放到data目录下。解压过程中注意一件事检查shp是不是一整套完整文件因为shp文件由.shp、.shx、.dbf、.prj等多个散文件组成zip里缺任何一个arcmap打开都会报无法识别或者属性表空白。发现有掉文件的情况重新解压一遍比手动补文件靠谱。3.2 用按掩膜提取Extract by Mask裁剪DEMarcpy脚本与参数说明这是我最推荐的处理方式。Extract by Mask会把shp面内的栅格像元保留面外的像元直接设为NoData输出形状紧贴shp边界不会留一个矩形大空白。前提是DEM和shp坐标系已经统一否则裁出来的边界会明显偏移。import arcpy from arcpy.sa import * arcpy.env.workspace rD:\lincang arcpy.env.overwriteOutput True # 让输出像元对齐原始DEM避免按掩膜提取后网格偏移半个像元 arcpy.env.snapRaster rD:\lincang\data\dem30m.tif in_dem rD:\lincang\data\dem30m.tif in_shp rD:\lincang\data\lincang_boundary.shp out_ras ExtractByMask(in_dem, in_shp) out_ras.save(rD:\lincang\output\dem30m_clip.tif) print(按掩膜提取完成)逻辑上ExtractByMask把shp的几何当作掩膜凡是落在掩膜外的栅格像元全部被过滤成NoData掩膜内的像元高程值保持不变不插值不重采样。代码里最关键的一行是arcpy.env.snapRaster它把输出栅格的像元对齐到原始DEM的网格上否则当shp边界压在像元中央时掩膜提取可能因为容差问题让边缘像元重排输出网格和原DEM错开半个像元。这个设置在批处理多个shp时尤其重要能保证每次裁出来的栅格在空间位置上是严格对齐的后面做镶嵌和叠加不会出现差半个像元的玄学错误。3.3 用Clip Raster裁剪和掩膜提取到底选哪个一张表讲透很多老手习惯直接用【数据管理工具】里的裁剪Clip这个工具默认行为坑了不少人它默认按shp的包络矩形来裁剪结果得到一个边界四四方方、大面积都是NoData的矩形。如果要用shp的真实形状裁剪必须手动把剪裁几何ClippingGeometry设成ClippingGeometry否则不生效。对比项Clip RasterExtract by Mask输出边界默认shp包络矩形可选真实几何直接输出shp几何形状NoData处理需单独指定NoData值掩膜外自动置NoData像元对齐默认跟随输入栅格受snapRaster控制适用场景批量分幅切图、规则矩形输出研究区边界裁剪ArcGIS Pro环境下的Clip Raster脚本写法如下import arcpy arcpy.Clip_management( in_rasterrD:\lincang\data\dem30m.tif, rectangle#, out_rasterrD:\lincang\output\dem30m_clip_rect.tif, in_template_datasetrD:\lincang\data\lincang_boundary.shp, nodata_value-9999, clipping_geometryClippingGeometry, maintain_clipping_extentMAINTAIN_EXTENT )这段代码里要特别注意clipping_geometry参数它接收两个值Extent表示用模板要素的包络矩形ClippingGeometry表示用要素的真实几何边界裁剪。不传这个参数工具就会按矩形裁剪输出的栅格范围比shp大一圈。maintain_clipping_extent设为MAINTAIN_EXTENT后输出栅格会保持裁剪后的完整范围而不是缩到数据实际覆盖范围。两个工具都能完成用面图层裁剪DEM这件事只是控制粒度不同我个人偏好是一次性的研究区裁剪用Extract by Mask批量分幅切图用Clip且走矩形范围效率高、结果可控。3.4 NoData与输出像元切完不黑边的小细节裁剪完成后在arcmap里打开如果边缘出现一圈刺眼的黑边先别急着骂数据。这个黑边有两种来源一种是shp边界本身超出了DEM数据覆盖范围边缘像元没有高程值所以显示成NoData另一种只是渲染问题NoData默认用黑色显示看起来像一圈黑线圈。判断方法很简单用识别工具点一下黑边处的像元如果属性显示NoData那就是真缺数据如果显示一个正常的高程值那只是符号渲染问题把图层背景色改成透明或者给NoData设一个透明色就行。针对真缺数据的边角常见做法是对shp做一个30米或一个像元大小的Buffer把裁剪范围略微向外扩让边缘落在有值区域里。注意Buffer会改变shp的几何形状涉及面积统计时要在文档里记录清楚。输出像元大小我一般不做特殊设置保持和原始DEM一致通过snapRaster来保证对齐这样后面做坡度、坡向计算时每个像元的高程邻域关系跟原DEM完全一致。4. 把30m DEM变成生产成果坡度、坡向、等高线、河网一条龙裁剪只是准备工作真正能交出去的是从DEM里派生出来的成果坡度、坡向、等高线、河网shp。这四个东西在临沧这种山区地形里都有实际用途比如坡耕地调查、山洪沟划分、地质灾害易发性分区。下面按生产流程走一遍。4.1 先把DEM投影到UTM 47N算坡度前的关键一步临沧市经度约在东经98度到101度之间UTM 47带正好覆盖这片区域。如果原始DEM是WGS84经纬度坐标水平单位是度、垂直单位是米直接算坡度会出现一个经典错误ArcGIS默认把水平单位度当成米来算坡度值会整体失真。所以第2章查出来的坐标系在这里发挥作用地理坐标系必须先行投影。import arcpy arcpy.env.workspace rD:\lincang arcpy.ProjectRaster_management( in_rasterrD:\lincang\output\dem30m_clip.tif, out_rasterrD:\lincang\output\dem30m_utm.tif, out_coor_systemarcpy.SpatialReference(32647), resampling_typeBILINEAR, cell_size30 )参数这里逐个解释out_coor_system写32647这是WGS 1984 UTM Zone 47N的EPSG编号临沧大部分区域都能被这个带覆盖如果shp用的是CGCS2000坐标就把输出坐标系改成CGCS2000对应的高斯投影但临沧的DEM绝大多数还是WGS84系列的UTM投影够用。resampling_type选BILINEAR双线性因为高程是连续表面用最近邻会破坏地形平滑性用三次卷积精度高但处理慢双线性是性价比最高的选择。cell_size指定30保持原始分辨率不变。如果原始DEM本身已经是投影坐标比如打开时看到WGS_1984_UTM_Zone_47N这一步可以跳过但要看一眼投影名称和UTM带号是否和临沧匹配不匹配的话面积和坡度计算结果会有系统性偏差。4.2 坡度、坡向的最小可用命令坡度是这组数据里最常用的派生产物。临沧山区坡面陡常用两种输出方式度DEGREE给成果报告用百分比PERCENT_RISE给工程土方计算用。坡向则反映坡面朝向对日照分析、植被分布研究有用。from arcpy.sa import Slope, Aspect import arcpy arcpy.env.workspace rD:\lincang arcpy.env.snapRaster rD:\lincang\output\dem30m_utm.tif dem rD:\lincang\output\dem30m_utm.tif slope_deg Slope(dem, DEGREE, 1) slope_deg.save(rD:\lincang\output\slope_deg.tif) slope_pct Slope(dem, PERCENT_RISE, 1) slope_pct.save(rD:\lincang\output\slope_pct.tif) aspect Aspect(dem) aspect.save(rD:\lincang\output\aspect.tif)代码逻辑很好理解Slope工具接收DEM、输出测量单位和Z因子三个关键参数Z因子在这里必须是1因为dem30m_utm.tif已经是水平垂直同单位的投影坐标系。z因子这个参数的用途是校正水平单位与垂直单位不一致时的比例关系很多人拿经纬度DEM直接算坡度还在Z因子里填1结果山区大片出现接近90度的假象就是这个比例没校正。Aspect工具更省心不需要Z因子输出范围是0到360度平地赋值为-1。参数上有一处值得提坡度输出是度还是百分比完全看下游用途。做灾害分区图一般用度做土地整治和道路设计习惯用百分比。两个都算一次成本极低我一般同时输出后面画图时再选配色方案。4.3 从DEM提取等高线shp等距参数怎么给等高线是很多政府报批项目里必须有的成果。从DEM生成等高线核心参数只有一个等高距。临沧的地形高差大20米等高距在坝区和平缓坡地会很密50米等高距在峡谷区又太稀我给全景图、县城周边这类成果一般用20米给全流域概览用50米。实际操作时还要一个起始高程能把等高线和已知高程点对上。from arcpy.sa import Contour arcpy.CheckOutExtension(3D) dem rD:\lincang\output\dem30m_utm.tif out_lines rD:\lincang\output\contour_50m.shp Contour(dem, out_lines, contour_interval50, base_contour500, z_factor1) arcpy.CheckInExtension(3D)Contour工具把连续高程表面分解成指定间隔的矢量等高线输出是shp线要素。contour_interval50表示每隔50米生成一根线base_contour500表示从海拔500米开始追踪这个值最好对应调查区域河谷最低点附近能让等高线落在整百或整五十的位置图面更整洁。z_factor1在前面的投影步骤已经保证了单位一致这里不要再填别的值。需要注意这个工具依赖3D Analyst扩展代码里用了CheckOutExtension来获取许可跑完释放许可。等高线提取出来往往会有锯齿感这是30mDEM原始精度决定的不是工具问题。要求不高直接用要求高的话先在DEM上做一步3x3的FocalStatistics均值滤波再提取平滑后的等高线代价是微地形细节会被抹掉一部分。4.4 提取河网shp水文分析五个工具的串联在临沧这种河流密集、山洪沟遍布的区域河网提取是高频需求。用30m DEM提取河网准确度能满足中小流域尺度分析单个沟谷也能识别出来流程是标准的五步填洼、流向、汇流累积量、河网栅格、转矢量shp。from arcpy.sa import Fill, FlowDirection, FlowAccumulation, Con, StreamOrder, StreamToFeature dem rD:\lincang\output\dem30m_utm.tif fill Fill(dem) fd FlowDirection(fill) acc FlowAccumulation(fd) threshold 1000 stream Con(acc, 1, 0, Value {}.format(threshold)) order StreamOrder(stream, fd) streams StreamToFeature(order, fd, rD:\lincang\output\streams.shp)一步步拆着说。Fill填洼是第一步也是很多人会跳过的一步DEM里总有因插值算法产生的封闭洼地不填掉会在洼地处断流河网直接断裂。FlowDirection计算每个像元的水流方向用D8单流向算法适合丘陵山区平坝地区多条流向算法效果更好但不展开。FlowAccumulation统计每个像元的上游汇流像元数结果是累计栅格。threshold这个参数是整个流程里最值得调试的它表示当汇流累积量大于某个值时被认为存在河流。1000像元在30m分辨率下意味着上游集水面积约0.9平方公里适合临沧这种中小流域密集的地形。想提取更多毛沟就降到500想只看主干河流就调到2000以上参数需要根据当地经验和流域面积试两三次。Con生成河网栅格Value大于阈值的像元赋1其余赋0。StreamOrder给河道分级StreamToFeature把栅格河网转成shp线文件落到磁盘上后面叠加到自己的地图里出图用。那步Con的赋值逻辑是true1、false0因为StreamOrder和StreamToFeature只识别有河流的像元0代表背景。如果你在ArcMap里操作用栅格计算器写一行SetNull(acc 1000, 1)效果也一样。5. 避坑临沧DEM数据处理最容易翻车的5个问题下面这篇数据踩过的坑早期基本一一翻过车按现象、原因、解决写出来省得再摸一遍。5.1 坐标系不一致shp和tif叠加后整体偏移几百米现象把zip里的边界shp拖到DEM上两者的边界错开一段距离有时几百米有时甚至跑到相邻乡镇去了。原因DEM是WGS84经纬度坐标系shp是CGCS2000高斯投影坐标系或者反过来。arcmap默认开启动态投影每个图层都在显示时重新投影所以肉眼看是叠上的但地理处理计算用的是原始坐标系掩膜提取、裁剪这类工具一跑结果就错位。解决处理前统一坐标系。以shp为准用ProjectRaster把DEM投到shp同一坐标系下或者反过来把shp投成DEM的坐标系。选哪一边以最终出图需要的坐标系为准没有特殊要求就统一到CGCS2000高斯投影成果拿去对接规划部门时不容易被质疑坐标系合法性。5.2 高程值出现负值或32767别当真实地形算现象直方图里出现一根孤零零的柱子DEM上对应区域显示成黑洞或者一块突兀的高台识别后高程值是-32768或32767。原因原始数据生产时为非数据区填充了极值这类值本质是NoData标记。SRTM、ASTER GDEM数据都这么干过压缩包里的数据如果是从全球分幅数据裁剪拼接的边缘或云覆盖区就会留下这种值。解决先通过栅格统计找出极值位置再用SetNull把这些值转成真正的NoData。处理时注意区分负海拔区域临沧全域没有负海拔所以把所有小于0的值和大于3000米以上的极端值先筛出来逐一确认位置再处理不要在没确认前就一刀切。5.3 提取坡度后大片接近90度Z因子和投影的坑现象坡度图在山脊和沟谷处出现成片的红色高值区很多地方直接顶着89.9度与实际地形完全不符。原因90%是坐标系问题。DEM还是经纬度坐标系就直接传给了Slope工具水平单位度被当成米处理实际地形坡度被极度放大。剩下的10%是Z因子填错或者填了0导致垂直方向被压缩或放大。解决回到第4章先把DEM投影到UTM 47N再算坡度Z因子设为1。如果必须保留原始坐标系才能算Z因子就需要填一个很大的数但算出来不直观别费这个劲直接投影。判断结果是否合理有一个简单标准临沧山地绝大多数坡面在15到60度之间如果坡度分布图上超过70度的像元占比超过百分之几流程里一定出了问题。5.4 等高线锯齿状且碎线多DEM里的噪声被放大现象提取出来的50米等高线弯弯曲曲像毛线团短的闭合碎线到处散落经不起细看。原因30m DEM本身的噪声在等高线追踪时会被线性放大陡峭地形尤其明显另外栅格转矢量时没有做简化像元级锯齿直接继承到线要素里。解决在提取等高线之前先对DEM做一次3×3的FocalStatistics均值滤波能显著削减噪声再提取等高线之后对线要素用SimplifyLine做一次简化去掉不必要的微小弯折。代价是滤波会让陡崖边缘轻微变钝对一般成果精度要求足够。注意滤波后的DEM只用于出等高线不要用它再算坡度和填方量否则地形被磨平工程量算出来会偏小。5.5 裁剪结果仍是矩形Clip默认行为没改现象明明用的是shp边界裁剪出来的DEM却是个四四方方的矩形大量像元没有数据文件大小也没变小。原因使用Clip Raster工具时没勾剪裁几何或者没把clipping_geometry参数设为ClippingGeometry工具默认按要素的包络矩形输出X方向和Y方向的最大最小值围成一个框把整个框都保留了下来。解决回到第3章的脚本ClippingGeometry传进去再跑一遍在ArcMap里操作则勾选使用输入要素裁剪几何选项。如果用的是Extract by Mask则不存在这个问题它天然按shp真实形状切。这步改完之后用属性或者范围工具验证一下输出范围的四个角应该紧贴shp轮廓而不是横平竖直的矩形。6. 进阶把这份30m DEM用得再细一点质量自检与更高精度预处理高程数据交出去之前我做两步快速自检不算复杂但能挡下大部分低级错误。第一步把dem30m_utm.tif和坡度图叠加上研究区shp目视检查山脊线和河谷形态看有没有突兀的矩形块、条带边界或者一大片平地这些特征通常指向填洼参数错误或残留NoData。第二步把原始tif拖进Global Mapper 14网上有汉化版界面顺手切换到3D视图把垂直夸张设到2倍沿着澜沧江或怒江的流向拉一条剖面线剖面出现断崖式跳变的地方基本就是条带噪声或者NoData残留区。没有Global Mapper的话用arcmap的3D Analyst里栅格剖面也能完成同样验证。还有一个概念在日常项目里经常混淆DSM和DEM。如果这份数据没有经过严格的植被和建筑剔除那它本质上是DSM临沧山区植被覆盖度高植被茂密区域的DSM高程可能比真实地面高出十几二十米做工程断面和土方估算时这个误差会直接体现在报表里。所以拿到数据先看它的产品说明无法确认来源时只把它用于宏观地形分析不用于单点高程结算。这几十个DEM项目跑下来我一个雷打不动的习惯是原始dem30m.tif永远保持只读所有裁剪、填洼、重采样结果都写到output和scratch子目录。中间产物翻车了就删掉重来原始数据始终留一副干净的底版比什么后悔药都管用。高程数据是很多分析的地基地基没动过后面才有回头路可走。希望帮到你。本文还有配套的精品资源点击获取