简介:河南省焦作市数字高程模型数据压缩包面向GIS学习者和地理空间分析人员,内含焦作市全境三十米分辨率数字高程数据,并附带市级行政区划边界矢量文件。数据可直接用于地形分析、流域划分、坡度坡向计算、可视性分析与区域规划等场景,是理解豫西北山地与平原过渡地带地貌特征的基础资料。压缩包共十二个文件,大小约十一点七九兆字节,核心为高程栅格数据文件,同时配有边界矢量文件及其属性表、投影信息、空间索引、坐标配准文件与元数据说明等辅助文件。各文件协同使用,可在常见地理信息系统软件中完整显示、叠加与编辑,无需额外转换即可开展后续分析。已有六百四十九人学习下载。数据属于中等分辨率,适合宏观地形研究、灾害风险评估、水资源管理和基础设施规划等应用。附带的市级边界便于按行政区划裁切统计或与人口经济等属性数据连接,是一份可直接上手、兼顾教学练习与实际项目的地理信息数据资源。
1. 焦作市30m DEM数据包:先看这里再动手
拿到这个"河南省焦作市DEM数字高程数据30m(含市级范围shp文件).zip",你大概率是要做地形分析、选址评估或者出图。先说结论:里面那一份30米分辨率的DEM栅格数据,配合市级边界的shp矢量文件,是你在ArcGIS或QGIS里能直接用的标准组合。30米分辨率意味着每个像元代表30米的地面真实尺寸,做市域尺度的坡度、坡向、流域分析足够,但你要是想看清一条沟的细节,就得换5米或1米的数据。这个包适合谁?做规划、测绘、环境工程的从业者,或者刚入门的GIS学生,需要焦作市范围的高程底图做分析,不想再满网找数据。常见做法是先解压、验证坐标系、再用shp裁剪DEM,最后生成坡度坡向。这篇文章就把这条链路拆开讲透。
2. 解压后先别急着开图:文件清单、坐标系统和分辨率的核实
2.1 先确认DEM和shp是不是同一套坐标系
解压zip后,你可能会看到一堆文件:dem栅格通常是tif格式,shp矢量则是.dbf、.prj、.shp、.shx一整套。很多人犯的第一个错,是直接把tif拖进ArcMap,再把shp拖进来,结果两者对不上——一个在河南,一个漂到海里。这不是数据坏了,多半是坐标系不认识对方。
先用ArcCatalog或者QGIS的图层属性面板查看投影信息。在ArcMap里右键图层,打开属性,切到"源"选项卡,看"空间参考"下拉。常见有两种情况:一种是GCGCS2000_3_Degree_GK_Zone_39,这是经纬度的投影坐标版本;另一种是WGS_1984_UTM_Zone_49N,把地球按6度带切开。焦作市位于东经112°到113°附近,投影带一般是39带或49带。如果shp和DEM坐标系不一致,先别急着裁剪,需要做投影转换。
一句经验:如果是CGCS2000的经纬度坐标,你直接在ArcMap里用"定义投影"强改,但不要乱选。正确做法是用"投影"工具(Data Management Tools → Projections and Transformations → Feature → Project)把shp或DEM统一到一个坐标系。栅格投影用"投影栅格"(Project Raster),会重采样,注意选"双线性"或"三次卷积"保留高程精度,别用最近邻。
2.2 用ArcMap或QGIS目视验证范围是否对齐
坐标系确认过后,把DEM和shp一起加载进去。一般能看到shp边界正好压在DEM的范围内。如果边界和栅格有明显偏移,哪怕几百米,说明数据源本身有问题——比如shp是焦作市最新的行政边界,但DEM是基于旧版本拼接的,两者边界有偏差。这时候不要硬来,先看shp范围是否超出DEM覆盖区。
我在实际项目里遇到过:DEM覆盖的是焦作市全境,但shp里多了一块"示范区"或"高新区"飞地。这时你需要用shp去裁剪DEM,但裁剪前确认shp是否包含多个要素。用"按属性选择"看看Feature有多少。如果shp有几十个面要素,裁剪时是否要合并边界?通常市级范围shp是一个整体面,偶尔也会有separate岛屿。
2.3 用gdalinfo快速读取元数据
如果你装了QGIS或GDAL环境,命令行里一条gdalinfo就能看到所有关键信息。打开cmd或终端,输入:
gdalinfo D:/jiaozuo_dem_30m/dem30.tif输出里重点关注这几行:Size is 1620, 1080(行列数)、Pixel Size = (30.0000000, -30.0000000)(分辨率)、Origin = (112.0000000, 36.2000000)(左上角起点)、Coordinate System is(坐标系)。Pixel Size必须是正的,负号表示行方向是向下,这是正常的。如果Pixel Size出现30.0006之类的小数,说明数据做过重采样或拼接,不是原生严格30米,后续坡度计算有一定误差,但市域尺度能接受。
对于shp,可以用ogrinfo查看。例如:
ogrinfo D:/jiaozuo_dem_30m/jiaozuo.shp -al -so这条命令会列出shp的要素数量、几何类型(Polygon)、投影信息。如果显示ERROR: Unable to open datasource,说明shp文件不完整——缺少.dbf或.shx。这时别急着骂数据,很多时候是zip解压时丢文件,重新解压到全英文路径下问题就没了。
3. 用市级范围shp裁剪DEM:面图层裁剪与掩码提取的正确姿势
3.1 常规裁剪工具和掩码提取到底什么区别
这是很多人混淆的点。ArcToolbox里有两个工具最常用:"裁剪"(Clip)和"按掩膜提取"(Extract by Mask)。前者在"数据管理工具 → 栅格 → 栅格处理 → 裁剪",后者在"空间分析工具 → 提取分析 → 按掩膜提取"。
区别在于:裁剪工具的速度快,且允许你在勾选"使用输入要素裁剪几何"时保留输入栅格范围外的区域(其实那样等于没裁)。它只做矩形裁剪,如果shp是圆形或不规则形状,裁出来的栅格仍然是矩形,只是范围被shp的包围盒限制,边界外的像元并不是空值,而是保留原始值,所以需要后续再做一步"设为空值"或设置背景色透明。而按掩膜提取则直接用shp的几何边界,把边界外的像元变成NoData,裁出来的栅格边缘就是shp的形状。这一点在出图时尤其重要——你不想整个矩形边界盖住焦作市外的区域,地图上会出现一块"白板"。
如果你只是想把数据范围缩小到焦作市周边,用普通裁剪就行;但要生成沿市界形状的DEM,必须用按掩膜提取。下面我以ArcMap 10.8为例操作。
3.2 ArcMap里的按掩膜提取参数怎么设
打开Spatial Analyst工具箱,找到"按掩膜提取"双击。输入栅格选DEM;输入掩膜数据选shp;输出栅格填一个路径,比如D:/jiaozuo_dem_30m/dem_clip_utm.tif。这时有个关键参数:"NoData值"和"提取范围"。
工具默认提取范围是"输入栅格数据与掩膜数据交集",一般不用改。但注意:如果shp有多个不相连的面,工具会按每个面元素分别裁剪,中间部分自动变成NoData。你要的是整个市域连成一块,所以最好在裁剪前用"融合"(Dissolve)工具把shp合并成一个要素。操作是打开"数据管理工具 → 要素类 → 融合",输入shp,融合字段随便选一个(比如"NAME"或"FID"),生成jiaozuo_dissolve.shp,然后用这个融合后的shp做掩膜。
如果你命令行更顺手,可以用GDAL的gdalwarp来做,一步到位。比如:
gdalwarp -cutline D:/jiaozuo_dem_30m/jiaozuo_dissolve.shp -crop_to_cutline -dstnodata -9999 -tr 30 30 D:/jiaozuo_dem_30m/dem30.tif D:/jiaozuo_dem_30m/dem_clip.tif这里-cutline指定shp,-crop_to_cutline表示裁到shp边界,-dstnodata -9999把原背景值设为NoData,-tr 30 30强制输出分辨率30米。注意,如果原图坐标系是投影坐标,输出会保持原投影;如果原图是地理坐标(经纬度),-tr 30 30代表30度而不是30米,这时要写-tr 0.000277 0.000277约等于30米。更稳妥的方式是在gdalwarp里加-t_srs EPSG:32649,先把坐标系转到UTM再设置分辨率。
3.3 用arcpy批量跑多个行政区
如果你手头有多个县的shp,比如要按山阳区、解放区分别裁,可以在ArcMap的Python窗口里用arcpy循环。代码如下:
# -*- coding: utf-8 -*- import arcpy arcpy.env.workspace = r"D:/jiaozuo_dem_30m" arcpy.env.outputCoordinateSystem = arcpy.SpatialReference("WGS 1984 UTM Zone 49N") dem = r"D:/jiaozuo_dem_30m/dem30.tif" city_shp = r"D:/jiaozuo_dem_30m/county.shp" out_dir = r"D:/jiaozuo_dem_30m/clip_result" # 先融合,避免多个小面碎? 这里直接用county.shp,要素字段为XZQMC fields = ["XZQMC"] with arcpy.da.SearchCursor(city_shp, fields) as cursor: for row in cursor: name = row[0] out_tif = out_dir + "/" + name + "_dem.tif" arcpy.gp.ExtractByMask(dem, city_shp, out_tif) print(name + " done")这段代码遍历shp的行政区名称字段,对每个行政区生成一个裁剪后的DEM。注意ExtractByMask里掩膜是shp字段的全部要素,不是当前循环的单个要素——如果你按属性选择,需要先用arcpy.MakeFeatureLayer加SelectLayerByAttribute。上面这个写法其实会把整个市域裁到每个县,只是文件名不同,是个坑。
正确做法是:
import arcpy arcpy.env.workspace = r"D:/jiaozuo_dem_30m" dem = r"D:/jiaozuo_dem_30m/dem30.tif" county_shp = r"D:/jiaozuo_dem_30m/county.shp" out_dir = r"D:/jiaozuo_dem_30m/clip_result" with arcpy.da.SearchCursor(county_shp, ["XZQMC"]) as cursor: for row in cursor: name = row[0] where_clause = "XZQMC = '{}'".format(name) temp_layer = "lyr_" + name arcpy.MakeFeatureLayer_management(county_shp, temp_layer, where_clause) out_tif = out_dir + "/" + name + "_dem.tif" # 掩膜提取,使用单个县图层 arcpy.gp.ExtractByMask(dem, temp_layer, out_tif) print(name + " done")参数说明:where_clause是按名称筛选;MakeFeatureLayer生成一个临时图层,只包含当前县范围;ExtractByMask输出结果就只有一个县的边界形状。这样批量操作时,每个县的高程范围都能准确裁出。
4. 数据落地时的四个高频坑:避坑排查手册
4.1 现象:裁剪出的DEM边缘发黑,像一圈黑框
原因:原始DEM有不少NoData区或背景值0值,裁剪后NoData被渲染成黑色。尤其在ArcMap默认拉伸符号里,NoData通常显示为黑色。
解决:在图层属性的"符号系统"里,把"NoData"颜色的勾选去掉,或者设置成白色背景。如果背景值不是NoData而是0,需要用"栅格计算器"把0值改为NoData:
# 在ArcMap栅格计算器里输入 Con("dem_clip.tif" == 0, -9999, "dem_clip.tif")把0变成-9999,再使用SetNull或直接在图层属性中设置。注意,焦作市范围内海拔最低也有几十米,不可能有0值,所以0通常是无效背景。如果计算坡度时黑边干扰了分析,一定要先做这一步。
4.2 现象:shp和DEM套在一起,shp跑到栅格外边几百米
原因:两个数据不是同一套投影。我见过最典型的情况是shp用的是CGCS2000 / 3-degree Gauss-Kruger CM 114E,而DEM用的是WGS 84 UTM 49N,两者虽然都是米制,但中央经线和椭球体不同,水平偏移能达到几百米。
解决:用"投影"工具把shp转换到DEM的坐标系。在ArcToolbox里选择"数据管理工具 → 投影和变换 → 要素 → 投影",输入shp,输出坐标系选择DEM的坐标系。转换方法选"PROJCS"里的已有方法,不要选"默认"。转换后重新加载,边界就能对上。如果你发现DEM是地理坐标(度),shp是投影坐标,思考一下:通常DEM不会以地理坐标存30米分辨率,因为30米在纬度上是约0.000277度,在经度上会根据纬度缩放,不利于分析。这种情况大概率是DEM的数据说明写错了,先检查gdalinfo里的Coordinate System is。
4.3 现象:裁剪出来的DEM范围比shp边界大,沿边界有一圈台地
原因:栅格像元是按行规则排列的,掩膜提取虽然把外部变为NoData,但边界像元只要有一部分落在shp内部,就会被保留,导致边界呈锯齿状,边缘往外扩半个像元(15米)。
解决:这是栅格数据固有的"像元对齐"问题。要减小误差,可以在裁剪前使用"重采样"工具把像元对齐到shp边界的精确网格。更简单的方法是接受这15米的误差,市域尺度分析无伤大雅。真正要命的是如果你用边界做高程剖面,剖面线正好压在边界上,会出现阶梯状跳变。这时可以给shp做一点点内缩(Buffer负值)比如-15米,再裁剪。但这个方法会把边界偏移,不推荐。更好的是在gdalwarp时使用-crop_to_cutline同时配合对齐:
gdalwarp -cutline jiaozuo.shp -crop_to_cutline -tr 30 30 -tap dem30.tif dem_clip.tif-tap参数会让输出栅格网格对齐到shp边界,减少锯齿。注意,-tap需要原始DEM的投影坐标和裁剪目标一致。
4.4 现象:zip解压后shp加载报错"找不到.dbf"
原因:你用某些国产解压软件解压时,把长文件名截断了,或者漏解出.dbf。shp文件通常由多个附加文件组成,少一个.dbf,属性表就空了,ArcMap直接报错。
解决:重新解压,解压到英文路径,比如D:\data\jiaozuo,不要放在桌面或中文目录下。如果还是报错,用ogrinfo测试shp是否能打开。能打开说明只是ArcMap缓存问题,重启软件;打不开说明shp本身就缺文件。需要向数据提供方重新索要。这里特别提醒:如果你准备只拷贝shp给别人,必须把同名的.shp、.shx、.dbf、.prj四个文件一起拷,缺一不可。很多"数据包坏了"其实是拷贝时漏了.prj,导致坐标系信息丢失。
5. 让30米DEM真正可用:坡度、坡向、山体阴影的派生参数
5.1 在ArcGIS里计算坡度和坡向:参数别选错
有了裁剪后的DEM,下一步通常要做坡度分析。在ArcMap中打开Spatial Analyst工具箱 → 表面分析 → 坡度。输入栅格选DEM,输出单位选"度"(Degree),"Z因子"(Z factor)到底填多少,是个玄学。
Z因子是在X、Y单位与Z单位不一致时,用来换算高程比例的参数。如果你的DEM是UTM投影(单位为米),高程也是米,Z因子设1就够了。但如果你用的DEM是地理坐标系(经纬度),那么X、Y单位是度,高程单位是米,Z因子不能直接设为1。标准做法是先定义投影(Project Raster到UTM),再计算坡度。否则一个纬度约等于111公里,你的坡度值会小得离谱。如果非要用经纬度DEM算坡度,常见做法是把Z因子设为111320(每度对应约111.32公里),但这只适用赤道附近,焦作在纬度35°左右,cos(35°)约为0.819,实际每度经度对应约91公里,更准确的是设Z因子为111320 * 0.819。所以我一般强烈建议:先投影再算坡度,别在Z因子上省事。
坡向工具(Aspect)默认输出0-360度,0代表北,90代表东。它输出的是一个浮点型栅格,你可以在"唯一值"里看到平地区域(-1)。计算后别忘了用"重分类"把坡向分成9类(北、东北、东、东南、南、西南、西、西北、平地),方便后面做光照分析或土地利用。
5.2 用GDAL生成山体阴影和彩色地形
如果你不想开ArcGIS,QGIS和GDAL同样能出效果很好的地形图。GDAL自带一个gdaldem命令,专治DEM可视化。先裁剪,再生成山体阴影:
gdaldem hillshade D:/jiaozuo_dem_30m/dem_clip.tif D:/jiaozuo_dem_30m/shade.tif -az 315 -alt 45 -z 1.2参数解释:-az是太阳方位角,315度表示西北方向来光,这是地形图默认;-alt是太阳高度角,45度比较自然;-z是垂直拉伸系数,可以放大起伏,焦作北部有太行山,如果你觉得山体阴影太平,把-z调到2,效果立刻立体起来。然后用color-relief给DEM上色,先准备一个色带文件color.txt:
300 153 204 255 600 204 255 204 900 255 255 153 1200 255 204 153 1500 204 153 102第一列是海拔,后三列是RGB颜色。执行:
gdaldem color-relief D:/jiaozuo_dem_30m/dem_clip.tif color.txt D:/jiaozuo_dem_30m/color_dem.tif最后在QGIS里加载shade.tif和color_dem.tif,把彩色DEM放在上层,把山体阴影放在下层,设置混合模式为"叠加",就能得到一张像卫星影像一样的地形图。这招在做项目汇报时特别好用,领导通常以为你花了很久做效果图,其实就是两条GDAL命令。
5.3 坡度分级和面积统计:用表格验证数据是否合理
算出坡度后,通常需要统计不同坡度级的面积占比。比如做建设用地适宜性评价,坡度<5%是平坦地,5%-15%是缓坡,>25%不适合建设。在ArcGIS里用"重分类"把坡度分成3级,再用"栅格转面"工具转成矢量,最后用"计算几何"算出每个级别的面积。这里有一个需要注意的点:面积要用投影坐标系下的栅格来算,如果你用经纬度DEM直接算,面积会严重失真。
具体参数:重分类时要确定临界值。焦作市北部太行山区坡度普遍大于25%,南部平原坡度小于5%,中间的丘陵地带在10%-20%之间。你可以先做坡度,再用栅格计算器统计百分比。用Python的arcpy也能直接算:
import arcpy from arcpy.sa import * arcpy.env.workspace = r"D:/jiaozuo_dem_30m" dem = r"D:/jiaozuo_dem_30m/dem_clip.tif" slope = Slope(dem, "DEGREE", 1) # 重分类:0-5, 5-15, 15-25, 25以上 reclass = Reclassify(slope, "Value", RemapRange( [[0, 5, 1], [5, 15, 2], [15, 25, 3], [25, 90, 4]]), "NODATA") reclass.save(r"D:/jiaozuo_dem_30m/slope_reclass.tif")这里RemapRange的每一行是[起始值, 结束值, 新值],注意区间是左闭右开(包含下界,不包含上界),所以5度会进到第一级还是第二级取决于你写的数字。如果想让0-5包含5,应该写[0, 5.0001, 1],避免浮点数临界。这个坑我踩过多次,一度以为分类结果错位,其实是被浮点精度坑了。
6. 最后一步:用shp范围验证DEM边界,和导出KML的取舍
这一章给你一个收尾技巧:用shp转KML,在手机上快速核对DEM覆盖范围。具体做法是在ArcMap的"转换工具"里选择"转为KML",输入shp,输出一个.kmz文件,用手机上的奥维地图或Google Earth打开,就能看到焦作市的边界。但这个操作有个前提:shp必须已经在WGS84经纬度坐标系下,否则转出来的KML位置会偏移。如果你的shp是CGCS2000投影坐标,先投影转换成WGS84,再转KML。
对于DEM本身,个人习惯是最后做一次"检验":用DEM在焦作市市区范围内随便取一个点,测一下高程。焦作市市区海拔大约在80-110米之间,如果你量出来是1000米,数据源肯定有问题。还有一个更实用验证方法:把DEM加载到ArcScene里,垂直拉伸3倍,看北部太行山的山脊线和山脚线是否符合常识。如果山脊线看起来像网格状,说明数据在拼接时用了低质量的重采样,需要谨慎使用。
最后分享一个教训:以前我拿DEM做流域分析,直接用了未投影的经纬度DEM,结果汇水面积差了将近30%,后来才意识到是坡度计算里的投影问题。从那以后,我拿到DEM的第一件事永远是gdalinfo看坐标,再决定要不要投影。你手里的焦作市DEM,如果shp和DEM坐标系不一致,别怕,按第3章的方法先统一坐标系再裁剪,后面所有分析都会顺畅很多。希望帮到你。
本文还有配套的精品资源,点击获取