简介:这份资源是面向地理信息系统学习者与研究者的山西省晋中市三十米分辨率数字高程模型数据包,覆盖全市范围并附带市级行政边界矢量文件,适用于地形分析、水文模拟、规划选址、教学演示等基础应用场景。包内共十二个文件,以GeoTIFF格式的DEM主文件、Shapefile格式的行政区划边界文件为核心,同时包含投影定义、属性表、空间索引、元数据及快速显示辅助文件,整体约77.49MB,便于在常用GIS软件中直接加载。目前已有三百二十八人学习使用,可作为区域尺度地形数据处理的实用样本。资源将高程栅格与边界矢量集于一体,用户既能提取坡度、坡向、高程剖面,制作地形晕渲图,也可结合属性数据进行范围裁剪、面积统计,或用于流域划分、选址叠加等分析练习,不仅提供数据基础,也有助于理解标准地理空间数据文件的组成逻辑。
1. 拿到晋中市30m DEM压缩包,先别急着解压导入ArcMap
做晋中市域的地形分析时,这个“山西省晋中市DEM数字高程数据30m(含本市级范围shp文件).zip”是典型的数据交付包:里面一份30米分辨率的DEM,一份晋中市边界shp,加起来就是一套能直接开工的地形底图。30m分辨率意味着一张像元对应地面上约30米见方,适合全市尺度的坡度、坡向、汇水分析,但别指望拿它做单体工程勘察,精度不够。拿到包以后,多数人的第一反应是解压、拖进ArcMap、直接开干,然后被边界锯齿、投影错位、NoData黑边反复折磨。这篇就围绕这个zip讲清楚一份30m DEM从解压到能用、再到验证的完整路径,也把最容易翻车的细节一次摊开。
2. 解压与文件体检:zip包里到底有什么,能不能直接进分析流程
2.1 解压前的三个确认:目录结构、路径编码、文件完整性
先别双击解压。压缩包在传输和打包过程中最容易出问题的地方不是DEM本身,而是目录嵌套、中文文件名和个别文件损坏。我一般会先在命令行里看一眼包结构,确认里面是“一个tif + 一组shp文件(shp/shx/dbf/prj)”还是套了好几层目录。如果套多层目录,后续无论是ArcMap还是GDAL命令行引用路径都会变长变脆,第一件事就是把它捋平。
# 先列包内容,确认有几个文件和目录层级,避免解压出一堆散文件 unzip -l "山西省晋中市DEM数字高程数据30m(含本市级范围shp文件).zip" # 解压到纯英文路径,比如 D:/gisdata/jinzhong,后续所有工具都省心 mkdir -p D:/gisdata/jinzhong unzip "山西省晋中市DEM数字高程数据30m(含本市级范围shp文件).zip" -d D:/gisdata/jinzhong/unzip -l 的作用是只列出内容不解压,能提前看到文件名是否乱码、是否有MacOS特有的__MACOSX目录、shp各附属文件是否齐全。解压路径一定要避开中文和空格,这看起来是小事,但GDAL在Windows下读中文路径偶尔会直接返回Null,ArcMap也会在要素类路径里把中文识别成非法字符。这个zip名字带着“山西省晋中市”字样,解压到英文目录下正好绕开这个坑。
如果你习惯用Python解压,还要处理一层隐藏编码问题。Windows上很多压缩软件打包中文文件名时用的是GBK,而Python的zipfile默认按UTF-8解出文件名,结果就是文件名变乱码,明明解压成功,文件却打不开。
import zipfile import os zip_path = r"D:/downloads/山西省晋中市DEM数字高程数据30m(含本市级范围shp文件).zip" out_dir = r"D:/gisdata/jinzhong" with zipfile.ZipFile(zip_path) as zf: for name in zf.namelist(): # 先按 cp437 读原始字节,再转成 gbk,解决中文文件名乱码 fixed_name = name.encode("cp437").decode("gbk") target = os.path.join(out_dir, fixed_name) os.makedirs(os.path.dirname(target), exist_ok=True) with zf.open(name) as src, open(target, "wb") as dst: dst.write(src.read())这里的逻辑是:zipfile拿到的name在Windows打包场景下多半是把GBK字节按cp437误解出来的,所以先encode("cp437")还原成原始字节,再decode("gbk")得到真正的中文名。如果包本身是用UTF-8打包的,这步也不会报错,只是decode出来还是原名。转码后逐文件解压,把层级目录自动建出来。实际项目里我还会顺手比对zipfile.namelist和os.listdir的数量,防止静默漏文件。
2.2 用gdalinfo和ogrinfo给数据做体检:不开ArcMap也能知道能不能用
解压完先别急着加载,用GDAL自带的小工具把DEM和shp的元数据读一遍。这一步花不了十秒钟,但能省掉后面在ArcMap里反复“图层属性→源”里翻找的时间。关键是GDAL读到的信息是文件层面的真实值,不会因为ArcMap的默认渲染或坐标系误判而骗你。
# 看DEM:分辨率、投影、NoData、值域 gdalinfo D:/gisdata/jinzhong/dem_30m.tif # 看shp:几何类型、要素数量、四至范围、字段列表 ogrinfo D:/gisdata/jinzhong/jinzhong_city.shp -so -algdalinfo的输出里优先看几个字段:Driver是GTiff还是HFA,决定你在ArcMap里是否需要额外的栅格扩展;Pixel Size如果是0.000277度左右,说明是地理坐标系下的1弧秒DEM,换算到晋中所在纬度约南北30米、东西27米;NoData Value常见是-32768或0,这直接关系到后面做坡度、填洼会不会出现黑边;Computed Min/Max能反映数据是否经过拉伸或分类处理,如果Min是0、Max是255,说明已经渲染成了8位影像,不是原始高程,不能用于坡度计算。投影信息一定不能略过,它决定DEM和shp能不能自动贴合到同一个位置。ogrinfo的输出重点看Geometry是Polygon还是MultiPolygon、Feature Count是否为1、Extent的范围是否和晋中市的经纬度范围对得上。如果shp的Extent显示成8位数的投影坐标,而DEM是经纬度,那坐标系对齐问题就摆在眼前了。
这里还要提一个容易被忽略的前提:30m DEM有可能是SRTM、ASTER GDEM这类开源数据重切片来的,也有可能是基于DSM产品做的地形化处理。DSM和DEM只差一个字母,前者包含了地表建筑和植被高度,后者是裸地面。如果发现gdalinfo输出里数据噪声偏大、等高线在上山地带出现异常锯齿,就要怀疑源头是不是DSM,后面用的时候得先做滤波或重处理,不能直接当DEM用。
提示:gdalinfo和ogrinfo是GDAL/OGR自带命令,安装QGIS或OSGeo4W时会一并装上,不需要单独配环境。如果电脑上已经有ArcGIS但不熟命令行,也可以用ArcMap里的“图层属性→源”看同样信息,只是慢一些。
3. 坐标系对齐:DEM和shp能不能贴合,全看这一步
3.1 先判断坐标系:prj文件缺失时怎么识别
拿到shp后第一件事是看旁边有没有prj文件。shp本质是“一组文件”,主文件存几何,dbf存属性,prj存坐标系。很多网上下载的行政区划shp只给了shp+shx+dbf,缺了prj,ArcMap就默认给它一个未知坐标系,加载后跑到莫名其妙的位置。判断方法很简单:用记事本打开prj文件看WKT字符串;没有prj就用ogrinfo -so -al的输出,如果Extent是类似(113.5, 36.2, 114.4, 37.8)的10以内小数字,那就是经纬度地理坐标系;如果是像(530000, 4020000, 580000, 4200000)这种大数字,就是高斯投影或UTM投影。晋中市大致位于东经111.5°到114.5°、北纬36.5°到38°,如果shp范围落在这一带,坐标系基本可以确认是CGCS2000或WGS84下的地理坐标系或3度带投影。
还有一个玄学但实用的判断技巧:把shp拖进ArcMap,同时加载一个已知坐标系的天地图底图或在线影像,看它贴在山西的哪个位置。如果跑到非洲或太平洋中间,基本就是投影坐标系被错当成地理坐标系使用了。这种错位不是“平移一下就行”的问题,必须先明确shp原本的坐标系定义,再决定是做Define Projection还是Project。
3.2 三种坐标转换方式:ArcMap手动、GDAL命令行、Python脚本
常见做法是把DEM和shp统一到同一个投影坐标系再做裁剪,因为投影坐标系下像元是等距的,面积和距离计算可直接进行。晋中市跨度不大,一般用CGCS2000 3度带,中央经线要看shp的prj里定义的是111°E还是114°E。ArcMap里的操作路径是:先给未知坐标系的shp用Data Management Tools→Projections and Transformations→Define Projection补上坐标系,再用Project工具把矢量转到目标投影;栅格则用Project Raster,重采样建议选Bilinear而不是Nearest,30m DEM经过重采样后值域更平滑。
命令行方式更适合需要反复处理或批量交付的场景。用gdalwarp一步完成投影转换:
# 将地理坐标系的DEM转成CGCS2000 3度带投影坐标 # 示例用EPSG:4547,对应中央经线111°E的3度带;实际以shp的prj为准 gdalwarp -t_srs "EPSG:4547" -r bilinear D:/gisdata/jinzhong/dem_30m.tif D:/gisdata/jinzhong/dem_4547.tifgdalwarp -t_srs指定目标坐标系;-r bilinear指定重采样算法为双线性内插。对于高程栅格,Nearest会保留原始像元值但容易出现锯齿状台阶,Bilinear更平滑,但要注意在陡峭地形上会略微削平峰谷。如果目标只是做山体阴影,Bilinear够用;如果要做等高线,有人坚持用Cubic或Cubic Spline,但平滑算法会压低真实谷底深度,30m数据本身已经不是原始测量值,不必过度追求平滑。转换完成后用gdalinfo再看一遍Pixel Size,确认已经从度变成米,例如接近(30,30)才说明投影正确。
Python方式适合把转换嵌入自动化流程:
from osgeo import gdal src = D:/gisdata/jinzhong/dem_30m.tif dst = D:/gisdata/jinzhong/dem_4547.tif ds = gdal.Warp(dst, src, dstSRS="EPSG:4547", resampleAlg="bilinear") ds = None这里gdal.Warp的参数和命令行一一对应,dstSRS等价于-t_srs,resampleAlg等价于-r。写成脚本的好处是可以对多个县区shp循环处理,也能把后续的裁剪、NoData赋值合并到一个流程里。还有一点容易被忽略:如果shp和DEM坐标系不一致,直接进入裁剪环节GDAL会按shp的坐标系对DEM再做一次隐式投影,输出栅格的像元尺寸可能会变成不规则数值。所以建议无论如何都先把两者统一到同一投影坐标系,再做切割,后面检查也能少一个变量。
3.3 转换后必须做的一致性检查
投影转换不是点一下“确定”就结束的。我见过太多人转换完不检查,最后交出去的数据像元尺寸变成30.0003145米,行数也变了,看起来没毛病,叠加起来就是差一点点。这里给一个快速检查清单:首先看Pixel Size是否仍是规则的(30,30)左右,如果出现(30.00004, 30.00003)这种值,说明输入的原始DEM在投影重采样时产生了微小形变,一般可接受,但如果要拼图或做变化检测就需要统一强制对齐网格;其次看Origin坐标是否为规则的米制数值,如果不是,说明gdalwarp按目标范围自动重新对齐了网格;最后看Min/Max值域是否和转换前基本一致,如果极值明显变化,重采样算法可能引入了异常像元。再打开ArcMap,把DEM和shp同时加载,缩放到市界范围,直接用肉眼确认边界贴合度。如果shp在DEM上浮空或错位,回到前面第一步检查prj,不要试图用“空间校正”硬拽。
4. 用shp范围文件裁剪DEM:ArcMap掩码提取与GDAL cutline的取舍
4.1 ArcMap依靠面图层裁剪DEM栅格tif:Clip和按掩膜提取的本质区别
投影统一之后,就到了整个交付包的核心操作:用晋中市范围shp把DEM裁出来。ArcMap里有两套工具经常被混用,一套是Data Management Tools→Raster→Raster Processing→Clip,另一套是Spatial Analyst Tools→Extraction→Extract by Mask。名字差不多,原理却不完全一样。
Clip工具是几何层面的裁剪,默认输出是裁剪框的外接矩形范围,如果想贴合市界,必须勾选“使用输入要素裁剪几何”选项;它直接按矢量边界切掉外部像元,输出范围贴合边界,但边界外的像元直接丢弃,不会写入NoData。Extract by Mask则是像元层面的提取,用掩膜面去“刷”栅格,落在掩膜外的像元会被写成NoData,输出仍然是矩形范围,边界内外的差别只体现在像元值上。对DEM来说,这个区别直接决定了后续分析时NoData区的形状:Clip出来的数据边缘是矢量边界的,Extract by Mask出来的数据在矩形四角会出现NoData块。
实际操作中我一般优先用Clip并勾选裁剪几何,因为DEM的NoData越少越好,四角的NoData块会让坡度分析时边缘多出一圈无效区。Extract by Mask也有它的用武之地:如果掩膜不是精确的市界,而是一片缓冲区,或者你要对多个面要素分别提取每个面的子DEM,用掩膜提取更顺手。不要把两者当成同一个按钮,选错的结果就是后续每次计算都在边角报“ERROR 010091: Not enough valid pixels”。
“依靠面图层裁剪dem栅格tif”和“依靠面图层掩码提取”的热搜词差别,本质就在这里。前者要的是结果边界贴合,后者更看重像元有效值的筛选逻辑。对晋中市这个zip来说,两种工具都能裁出市域DEM,但如果你后续还要做填洼和水文分析,Clip更合适。
4.2 用gdalwarp -cutline命令行裁剪:更适合批量与自动化
如果你手上的数据不只一个市,或者要写进交付脚本,ArcMap的手工点选就不够了。GDAL的gdalwarp自带矢量裁剪能力,而且-crop_to_cutline参数正好等价于Clip工具里“使用输入要素裁剪几何”的勾选。
# 用晋中市边界shp裁剪DEM,输出贴合市界边界 # 前提:DEM与shp已经统一坐标系 gdalwarp -cutline D:/gisdata/jinzhong/jinzhong_city.shp \ -crop_to_cutline \ -dstnodata -32768 \ D:/gisdata/jinzhong/dem_4547.tif \ D:/gisdata/jinzhong/dem_jinzhong_clip.tif-cutline指定矢量边界文件,支持shp、GeoJSON等GDAL能读的矢量格式;-crop_to_cutline让输出范围严格贴合边界,不写它的话输出还是整幅矩形;-dstnodata给输出栅格指定NoData值,DEM常用-32768。这里有个细节:如果原始DEM的NoData不是-32768,加不加-dstnodata输出都可能把原始NoData保留下来,但显式指定能保证文件头里写入统一的NoData标识,后续做填洼和坡度分析不用再猜。裁剪完最好顺手用gdalinfo看一眼Size和Pixel Size,确认没有因为投影变化把像元变成非整数尺寸。
4.3 用山体阴影和gdalinfo检查裁剪结果
裁剪完不要直接交付,先做两步自检。第一步是数值自检,运行gdalinfo看Statistics里的Min和Max,和裁剪前的DEM对比,值域应基本一致。如果出现极端值,可能是边界像元被插值污染,这类问题在ArcMap的Pixel Inspector里能看到,但命令行更快。第二步是视觉自检,用gdaldem生成山体阴影,然后连边界shp一起加载:
# 生成山体阴影,用来目检边界形态和NoData黑边 gdaldem hillshade D:/gisdata/jinzhong/dem_jinzhong_clip.tif D:/gisdata/jinzhong/hillshade.tif山体阴影的作用是让地形起伏肉眼可见,边界是否贴合市界、边缘有没有一圈黑色NoData、内部有没有出现异常的条带状亮度变化,一眼就能看出来。山体阴影的默认参数是方位角315度、高度角45度,主要看宏观地形,不必纠结参数。如果发现边界有一圈黑,回到上一步检查-dstnodata和原始NoData是否一致;如果边界出现锯齿状白边,说明cutline的shp边界和DEM像元边缘没对齐,常见原因是shp与DEM分辨率差异大,30m像元在边界处只能近似表达,这种情况不影响内部分析精度,但出图不好看,可以考虑先对shp做一点平滑再裁。
5. 避坑指南:DEM处理中最容易翻车的五个细节
5.1 裁剪后边界出现黑边或锯齿
现象:裁剪后的DEM加载进ArcMap,沿市界边缘看到一圈黑色区域,放大后黑边宽窄不一,有的地方甚至呈锯齿状。
原因:绝大多数情况是NoData值不统一。原始DEM的NoData是0或-32768,裁剪时没有指定-dstnodata,GDAL按原文件写了NoData,但ArcMap渲染时按文件头读到的是-32768,便默认显示为黑色。另一种原因是掩膜边界与像元不对齐,边界穿过了像元中心,导致边缘像元要么保留要么丢弃,形成锯齿。
解决:裁剪时统一显式指定NoData,用gdalinfo确认输出栅格的NoData值和原始值一致;如果已裁完,用“栅格计算器”把NoData区域重新赋值为统一值,或用栅格转ASCII后批量替换。锯齿是栅格化边界的固有现象,30m像元在平滑曲线上必然有阶梯感,只要边界像元的中心在shp范围内,就不算数据缺陷,出图时可以用shp边界线遮一下。
5.2 shp缺少prj文件导致坐标错位
现象:DEM和shp明明都在晋中范围,加载后却一个在山西、一个在海外,或者重叠位置差了几百公里。
原因:shp没有prj文件,ArcMap把它当成未知坐标系,默认加载到WGS84经纬度原点上。更危险的是,shp的坐标值本身是投影坐标,却被当成经纬度解释,位置就完全错乱了。
解决:先确认shp坐标值用数字特征判断投影类型,再用Define Projection手动指定正确的坐标系,最后用Project转到与DEM一致的系统。这里强调一句:Define Projection只是修改坐标系标记,不改坐标数值;Project才是真正做坐标变换。很多人拿错工具,结果坐标依然错位。
5.3 坡度计算后NoData区域异常扩大
现象:对裁剪好的DEM做坡度分析,输出栅格的NoData区域比原始DEM大了整整一圈,山脊和谷底出现碎小的无效点。
原因:坡度计算需要中心像元与周围8个像元共同参与,边界上的像元缺少邻居,被算法标记为NoData。这是栅格分析的正常行为,但裁剪时如果NoData值设置不统一,原本有效的边界像元也会被当成无效值处理,无效区就会扩大。
解决:在计算坡度前,先确认DEM的NoData统一为-32768,并且填补内部空洞。如果需要保留边界像元,可以对DEM做一次低通滤波或使用“焦点统计”填充边缘,细节多但能保住边界。实际项目中如果用边缘效应明显的分析,我一般直接把分析范围外扩几个像元,裁到市界后再做坡度,输出再以市界shp裁剪一次。
5.4 没有裁剪几何导致输出范围比shp大一圈
现象:明明用市界shp做了Extract by Mask,输出DEM的四角还留着大片NoData区域,范围仍是矩形,比市界大一圈。
原因:Extract by Mask的默认输出范围是掩膜要素的外接矩形,不是掩膜本身。没有勾选“裁剪几何”的Clip工具也一样。很多人误以为用了面要素就会按面裁剪,但这个“面”只是用来筛选像元,输出范围还是矩形。
解决:确定你需要的交付形态。如果接受NoData填充矩形,直接用Extract by Mask;如果要贴合市界,改用Clip工具并勾选“使用输入要素裁剪几何”,或者用gdalwarp加-crop_to_cutline。两种结果的差别在文件大小和后续分析范围上都有体现,提前想清楚比事后补救省事。
5.5 大面积填洼填出平地,水文流向失真
现象:为了做水文分析对DEM执行Fill工具,结果地形图上出现一片片完全平坦的区域,汇流累计栅格在平地上随机乱流,河道走向和真实水系对不上。
原因:30m DEM里本来就有一些真实或虚假的洼地,Fill工具的默认阈值会把这些洼地全部填平。但如果在黄土高原这类地形起伏大的区域,宽浅的干沟和真实洼地被错误填平,平原部分就直接变成零坡度,流向计算失去依据。
解决:填洼前先检查NoData和空洞,不要直接对整个市域一次性Fill。可以用“填洼”工具的“限制”参数,或者先做两次不同阈值的填洼对比,选阈值更小的一次。对30m DEM来说,阈值设在1到20之间比较常见,具体看研究区的微地形起伏。填完后必须做汇流累积验证,如果河道出现大段平直或断裂,及时回退,这一点放在下一章展开。
6. 验证DEM能不能用:用汇流累积抽检地形连续性
裁剪完成、避坑到位,最后一步是确认这块DEM真的能支撑地形分析。我最常用的验证方法不是看颜色渲染,而是做一次快速的水文分析抽检。流程不长:先对裁剪后的DEM执行Fill填洼,再算Flow Direction流向,再算Flow Accumulation汇流累积。在ArcMap里依次走Spatial Analyst Tools→Hydrology的三个工具即可,参数用默认值。
汇流累积结果出来后,用栅格计算器把累计流量大于1000或2000的像元提取成河流网络,再叠加晋中市范围内的已有河流shp做目检。判断标准有两条:主河道应该连贯,不能出现中断;河道走向与实际水系要大量重合。如果出现长距离平行河道或大块平地上的随机乱流,说明DEM内部还有大量未修复的洼地,需要回到填洼步骤调整阈值。这个验证虽然不能量化到毫米级,但对30m分辨率的数据来说,已经足够暴露绝大多数问题。
注意:这个验证方法对平坦地区尤其敏感。晋中市东部太行山区和中西部盆地地形差异大,如果只验证山区,西部平原因填洼过度造成的水系错乱很容易被漏掉。
验证通过后,这份“山西省晋中市DEM数字高程数据30m(含本市级范围shp文件).zip”才算真正能落地使用。无论是做坡度坡向、生成等高线,还是后续把市界shp转成三维服务的边界,都不会再被底层数据缺陷拖后腿。30m DEM的数据源本身不可能完美,但裁剪、坐标、NoData这些环节做扎实,后面的分析才谈得上可靠。这套流程我从第一次处理同类数据翻车之后就一直沿用,先体检再对齐,最后用汇流累积收尾,已经成了习惯。希望帮到你。
本文还有配套的精品资源,点击获取