简介:山西省晋中市三十米分辨率数字高程模型数据包,面向地理信息系统学习者、规划人员与科研工作者。数据覆盖晋中市全域,并附有市级边界矢量文件,可导入常用地理信息软件,用于地形渲染、坡度坡向分析、流域水文模拟、道路选线与环境影响评估等任务,解决区域精细地形数据获取困难的问题。压缩包共十二个文件,以高程栅格为核心,搭配边界矢量、投影参数、空间索引与元数据等辅助文件,整体约七十七点四九兆字节,结构清晰,加载方便。已有三百二十八人学习下载。数据源自遥感影像与测绘成果,海拔精度达到三十米网格标准,用户可结合边界文件快速裁剪出晋中市精确范围,为学术研究、课程设计或工程初步评估提供可靠的地形底图。
1. 晋中市DEM数字高程30m:数据包能干什么,不能干什么
前一阵有个做光伏踏勘的朋友让我帮忙算一块山坡的坡度,我第一反应是让他下载晋中市DEM数字高程30m这份数据包。现在30米分辨率的高程数据不算罕见,但多数是分幅下载,拿回来要自己拼接、裁剪、统一坐标系,一套流程下来半天没了。这个包把晋中市本市级范围的shp文件和DEM栅格放在同一个zip里,解压后直接就能拖进ArcMap做坡度、坡向、汇水区、等高线这些基础地形分析。适合谁用?搞测绘、国土空间规划、农林调查、风电光伏前期踏勘的人,尤其适合那种“只想分析晋中市这一个市域范围、不想自己拼图”的情况。
但也要把丑话说在前面:30m分辨率在市级尺度看宏观地形趋势是够用的,像山地和盆地界限、河流沟谷走向、坡度分级都能应付;可你要是做单栋建筑排水、小区级内涝模拟,或者需要识别两米宽的田坎,那还是得去找1m或5m的LiDAR数据。这份数据定位是区域分析底图,不是精细设计底图。明白这个定位,后面用起来才不容易翻车。
2. 解压与坐标系检查:用shp裁剪前必须做对的两件事
2.1 文件清单核对:别让“李鬼”shp混进来
拿到zip先别急着双击。常见做法是先解压到一个干净目录,再看文件列表。我用Linux终端习惯这样做:
unzip 晋中市DEM数字高程30m(含本市级范围shp文件).zip -d ./jz_dem ls -lh ./jz_demunzip -d指定解压到当前目录下的jz_dem文件夹,避免zip里自带多层目录造成路径混乱;ls -lh用易读格式列出文件大小。这一步重点关注两个东西:第一,DEM栅格通常是tif格式,文件名里一般带dem或dsm字样;第二,shp矢量文件不是单文件,必须有.shp、.shx、.dbf、.prj这几个基础文件同时存在。.shx存几何索引,.dbf存属性,.prj存投影信息。缺了.prj,后面ArcMap会把它当成未知坐标系,裁剪时直接导致范围错位。
如果列表里发现shp是零字节或者.prj缺失,别继续往下做,先重新获取数据。很多人在这时候选择“自己手动定义投影”,一旦猜错,整个分析结果都会偏,而且是那种很难发现的偏移。我一般会先用gdalinfo把坐标系统认一遍,再决定要不要修。文件清单这步花不了两分钟,但能省后面两个小时的排查时间。
2.2 用gdalinfo快速判断坐标系和像元大小
数据解压正常后,先确认坐标系,再谈裁剪。因为晋中市跨UTM 49N分带,DEM可能来自SRTM,原始是WGS84经纬度;而shp可能是CGCS2000高斯投影。如果一上来就裁剪,两个数据叠加后可能错位几百米,讨论裁剪就没有意义。在命令行里看基础信息:
gdalinfo DEM_30m.tif | grep -E "Size|Coordinate System|Pixel"Windows下用findstr替换grep。输出里如果看到GEOGCS["WGS_1984",说明DEM是地理坐标系,单位是度,像元大小大概是0.00027度左右;如果看到PROJCS["CGCS2000_3_Degree_GK_CM_111E",说明它已经是高斯投影,像元大小就是30。这个信息直接决定了待会怎么设置snapRaster和cellSize。另外还可以顺便看一眼Metadata,有些TIFF头会写清楚是SRTM还是ALOS,这能帮你判断数据源的高程异常值特征。
为什么要强调坐标系?因为30m分辨率在WGS84地理坐标系下,一个像元的实际地面宽度在不同纬度不一样,晋中市纬度约37°N,经度方向一个0.00027度的像元实际不到24米,而纬度方向约30米,所以做坡度或面积分析之前,最好先投影成高斯-克吕格投影或UTM投影。当然,这份数据如果已经自带投影,那这步可以跳过;如果没有,就要用“投影栅格”工具重采样。注意重采样方法不要选“最邻近”,地表高度是连续面,用“双线性”或“三次卷积”更平滑。常见做法是在投影的同时把像元大小锁定为30。
2.3 用arcpy做一致性检查与投影统一
确定坐标系之后,做一个自动化的空间参考对比。我习惯写一段小脚本,把DEM和shp都读进来,不一致就自动投影shp,这样后面在ArcMap里操作就不会因为坐标系问题翻车。
import arcpy dem = r"C:\data\jz_dem\DEM_30m.tif" shp = r"C:\data\jz_dem\晋中市本市级范围.shp" dem_sr = arcpy.Describe(dem).spatialReference shp_sr = arcpy.Describe(shp).spatialReference print("DEM:", dem_sr.name) print("SHP:", shp_sr.name) if dem_sr.name != shp_sr.name: out_shp = r"C:\data\jz_dem\晋中市_投影.shp" arcpy.Project_management(shp, out_shp, dem_sr) print("已投影shp到DEM坐标系:", out_shp)这段代码先通过Describe拿到两个数据的空间参考名称,然后比较。如果不一致,就调用Project_management把shp投影到DEM的坐标系。我一般选择投影shp而不是投影DEM,因为重投影栅格会重采样,会改变像元值;矢量只是重新计算坐标,属性不变。如果你的DEM本身不是目标坐标系,而你后续又要等高线,那也可以先把DEM投影到当地高斯坐标系,这时要把arcpy.env.snapRaster设成原始DEM,防止投影后网格发生位移。
这里要提醒一句:arcpy脚本跑不跑得起来,取决于你有没有ArcGIS Pro或Desktop的Python环境。如果没有,直接用ArcMap工具箱里的“投影”窗口操作一样,原理完全一致。判断标准很简单:裁剪之前打开地图属性,看两个图层的坐标是否在同一坐标系下,如果地图单位是度,那就说明至少有一个不是投影坐标。
2.4 多幅DEM拼接:先Mosaic再裁剪,避免边界像元打架
市级范围的DEM数据源通常不是一幅完整tif,而是按标准分幅或按图幅切好的多块。如果解压后看到好几个tif,先用gdalinfo看一眼它们覆盖范围是否有重叠,再决定拼接方案。我常用的命令是这样的:
gdal_merge.py -o jz_merge.tif -n -9999 -a_nodata -9999 tif1.tif tif2.tif tif3.tif-n -9999表示输入中等于-9999的像元视为NoData并忽略,-a_nodata -9999给输出tif统一写入NoData为-9999。这么做能避免两幅图之间出现虚假的“接缝低值”或“接缝高值”。在ArcGIS里对应的工具是“镶嵌至新栅格(Mosaic to New Raster)”,像素类型建议选“16位有符号整型”或“32位浮点”,NoData设为-9999。
先拼接再裁剪和先裁剪再拼接,我的结论很明确:先拼接后裁剪。如果你拿着shp把每一幅都裁剪一遍再拼,边界处重叠像元可能因为原始数据源版本不同而出现值跳跃,拼出来的等高线会在行政边界内侧一圈出现“台阶”;先拼接成一整幅,再按shp裁剪一次,你只需要处理一次边界,NoData也统一。这个顺序能省掉很多后续查错的时间。
3. 掩膜提取与栅格裁剪:ArcMap里哪个更适合晋中市DEM
3.1 区别:按掩膜提取和按掩膜裁剪的底层行为差异
ArcMap里裁剪DEM最常被问的问题就是“Extract by Mask”和“Clip”到底有什么区别。简单说,**按掩膜提取(Extract by Mask)**是Spatial Analyst工具箱里的工具,它先把面shp栅格化成掩膜,然后把掩膜范围内对应的原始像元值提取出来,范围外一律设成NoData。**栅格裁剪(Clip)**是数据管理工具箱里的工具,直接按面要素的几何边界去“切”栅格外接矩形,切下来的像元保持原始值;裁剪范围外的部分默认不会自动变成NoData,需要你指定NoData值。
这个差异在实际使用中非常关键。举个例子,晋中市的shp边界是弯弯曲曲的行政区界线,如果你用Clip的默认矩形模式,裁出来是一个矩形,shp边界外还带着一圈高程数据;如果你勾选了“使用输入要素裁剪几何”,边界外的数据会被切掉,但切掉的边缘可能保留半个像元,而且背景值往往是0而不是NoData。而Extract by Mask是面内保留、面外直接变成NoData,做坡度、汇流分析时NoData不会被当成0值参与运算,数据更干净。所以我的判断标准是:做分析用Extract by Mask,做快速出图可以用Clip,但Clip之后必须检查背景值。
另外一个容易被忽略的区别是性能。Clip因为不需要把shp转栅格,运行速度快很多;Extract by Mask要先跑一次栅格化,对于像晋中市这种几千平方公里的市级范围,差别也就是几十秒,不值得为速度妥协导致后续统计错误。
3.2 操作步骤:用按掩膜提取输出晋中市DEM
我一般会在ArcGIS Pro里跑这样一段Python脚本,用面要素做掩膜提取:
import arcpy from arcpy.sa import ExtractByMask arcpy.env.workspace = r"C:\data\jz_dem" arcpy.env.snapRaster = "DEM_30m.tif" arcpy.env.extent = "晋中市本市级范围.shp" dem = "DEM_30m.tif" mask = "晋中市本市级范围.shp" arcpy.CheckOutExtension("Spatial") out_raster = ExtractByMask(dem, mask) out_raster.save(r"C:\data\jz_dem\晋中市_DEM_cut.tif") print("提取完成")这段脚本最关键的是两个环境设置。arcpy.env.snapRaster指定捕捉栅格为原始DEM,这是保证输出像元网格与原始DEM完全对齐的关键;如果你不设置,ArcMap可能按默认网格重新对齐,结果面和原始数据之间错半个像元,也就是大约15米的位移,后面叠加shp时肉眼不一定看得出,但做剖面线、高程点采样时误差会暴露。arcpy.env.extent指定处理范围为shp范围,比直接让工具自己计算省内存,也防止输出范围富余过多。
如果不想写脚本,界面操作路径是:工具箱 → Spatial Analyst 工具 → 提取分析 → 按掩膜提取。输入栅格选DEM,输入掩膜数据选晋中市shp,输出栅格路径选到jz_dem目录下,然后环境设置里把“捕捉栅格”设为DEM,“像元大小”设为与DEM相同。工具箱界面里这些环境变量藏在“环境”按钮里,不是工具参数页,很多人找不到。
参数可以按这个表来设:
| 环境项 | 推荐值 | 原因 |
|---|---|---|
| 处理范围 | 晋中市本市级范围.shp | 只计算市域内,减少冗余 |
| 捕捉栅格 | DEM_30m.tif | 保证像元对齐 |
| 像元大小 | 与捕捉栅格相同 | 防止被默认值改大/改小 |
| NoData值 | 保留原值 | 便于后续坡度计算 |
3.3 如果一定要用Clip:如何按要素几何裁剪而不是外接矩形
有些项目对数据格式有要求,非要用Clip不可。比如你需要把裁剪结果直接转成WMS发布,Clip的输出更直接。工具语法是:
arcpy.Clip_management( "DEM_30m.tif", "晋中市本市级范围.shp", r"C:\data\jz_dem\晋中市_DEM_clip.tif", "#", "ClippingGeometry")第4个参数传#表示保留原始NoData,第5个参数ClippingGeometry表示使用输入要素的几何作为裁剪边界。如果这里不写ClippingGeometry,ArcMap就会把shp的外接矩形当成裁剪范围,输出还是矩形,shp没起作用。很多用户勾选“使用输入要素裁剪栅格”其实在这个工具箱里就是设置这个参数。需要特别指出,Clip的输出栅格在范围外会填0,而不是NoData,所以跑完后我通常会紧接着执行一次栅格计算器,把0值设为NoData,否则统计结果里会出现大量假0。
这里还有个小细节:Clip的输入要素可以是面,也可以是线。如果shp是线,ClippingGeometry会用线的外接矩形,但线的宽度是零,裁出来的栅格只有一行像元,这在某些工程上反而是需要的。但如果你想要的是晋中市市域内的完整DEM,一定要用面shp,而不是线shp。
3.4 QGIS用户:用Clip raster by mask layer代替
如果你用的是QGIS而不是ArcMap,不需要担心Spatial Analyst许可问题。在菜单栏选择“栅格 → 提取 → 按掩膜图层裁剪栅格”,输入栅格选DEM,掩膜图层选shp,勾选“裁剪到掩膜范围”,后台命令会自动把图层转为栅格掩膜。有一个额外技巧:在“其他参数”里加上-dstnodata -9999,可以让输出范围外的像元写成-9999,便于统一管理。QGIS的这个工具本质与Extract by Mask一致,而且基于GDAL,裁剪速度通常比ArcMap的Spatial Analyst工具更快。但要注意,QGIS的掩膜图层如果是多边形,线划边界的锯齿化算法和ArcMap不同,输出结果比ArcMap更贴近要素边界,但也导致部分边界像元被裁掉,适合做制图不适合做定量分析。我自己的习惯是:做分析还是用Extract by Mask,QGIS的裁剪只用来画预览图。
3.5 裁剪结果的验证:范围、像元数、NoData占比一眼看穿
裁剪完成不代表万事大吉。我会立刻用gdalinfo检查输出:
gdalinfo -stats 晋中市_DEM_cut.tif重点看三行:Size is后面的宽高像元数,STATISTICS_MINIMUM和STATISTICS_MAXIMUM,以及末尾的NoData Value。如果NoData Value显示为undefined或没有,背景0很可能被当成有效高程了;如果最小值和最大值都在晋中市合理海拔范围内,说明底色没有被误统计。还有一种验证法:直接在ArcMap里把DEM和shp叠在一起,设置DEM符号系统为“唯一值”或“拉伸”,看边界外是否有黑色区域。我在项目里遇到过裁剪结果范围正确但栅格几何偏了半个像元的案例,肉眼根本看不出来,最后就是靠gdalinfo里Origin坐标和shp的extent比对发现的。所以这个验证步骤我从不跳过。
4. 避坑/常见问题:裁剪DEM时我踩过的五个坑
4.1 zip解压报错“invalid zip archive: could not find EOCD”
现象:下载的zip在解压时提示invalid zip archive: could not find EOCD,或者用Python的zipfile解压到一半中断。问题出在下载文件没有下完整,EOCD是zip文件末尾的中央目录记录,缺失说明文件不完整。解决:先别急着重新下载,用7-Zip测试一下:
7z t 晋中市DEM数字高程30m(含本市级范围shp文件).zip如果输出里有ERROR,说明文件损坏,重新下载,并且下载时选择支持断点续传的客户端。如果测试通过,再解压。这个坑很典型:很多人看到zip能解压出部分tif就认为文件没问题,结果DEM在ArcMap里显示到一半变成白块。血泪经验:zip下载完成后先校验大小和哈希值,再动手。
4.2 裁剪出来还是整个矩形,shp像没生效
现象:运行Clip工具输出tif,打开一看还是覆盖整个矩形范围,shp边界完全没起作用。原因:没有在裁剪参数里选择ClippingGeometry,工具默认用的是面要素的外接矩形。很多人在界面里只看到输入栅格、输出范围,没注意到“使用输入要素裁剪栅格”那个复选框。解决:在Clip_management里显式传递ClippingGeometry,或者改用Extract by Mask。也可以用另一种判断方法:裁剪前在ArcMap里选中晋中市shp,再运行工具,如果工具窗口里“输出范围”字段自动填了shp范围,说明已经认为你选了外接矩形,需要手动把裁剪几何改成ClippingGeometry。
这个坑还有另一个变体:shp本身有多个面要素,你只选中了其中一个县域,裁剪结果就只剩一个县城。所以在运行前打开shp属性表,确认选中的要素数量。如果你希望保留全部市域范围,清除所有选中项,或者用Select by Attributes把“NAME = '晋中市'”筛选出来再跑。我踩过一次,那个结果让甲方误以为只有榆次区有DEM,害得我解释了半天。
4.3 裁剪结果四周黑边,最小统计值为0
现象:裁剪后的DEM在ArcMap里显示一圈黑色背景,用“识别”工具查看,像元值等于0,但原始DEM里并没有这些0。原因:Clip工具默认把范围外填充为0,而不是NoData;或者说原始DEM在边界外本来有0值的无效背景,被一起剪了进来,而ArcMap的默认符号系统把0渲染成黑色。解决:用栅格计算器把0转成NoData,或者直接在符号系统里设置:
Con("晋中市_DEM_cut.tif" == 0, NoData, "晋中市_DEM_cut.tif")如果你用的是arcpy,可以写成:
import arcpy from arcpy.sa import SetNull arcpy.CheckOutExtension("Spatial") out = SetNull("晋中市_DEM_cut.tif", "晋中市_DEM_cut.tif", "VALUE = 0") out.save("晋中市_DEM_no0.tif")注意:如果原始DEM在低海拔地区本身就存在0值(比如水面),这一步会把真实高程也变成NoData,建议先用属性统计确认最小值是多少。晋中市最低海拔肯定不是0,所以通常可以放心。但如果数据源是沿海或湖泊区域,就需要用VALUE < <某个阈值>来过滤,而不是直接等于0。
4.4 shp与DEM坐标系不一致,裁剪后整体偏移几百米
现象:裁剪结果边界与shp在图上显示错开,有时边界在山区,裁剪出来的范围却落在盆地边缘,偏移距离约几百米。原因:DEM是WGS84地理坐标系,shp是CGCS2000高斯坐标系,两者参考框架不同,不做统一直接分析就整体平移。也有可能是shp被误定义成了西安80或北京54,而不是CGCS2000。解决:裁剪前运行第2.3节的脚本,对比坐标系。如果shp没有.prj文件,不能用“定义投影”盲猜,要先看属性表里的经纬度值范围:经度在111到114之间,纬度在36到38之间,基本可以判断是WGS84或CGCS2000经纬度;如果坐标变成了三维的(如38530765, 4123456),那就是高斯投影。确认后再投影统一。
这里分享一个定位技巧:在ArcMap里打开“视图 → 数据框属性”,把“显示”设为“已知坐标参考”,如果地图单位是度,说明至少有一个图层没有投影。如果两个图层都显示度,但叠上去依然对不齐,那就是参考椭球体不同,需要用一个公共转换参数做七参数转换。这种情况在山西的CGCS2000历史数据里偶尔出现,建议直接使用“创建自定义地理变换”工具,参数一般从当地测绘局拿。
4.5 像元分辨率被悄悄改掉:裁剪后文件大了30倍
现象:裁剪前后DEM像元大小完全变了,原来0.00027度变成了0.00002度,或者高程值变得很奇怪,文件体积暴涨。原因:在环境设置里,“像元大小”没有固定为与输入相同,ArcMap会根据范围自动计算一个默认真像元大小,导致重采样。解决:运行前在环境设置里把“像元大小”设为“与栅格图层相同”,并在arcpy.env.snapRaster里指定原始DEM。这个坑最隐蔽,因为结果肉眼看上去没问题,但做坡度时角度变了,做面积时距离也变了,属于玄学级错误。我现在每次跑提取之前,都会先记录原始像元大小,跑完再用gdalinfo检查一遍输出,不一致就直接删掉重跑。
另外补充一种情况:如果原始DEM是浮点型,裁剪后变成了整型,高程精度直接丢掉了。原因是Clip工具的输出像素类型可能默认与输入相同,但如果你在“环境设置”里改了“像元类型”或“位深”,就会造成截断。这个坑我见过同事让DEM从30.5变成31,坡度分析全废。
5. 进阶:从30m DEM生成等高线shp并导出DXF
5.1 用Contour工具生成等高线shp
裁剪好的DEM不再需要插件,直接用Spatial Analyst的“等值线”工具就能提取等高线。对于30m数据,我一般设置等高距50m,如果是山地地形可能用25m辅助线,太细会在ArcMap里变成一团毛线,无法制图。arcpy调用如下:
import arcpy from arcpy.sa import Contour arcpy.CheckOutExtension("Spatial") Contour( r"C:\data\jz_dem\晋中市_DEM_cut.tif", r"C:\data\jz_dem\等高线.shp", 50, 0, 1)参数含义依次是输入栅格、输出线要素类、等高距(50m)、起始等高线值(0表示从0开始)、Z因子(默认1)。如果DEM已经投影成高斯投影,Z因子保持1,水平距离和垂直距离都是米,生成的等高线是准确的地面高度。如果还在经纬度下跑,Z因子需要调成一个大约111320的值,否则高宽比完全错误。所以务必先做好第2章的投影统一。
5.2 在Global Mapper里交叉验证边界和异常值
等高线生成后,我会用Global Mapper打开裁剪后的DEM和shp,开启“控制中心”查看两个图层的投影信息;如果叠加显示边界和影像完全贴合,说明坐标系没问题。再用“路径剖面”工具沿河谷方向画一条线,看剖面线的每个拐点高程是否合理。比如晋中市东部太行山区的一个沟谷,剖面高程若出现瞬间从1200跳到0再跳回1200,十有八九是0值NoData没处理干净,需要回到4.3重跑。Global Mapper还有个优势是能直接输出Blend格式的三维模型,不过不是今天重点。
5.3 导出DXF给设计院
如果项目需要把等高线交给规划或设计院,最通用的方式是导出DXF。在Global Mapper里选择等高线图层,File → Export Vector → DXF,坐标系选择“CGCS2000 / 3-degree Gauss-Kruger CM 111E”。用ArcMap的话,可以直接右键图层 → 数据 → 导出至CAD,生成DWG/DXF。这里有一个坑:导出DXF时如果坐标系选错,CAD里打开位置偏移,设计院会拿着纸质地形图来找你核对。我现在导完都会在CAD里随手画一个点,点的坐标和shp的某个路交点对照一下,确认没问题再发出去。
第一次拿到这份晋中市DEM数据包时,我也跳过坐标系检查直接拖进ArcMap,结果生成的等高线整体北偏了一百多米,现场复测时才发现问题,重新跑了一遍投影才找回来。从那以后,我每次拿到DEM都会先跑一遍gdalinfo记录坐标系和像元大小,裁剪后看一眼背景值和范围,再进等高线流程。这套流程虽然多花五分钟,但再也没有因为NoData或投影问题返过工。希望帮到你。
本文还有配套的精品资源,点击获取