简介吉林省吉林市30米分辨率DEM数字高程数据包适合GIS学习者、城乡规划从业者及从事环境分析、地质灾害评估、遥感图像匹配等工作的技术人员使用。数据以TIFF格式存储核心高程栅格并附带覆盖吉林市本级的Shapefile矢量边界配合PRJ投影信息、TFW坐标参考以及XML/SBN/SBX等索引与元数据文件可在ArcGIS、QGIS中直接加载并开展地形分析、坡度坡向提取、可视域分析与地图制图等工作。压缩包共12个文件整体大小约81.88MB结构清晰既有高程模型也有行政边界能保证多源数据正确叠加与裁剪。目前已有712人浏览学习适合需要获取吉林省吉林市地区高精度地形底图的专业人士和学生作为基础地理数据使用。借助配套的矢量边界文件可快速完成市域范围的DEM掩膜提取与后续空间统计分析为区域规划与灾害风险评估提供可靠的数据支撑。1. 拿到吉林市DEM zip包第一件事不是解压“吉林省吉林市DEM数字高程数据30m含本市级范围shp文件.zip”这个文件名已经把三个关键信息贴在脸上30米分辨率的数字高程模型、zip压缩包交付格式、一个市级行政边界shp。大多数人会先双击解压然后直接把tif拖进地图软件里看等到发现裁出来的范围对不上、行政边界歪到隔壁市才回头找坐标系的麻烦。接下来按拿到包之后最稳的操作顺序走一遍验证压缩包完整性与文件结构、对齐DEM和shp的坐标系、用市级边界做像素级裁剪、从DEM派生坡度坡向等地形因子、最后把成果压缩交付。这套流程对ArcGIS、QGIS、纯Python/GDAL用户都适用区别只在命令写法上。2. zip包检查文件结构、数据完整性与坐标系对齐2.1 先看压缩包内部结构用zipfile列出每一个文件用unzip -l可以不解压就列出zip内容但我更常写一段Python的zipfile脚本。原因是这一步顺手把每个文件的压缩比、CRC完整性一起查了省得解压到一半才发现tif损坏。from zipfile import ZipFile from pathlib import Path pkg Path(吉林省吉林市DEM数字高程数据30m含本市级范围shp文件.zip) with ZipFile(pkg) as zf: for info in zf.infolist(): ratio info.compress_size / info.file_size if info.file_size else 0 print(f{info.filename:55s} 原始 {info.file_size / 1024:9.1f} KB f压缩后 {info.compress_size / 1024:8.1f} KB 压缩率 {ratio:6.1%}) bad zf.testzip() print(损坏的文件:, bad if bad else 无包完整)testzip()会逐个读取密封项做CRC校验能提前发现下载中断导致的文件内部错误。如果zip这一层损坏后面GDAL打开tif时往往报“不支持的TIFF压缩类型”这类误导性信息排查方向就偏了。压缩包里一般有一个tif格式的DEM栅格文件和一个shp边界文件外加readme.txt或元数据xml这几个文件按表2-1区分处理顺序。zip内文件类型作用拿到后的处理xxx_dem.tif高程栅格30m格网后续全部操作的主数据市级范围.shp矢量边界用于裁剪先看坐标系再决定是否转换.prj / .xml投影定义与元数据与tif的坐标系比对readme.txt数据来源与说明确认数据是否原始裸地DEM而非DSM提示把解压后的tif和shp都改成纯英文文件名再开始操作Windows下GDAL读取中文路径偶尔会失败报错又看不懂纯属浪费时间。2.2 用gdalinfo看tif头信息分辨率、nodata与概览拿到栅格后先用gdalinfo读一遍头信息它输出的信息比ArcGIS属性框更直接。解压后把tif重命名为jilin_dem_30m.tif执行gdalinfo jilin_dem_30m.tif重点看Size is行列数、Origin和Pixel Size三行。30m数据有两种常见存储一是Pixel Size为 (0.00027778, -0.00027778) 这样的经纬度度数二是已经投影成 (30, -30) 米。这两种情况在算坡度时处理方式完全不同度数坐标必须先投影成米制再计算否则坡度公式里分母单位是度算出来全是废值。还要确认NODATA是不是0或-9999水体区域如果填了0直接进坡度计算会把湖面变成大片平坦假地形。用rasterio读同样的信息更方便接进后续Python工作流import rasterio as rio with rio.open(jilin_dem_30m.tif) as src: print(坐标系:, src.crs) print(尺寸:, src.width, x, src.height) print(像素尺寸:, src.transform.a, src.transform.e) print(nodata:, src.nodata) print(边界范围:, src.bounds) print(概览比例:, src.overviews(1))概览比例如果为空说明tif没有预生成金字塔后面裁剪时每次读取全分辨率数据会比较慢可以先用gdaladdo或转COG解决。另外要看tif到底是DEM还是DSM。DEM是裸露地表高程DSM包含建筑物和树冠。30m级的DSM在城市区域会严重抬高地表文件名带DEM的包一般已经处理过但遇到可疑的尖峰值分布建议用gdalinfo -hist看一眼高程直方图城市区域出现大量突兀的非连续高值就要警惕。2.3 shp与DEM坐标系不一致时的两种对齐方案shp的坐标系写在同名.prj文件里GIS软件读取shp时会自动加载。用geopandas读取shp并打印crs就能和tif的坐标系做对比import geopandas as gpd shp gpd.read_file(吉林市范围.shp) print(shp.crs)吉林市跨域范围较大包里的shp大概率是WGS84或CGCS2000地理坐标。遇到tif和shp坐标系不一致常见做法有两个方向把tif重投影到shp坐标系或者把shp转换到tif坐标系。多数情况下我选后者因为裁剪不改变栅格像素网络不会触发插值完整保留原始高程值而重投影tif会让30m格网边缘产生拉花和重采样误差。代码上就是一行shp shp.to_crs(src.crs) # 用裁剪目标tif的crs作为统一基准注意EPSG:4490和EPSG:4326都属于CGCS2000/WGS84体系差异小于一个像素直接统一就行。但如果tif是老西安80或北京54数据就不能用简单crs变换需要七参数转换。动手裁剪前先打印shp的total_bounds和tif的bounds确认两者相交面积占shp的95%以上再继续这一步能避开“裁出来全是黑”的一半原因。3. 用吉林市市级shp裁剪30m DEM跑通像素级mask3.1 geometry_mask完整流程读取shp、统一坐标、numpy掩膜写出把shp作为裁剪掩膜tif作为输入输出只保留吉林市行政区范围内的DEM。核心不是代码量而是理解geometry_mask返回值的含义out_shape对应栅格行列数transform直接取自tif本身shp几何在传入前必须已经统一到tif的坐标系。import numpy as np import geopandas as gpd import rasterio as rio from rasterio.features import geometry_mask dem_path jilin_dem_30m.tif shp_path 吉林市范围.shp # 读取shp并统一到tif的坐标系 boundary gpd.read_file(shp_path) with rio.open(dem_path) as src: if boundary.crs ! src.crs: boundary boundary.to_crs(src.crs) # 生成掩膜maskTrue表示像素中心落在shp内部 mask geometry_mask( list(boundary.geometry), out_shape(src.height, src.width), transformsrc.transform, all_touchedTrue, invertTrue, ) # 行政区域外的像素全部置为nodata dem src.read(1).astype(float32) nodata src.nodata if src.nodata is not None else -9999 dem_clipped np.where(mask, dem, nodata) profile src.profile.copy() profile.update(dtypefloat32, nodatanodata, compresslzw) with rio.open(吉林市_dem_clip30m.tif, w, **profile) as dst: dst.write(dem_clipped, 1)这段代码里最容易出错的是invertTrue。默认的geometry_mask返回True给几何外部设置invertTrue后True表示几何内部这样np.where(mask, dem, nodata)才能正确保留市区像素。另外dem src.read(1).astype(float32)这一步把原始int16转成了float32目的是让nodata-9999不会和真实高程0值混淆。整个吉林市范围的30m网格float32数组大约占用300MB内存普通电脑跑得动但若换成全省数据就必须分块处理。裁剪后检查输出范围直接用gdalinfo看Size is和Corner Coordinates确认外接矩形没变、内部的nodata区域占据了边界外空白区。更可靠的验证是统计非nodata像素数乘以900平方米30m×30m和已知市域面积对比误差在1%以内算正常。超出太多就说明坐标系统一或mask生成有问题。3.2 all_touched与临界像素边界像素怎么取舍geometry_mask的all_touched参数值得单独讲。默认False表示只有像素中心落入shp内部的格子才会被保留True表示只要像素与shp边界有接触就保留。对30m栅格两者的视觉差异只在边界的一圈像素上但下游分析结果差很多做坡度和水文分析建议用True边界多保留一圈能保证山谷和水系在边界处不截断做面积统计和出图建议用False避免边界外多出来的邻市像素混进统计值。边界像素的另一个典型问题是shp边界与tif边缘重合不准裁出来的高程出现沿海岸线或沿边界的锯齿条带。这通常是shp几何做了简化或缓冲造成的解决办法是先把shp栅格化用与tif相同的transform生成像素级掩膜再相乘这样不依赖像素中心与几何边界的空间关系判断gdal_rasterize -burn 1 -init 0 -ot Byte -of GTiff 吉林市范围.shp 吉林市_mask.tif gdal_calc.py -A jilin_dem_30m.tif -B 吉林市_mask.tif --outfile吉林市_dem_calc.tif \ --calcA*B --NoDataValue-9999栅格化方案在纯命令行环境里最顺手不依赖Python的geopandas逻辑也更直白shp内部burn成1tif和mask逐像素相乘。注意这套方案里tif的nodata若本来就是-9999相乘后需要复查行政区内是否混入了由于原nodata导致的异常乘积0值。做完后按shp外接矩形用gdal_translate -projwin切一下减少后续处理的数据量。3.3 裁剪失败或结果全黑按顺序查三个检查点裁剪输出全黑或全是nodata时先不要怀疑mask代码按顺序排查三处。第一shp读进来是不是真正的Polygon有的“边界”是用LineString做的轮廓线geometry_mask不接受非面几何。第二boundary.crs ! src.crs的比较是严格的对象比较EPSG:4490和EPSG:4326即使坐标几乎一致也判定不等转换后如果to_crs没生效shp的位置会完全错开。第三shp的total_bounds和tif的bounds是否真的相交两者范围写反或坐标系统一方向写反时geometry_mask不会报错只是返回全False。剪出来有横向条纹时用shapely的make_valid处理shp自相交再重新跑一遍from shapely.validation import make_valid boundary.geometry boundary.geometry.map(lambda g: make_valid(g) if not g.is_valid else g)这套检查顺序能覆盖九成以上的裁剪失败场景。剩下的尾数原因几乎都是数据本身问题比如tif里有NaN、源DEM缺失值写成了32767这种极大值需要回到数据头信息做直方图检查。4. 从DEM派生坡度坡向30m数据算出来的地形因子怎么用4.1 为什么先填洼再算坡度填洼参数怎么定DEM拿到手直接算坡度是新手的典型错误。30m分辨率的DEM里既有真实洼地也有采集噪声不填洼直接算平地区域会出现坡度接近0的假平地夹杂尖点山谷地带的坡向方向乱成一团后续日照分析或汇水分析基本不可用。填洼在ArcGIS Spatial Analyst里对应Fill工具开源生态里用richdem库最直接import richdem as rd dem rd.LoadGDAL(吉林市_dem_clip30m.tif) rd.FillDepressions(dem, epsilon0, topologyD8) rd.SaveGDAL(吉林市_dem_filled.tif, dem)epsilon0表示把所有洼地填平到水流能顺畅流出的程度适合水文分析。做道路和场地平整分析时我只填明显低于周围的地形坑把epsilon调到5到10这样能保住人工沟渠和小型汇水洼地不被填平。吉林市的地形特征是河谷与山地并存市区沿松花江河谷的高程在180米到250米左右东部山地区域可以超过1400米这种大起伏区域填洼要保守填过头会压掉真实山谷形态。填洼完成后对比高程标准差下降超过5%就要检查参数。richdem在Windows下直接conda install -c conda-forge richdem安装比裸pip稳。4.2 gdaldem的坡度、坡向、山体阴影参数与输出对照填洼完成后用gdaldem一口气生成坡度、坡向、山体阴影三个文件# 坡度-p表示以度为单位输出s1表示水平垂直单位一致 gdaldem slope 吉林市_dem_filled.tif 吉林市_slope_deg.tif -p -s 1.0 # 坡向0-360度平地标记为-9999 gdaldem aspect 吉林市_dem_filled.tif 吉林市_aspect.tif -zero_for_flat # 山体阴影光源来自西北方向适合东北地区地貌渲染习惯 gdaldem hillshade 吉林市_dem_filled.tif 吉林市_hillshade.tif \ -z 2.0 -az 315 -alt 45-s参数是最容易踩的坑当DEM还是经纬度坐标时水平单位是度、垂直单位是米必须先把DEM投影成高斯或UTM再计算或者在-s里填入每度对应的米数。吉林市位于东经126度附近最稳的做法是回到第2章把tif和shp统一投影到中央经线126度的CGCS2000高斯-克吕格投影让坡度计算的水平垂直单位都是米-s 1.0才成立。山体阴影的-az 315是光源方位角、-alt 45是太阳高度角这两个参数纯出图需求时按需调。DOM影像叠加场景里把山体阴影以50%透明度叠在正射影像上能直接看出高差变化和地形断裂带。坡度输出后我一般会做一次分级0-2度、2-6度、6-15度、15-25度、25度以上五级对应建筑选址和道路建设的坡度门槛。对吉林市东部和南部的山地区域25度以上的坡基本不适合作为建设用地这张分级图可以直接进规划分析报告。4.3 用zonal统计算各县级单元的平均高程与坡度shp里如果已经按县级区划分了多个面可以用zonal statistics把每个面的平均高程、最大最小高程、分位数统计出来再合并回shp属性表。rasterstats是最短实现路径from rasterstats import zonal_stats stats zonal_stats( 吉林市_dem_clip30m.shp, 吉林市_dem_filled.tif, stats[mean, max, min, percentile_90], nodata-9999, prefixdem_, ) print(stats[0])percentile_90比max稳健得多能反映一个县级单元内真实的高值地形区而不被单一点位异常值带偏。rasterstats内部对每个面都做一次全图裁剪统计如果shp有十几个面效率偏低更快的做法是先把tif按每个面的外接矩形切小再统计。统计结果用gpd.GeoDataFrame.merge合并回shp属性表就能直接做分级设色图。沟谷高程差大的区域主要看mean和percentile_90的差值差值超过300米说明这个区域地形切割强烈对基础设施选址有直接提示作用。5. 把裁好的DEM与shp打包交付压缩率、COG转出与能落地的周边技巧5.1 将裁好的DEM转为COG再打zip包裁剪完的tif经常比原始包还大原因是行政边界外被写入了大量nodata如果输出时没设置压缩或者nodata被填成了0打包体积会很难看。我一般会把最终成果转成Cloud Optimized GeoTIFF内置金字塔后续发布到GeoServer或Web端都不用再单独建overviewgdal_translate -of COG 吉林市_dem_clip30m.tif 吉林市_dem_cog.tif \ -co COMPRESSDEFLATE -co OVERVIEWSIGNORE_EXISTING \ -co RESAMPLINGNEARESTRESAMPLING必须用NEAREST双线性或三次卷积会让山顶高程在概览层被平滑掉量测时出现5到10米的偏差。DEFLATE对DEM这种连续渐变数据压缩率通常比LZW高文件更小。最后打包时用zip -9压紧zip层虽然对已压缩的tif收益有限但tif是DEFLATE转出来的前提下再收几个百分点也是赚。交付zip的同时在旁边生成一个.sha256文件让接收方校验这比口头说“我传完了”靠谱得多。5.2 打包前用gdalinfo和渲染做正确性抽检交付前的抽检内容固定四样tif能正常打开、nodata排除后最小值正常、坐标系没变、渲染不出现整片黑色。gdalinfo -stats 吉林市_dem_cog.tif | head -n 30-stats输出里如果Minimum显示-9999说明nodata没有被统计排除加载的人会把-9999当地形高点渲染出一座黑色假山。更直观的验证是直接用QGIS拖进去看一眼把shp也一起加载确认市界和DEM边界贴合松花江沿线没有飘出边界的大块黑色区域。这些验证必须在压缩前做完zip包一旦交付接收方解压后发现问题再回传时间成本高出一截。5.3 两个实用技巧外边界简化和shp属性转txtzip里附带的本市级shp后续还要继续用交付前可以对shp外边界做一次简化去掉与DEM对齐时产生的碎小毛刺让出图更干净。简化容差取10米正好是30m DEM单像素的三分之一左右视觉上基本无感但几何体量能小不少from shapely.geometry import shape import json with open(吉林市外轮廓.geojson, r, encodingutf-8) as f: feat json.load(f) geom shape(feat[geometry]).simplify(10.0, preserve_topologyTrue)如果接收方只需要边界坐标点列表不装GIS也能拿数据干活直接用ogr2ogr把shp属性表转成文本格式省去对方装ArcGIS再导出的环节ogr2ogr -f CSV 吉林市范围.txt 吉林市范围.shp最后交付时zip里除了tif、shp和验签文件还应附一份简短readme写清楚数据源类型、坐标系、填洼参数、裁剪时间和nodata值。收包的人不打开GIS就能判断数据能不能直接用。30m DEM本身不是高精度地形数据但配齐投影信息、裁剪边界和派生因子后它足以支撑市级尺度的坡度分级、流域分析和三维地形展示。本文还有配套的精品资源点击获取