★ FEATURED ARTICLE 简介全国DEM数据TIF数据.zip提供覆盖全国范围的数字高程模型栅格数据面向地理信息、遥感及相关地学研究人员与工程技术人员可用于地形特征提取、流域划分、土地利用适宜性评价以及三维地形可视化等多种分析场景。压缩包内共有七个文件核心数据为TIF格式的高程栅格同时附带地理配准文件、投影参数定义、金字塔文件、属性表以及文本格式的补充说明可辅助在ArcGIS、ENVI等软件中实现坐标加载、快速显示和基础查询。数据覆盖完整且带有空间参考信息使用时需确认投影坐标系和分辨率是否匹配任务要求必要时可进行填洼、裁剪和误差校正等预处理。整体压缩包约三十一兆字节数据量适中便于存储和交换。目前已有10778人学习浏览对于需要全国尺度高程底图的区域分析、科研教学与项目开发而言是一份实用且易用的基础数据。1. 全国DEM数据TIF数据.zip拆开之前先想清楚你要拿它干嘛拿到一个“全国DEM数据TIF数据.zip”第一反应通常是解压之后直接拖进GIS软件然后等着图层出现。实际做过的都知道这个操作翻车概率极高要么打开一片黑要么坐标跑到海里要么坡度算出来几十度起步。问题不在数据本身而在大多数公开DEM产品是按分幅组织的原始采样格网并不是为“桌面直接出图”准备的你需要先把它当成原材料走一遍检查、拼接、重投影、裁剪的预处理链。这个zip能解决的是全国尺度的地形底图问题坡度、坡向、山体阴影、水文分析、选址通视、淹没模拟凡是依赖高程栅格的生产任务都从它开始。适合两类人一类是用GIS做分析但被数据组织方式搞烦的从业者另一类是写自动化脚本处理栅格的开发工程师。后者如果直接拿一个分幅TIF就开跑后面会被坐标系和NODATA折磨到怀疑人生。先说结论这份数据包本身不是拿来直接出图的它是“半成品原料”。能不能用得顺取决于你愿不愿意在动手前花十分钟搞清楚它的分幅规则、高程基准和无效值设置。这三件事搞明白后面每一步都有后悔药可吃搞不明白每一步都是黑匣子。2. 先读懂分幅TIF分辨率、高程基准与坐标参考决定后面每一步全国尺度的DEM数据通常不会拼成单个大文件分发而是按经纬度格网或标准图幅切成大量TIF分幅每个分幅覆盖一个矩形区域。这带来一个直接后果你看到的是一个几百个文件组成的zip不是一个“图层”。在动手拼接之前必须先回答三个问题一个像元在地面上到底占多大面积、高程值相对于哪个基准面、每个分幅的坐标参考是什么。这三个问题决定了后面每一步的操作参数。分辨率影响重采样设置基准面影响高程比较坐标参考影响投影转换。很多人拿到数据直接做坡度分析算出来的结果离谱到没法看根因就是没搞懂这三个属性之间的关系。2.1 “30米分辨率”不是你想的那种均匀网格所谓30米分辨率在大多数公开DEM产品里指的是赤道附近的像元尺寸。实现方式通常是在经纬度方向上按固定角秒切分比如1角秒约等于30米。这样在地球高纬度地区南北方向还是30米但东西方向的实际地面距离会明显小于30米因为经线在向两极收敛。换句话说同一份DEM在低纬和高纬地区的地面采样密度并不一致。这个特性直接影响到两件事。第一如果你把经纬度坐标的DEM直接当成平面数据算面积、算坡度高纬度地区会整体变形第二如果同一个项目要拼接不同纬度的分幅你没有必要把每个分幅都重采样到同一个像元尺寸重投影时会统一处理。另外高程值本身还要区分基准有些产品的高程是相对于椭球面的椭球高有些是相对于大地水准面的正高两者差几米到几十米不等。做淹没分析或者土方计算时这几米可能就是致命的误差来源。判断依据一般藏在zip内的说明文件里文件名里带版本号或者源数据标识的要优先按随附元数据来解读。2.2 用gdalinfo在30秒内识别分幅TIF的关键信息不管后续用桌面GIS还是脚本处理第一步都建议用GDAL自带的gdalinfo把分幅TIF的元数据打出来看看。它对GeoTIFF的读取是行业最稳的不会像某些桌面软件那样自动做“看不见的修正”。gdalinfo dem_30m_n28_e105.tif这是命令行下最朴素也最有用的动作。输出信息量很大但真正要读的字段就下面几个我做了个快速对照表。字段含义判断要点Size is 1201, 1201栅格宽度和高度像元数确认分幅行列数是否符合预期Origin (104.9999999, 29.0001389)左上角坐标和文件名里的经纬度对得上才对Pixel Size (0.0002777778, -0.0002777778)每像元对应的经纬度跨度转成米时要除以纬度余弦Coordinate System is坐标参考描述缺这个说明坐标信息丢失要警惕NoData Value -9999无效值标记没有这行说明后面渲染几乎必黑屏Compression / Interleave压缩方式和波段存储有LZW压缩说明分幅设计合理BIGTIFF是否启用大文件格式没看到但你处理全国数据预算好4GB上限这些字段在拼接之前就值得逐一过一遍。尤其是NoData Value如果没有输出后面所有分幅拼到一起后无效区域会被当成真实高程参与计算渲染和坡度结果都会出现大面积异常。遇到缺NoData的TIF用gdal_translate补一个再进流程。2.3 压缩、分块与BIGTIFF加载性能的三个隐藏因素TIF后缀只说明容器格式不说明内部组织方式。同样是TIF有的加载飞快有的拖动一下转半分钟差别就在压缩、分块和文件大小上限三个设置上。LZW或DEFLATE压缩的TIF体积更小读取时GDAL会解压CPU开销略高但磁盘IO压力大幅降低。分幅DEM一般用LZW就够DEFLATE的压缩率更高但解压更慢。TILEDYES意味着数据按小块如256×256组织存储而不是按行扫描这对切片服务和局部窗口浏览非常友好全国范围的大栅格更应该开。至于BIGTIFF是因为经典TIFF格式单个文件超过4GB就无法写入全国DEM拼接后动辄几十GB如果你直接gdal_translate出去而不开BIGTIFFYES文件会写到一半报错或者被迫拆成多个文件。预处理阶段就把这三个参数想清楚后面能省掉大量返工。3. 把全国分幅DEM TIF整理成可用的底图解压检查、虚拟拼接与重投影知道了分幅TIF的底细接下来就是预处理。我一般会遵循一套固定顺序先解压做完整度检查再建立虚拟拼接然后一次性做重投影和重采样最后才根据分析目标去裁剪和派生地形因子。这个顺序不是随便定的。直接对几百个分幅跑gdal_merge.py不仅慢而且内存经常撑不住直接在每个分幅上分别做投影转换又会产生海量中间文件。正确做法是先用gdalbuildvrt“假装拼好”再用gdalwarp一次性物化这样只产生一份大TIF效率最高也最容易排查问题。3.1 解压后先做三件事检查文件命名、配套文件与NODATA解压不是双击完事。打开目录后马上确认三件事分幅文件的命名规律、有没有配套的坐标文件、总数据量是否和zip标注一致。这些信息决定你后面能不能用通配符批量处理。ls -lh dem/ du -sh dem/ls看一眼命名du看一眼总量。命名规律很关键常见的如dem_30m_n28_e105.tif其中n28代表北纬28度e105代表东经105度这种命名可以直接推算出分幅覆盖范围也可以用来判断数据是否存在空洞。配套文件方面如果每个TIF旁边有.tfw、.prj或者同名xml说明坐标信息齐全如果只有一个裸TIF说明坐标信息可能内嵌在GeoTIFF头里也可能已经丢失需要用gdalinfo逐一确认。我见过最坑的情况是zip里几百个文件解压后有一半分幅的.tif文件缺失但没有任何报错提示直到拼接完成后才发现北部区域少了一大块。所以进入下一步之前最好用文件名正则把该有的分幅清单列出来和实际文件做一次差集。这步不复杂但能省下后面排查位置错漏的大把时间。3.2 用gdalbuildvrt做虚拟拼接避免一次合并全国数据分幅数量多的时候直接合并不仅慢而且一旦参数不对就要重新跑一遍。虚拟拼接VRT是更好的中间状态它不复制任何像元数据只是记录“哪些TIF、按什么规则叠在一起”的XML描述。加载时GDAL会按需去读原始文件速度取决于原始文件的读取效率。gdalbuildvrt -o national_dem.vrt dem_*.tif这条命令把dem目录下所有TIF按空间位置拼进一个VRT。-o指定输出文件名dem_*.tif是shell通配符会展开成所有分幅。如果分幅之间有重叠VRT默认按顺序取第一个文件的数据所以d_e排序要留意——建议命名规律能保证按行列顺序排列。遇到个别分幅坐标缺失时VRT会失败并提示找不到坐标信息这就是前面检查配套文件的价值所在。VRT还有一个好处是方便做临时验证你可以直接加载VRT进QGIS看全局渲染效果不需要等物化完成。发现问题时删掉VRT重做即可根本不伤原始分幅。3.3 坐标与投影重采样gdalwarp的五个必设参数VRT本身不解决坐标系问题。全国DEM原始数据通常是WGS84经纬度坐标而地形分析需要平面坐标下的距离和面积才有意义。这一步用gdalwarp把整个VRT重投影到目标投影同时完成像元重采样直接输出一份正式TIF。gdalwarp -t_srs EPSG:32650 -tr 30 30 -r bilinear -dstnodata -9999 \ -co COMPRESSDEFLATE -co TILEDYES -co BIGTIFFIF_SAFER \ national_dem.vrt national_dem_utm.tif五个关键参数拆开来看。-t_srs指定目标坐标系EPSG:32650是UTM 50N投影适用于东经117度到123度之间的区域你在其他区域要换成对应的分带。-tr 30 30表示重采样后的像元尺寸为水平方向和垂直方向各30米重投影后像元不再以度为单位。-r bilinear选择双线性插值地形数据用bilinear比nearest平滑比cubic快是兼顾效率和精度的选择。-dstnodata -9999把无值区域统一写成-9999防止不同分幅各自为政。最后的-co参数控制输出格式DEFLATE压缩、TILED分块、BIGTIFF按需开启这三件套在大范围栅格里基本属于标配。有一点值得反复强调分析用的投影别选Web墨卡托EPSG:3857。墨卡托在高纬度会把面积拉大好几倍坡度、面积、通视分析全部失真。做省级或更小尺度分析时用高斯-克吕格或UTM分带投影都行具体分带要看项目区经度范围。这个选择比任何花哨算法都重要。3.4 什么时候该把VRT“物化”成正式TIF有人会在VRT上直接做后续的裁剪和坡度计算这当然可以但性能很差。VRT每次读取都要临时去解析原始文件的位置和格式反复分析时等于把IO压力放大好几倍。我一般会在重投影完成后就把结果物化成正式TIF也就是上面那条gdalwarp命令输出的national_dem_utm.tif。之后的所有裁剪、坡度、填洼操作都基于这份文件。这样做的另一个好处是这份TIF已经是统一坐标、统一分辨率、统一NODATA的“干净底图”分发给别人或做缓存都用它不会因为原始分幅存在差异而产生理解偏差。4. 全国DEM栅格生产流程正式拼接、按行政区裁剪与地形因子计算底图准备好之后进入生产阶段。这个阶段的任务很明确把全国底图变成真正用于分析的成果——拼接好的全域栅格、裁剪到工作区的DEM、以及坡度、山体阴影、填洼后的水文修正栅格。每一步都有对应的GDAL工具或轻量Python库参数不复杂但值得逐个确认因为地形成果不像普通遥感影像“看得见就算对”它需要数值尺度和空间范围都正确。4.1 用gdal_merge.py做正式拼接BIGTIFF与Float32VRT是虚拟拼接物化后的正式拼接用gdal_merge.py最省事。它的逻辑是把多个输入栅格按空间位置合并到一块画布上重叠区默认取第一个文件的值。如果你的分幅之间严格无重叠这个工具几乎零成本。gdal_merge.py -o dem_merge.tif -ot Float32 -a_nodata -9999 \ -co BIGTIFFYES -co TILEDYES -co COMPRESSDEFLATE \ national_dem_utm.tif-ot Float32把输出数据类型强制为32位浮点确保高程精度不丢失尤其是后续要做浮点运算的坡度或填洼。-a_nodata统一无效值。这里注意如果输入文件已经做过gdalwarp的统一NODATA设置-a_nodata可以省略但写上更保险。输出的-co和前面一样用三件套全国范围的数据量非常大BIGTIFFYES是必不可少的否则超过4GB直接报错。4.2 用gdalwarp按行政区裁剪cutline的坐标一致性把全国底图裁剪到项目区属于高频需求。最稳妥的做法不是“按范围框裁”而是“按矢量边界裁”避免把边界外的高程带进来。gdalwarp -cutline county.shp -crop_to_cutline -of GTiff \ dem_merge.tif county_dem.tif-cutline指定行政区矢量文件-crop_to_cutline让输出栅格的范围完全贴合矢量边界而不是保留矩形外框。这两个参数组合起来就等于“无锯齿裁剪”。一个常见的坑是矢量边界坐标系与栅格不一致时GDAL会做动态转换但转换精度取决于矢量本身的精度。建议在裁剪前先确认两者的坐标描述若一个WGS84一个国测局坐标先统一再裁不要指望gdalwarp自动收拾残局。边界处的高程值有时会出现异常条纹尤其是沿切割线的像元。可以在裁剪后加一道腐蚀处理把边界外一圈NODATA写入。做法不复杂用gdalwarp的-dstnodata再跑一次但这一步我通常留到坡度计算之后再做避免过早丢边。4.3 地形因子计算slope、hillshade与参数选择有了裁剪后的DEM坡度、坡向、山体阴影这类派生栅格就是一条命令的事。gdaldem是专门干这个的工具参数不复杂但容易因忽视比例因子而翻车。gdaldem slope county_dem.tif county_slope.tif -s 1 -alg ZevenbergenThorne gdaldem hillshade county_dem.tif county_hillshade.tif -z 2 -az 315 -alt 45gdaldem slope里-s 1表示高程单位与水平单位都是米如果高程是英尺这里要填0.3048。-alg ZevenbergenThorne是可选算法适合相对平滑的自然地形默认的Horn算法在陡峭区域更稳健但没有绝对优劣。坡度输出以度为单位范围0到90。gdaldem hillshade的-z 2是垂直夸张系数山地地形用1或2都行-az 315是光源方位角315度对应西北光-alt 45是光源高度角。这套参数出来的渲染效果比较接近传统纸质地貌图的立体感。坡度计算前一定要确认输入栅格已经是平面投影坐标否则水平单位是度高程单位是米算出来的坡度值每个都像“悬崖”。这个坑在后面章节会展开讲但这里先记住投影先行坡度后算。4.4 水文分析前的填洼GDAL不提供用pysheds补上水文分析是DEM最重要的应用方向之一但它有一个前置动作是GDAL本身没提供的填洼。真实地表有很多凹陷区域如果不处理水流会在这些地方积聚而无法继续向下游流动流向网络和汇流累积量都会失真。所以做水文分析前必须先把DEM中的洼地填平。常见做法是用开源库pysheds它在处理流向、填洼、汇流累积上比桌面GIS更轻量也更容易嵌入脚本流程。from pysheds.grid import Grid grid Grid.from_raster(county_dem.tif) dem grid.read_raster(county_dem.tif) filled grid.fill_depressions(dem)这段代码先用Grid.from_raster读取DEM文件完成网格初始化然后read_raster把高程数组载入最后fill_depressions返回填洼后的数组。填洼结果的数组可以直接传给D8流向计算也可以另存为TIF看后续需求。pysheds的fill_depressions默认会把所有内部洼地填到最低出水口高度参数很少基本不用调。注意填洼只用于水文分析不要为了做坡度或渲染去做填洼它会抹掉真实的地形凹陷信息让坡度图和山体阴影变得失真。5. DEM数据落地五个常见坑现象、原因与解决办法从拿到zip到跑出地形因子这中间每一步都有对应的坑。我按踩坑频率找几个最高发的列出来每条都先说现象再给原因最后是解决办法。这一章可以直接当排查手册用。现象根因一句话解决TIF拖入GIS一片黑NoData被当成有效值参与渲染检查并重写NoData坡度值高得离谱经纬度坐标直接算坡度先重投影再加比例因子水文线在湖面断开水体区域高程是空洞填洼或先补洞拼接时内存溢出一夜白跑输出超4GB且未开BIGTIFF先buildvrt再warp拼好后有明显接缝台阶分幅间高程基准不一致检查数据源统一基准后重拼5.1 打开一片黑NODATA没被识别现象分幅TIF在GIS里加载后图层区域大面积是全黑或全白拉伸符号化时影像看起来像一张噪点图。很多人这时候会认为是文件坏了重新去解压一份结果一样。原因TIF的NoData值没有被GIS读取到。常见的NoData值是-9999或-32768如果这个标签丢失软件会把无效区域的像元当成正常的高程值参与统计与拉伸整个色带被无效值霸占。尤其是高程类型为Int16时-9999看起来更像真实值黑屏概率更高。解决用gdal_translate重写一次NoData标签具体做法是gdal_translate -a_nodata -9999 -co COMPRESSDEFLATE input.tif output.tif。如果你的数据源本身没有NoData概念却用0表示海洋或空洞那就把0也设成NoData但前提是你确认0不是真实地形。5.2 坡度值离谱直接在经纬度坐标上算slope现象从同一份TIF跑出来的slope绝大多数像元值直接在50度以上甚至接近90度整张坡度图全是红色。换一个分幅还是这样怎么看都不对。原因输入栅格是经纬度坐标水平单位是度高程单位是米两者在数值上差了约10万倍。gdaldem按单位一致处理自然把每个像元的高差当成在“度”的水平尺度上变化坡度值必然离谱。解决回到第3章先用gdalwarp重投影到平面坐标再跑坡度。如果因为某种原因必须保留经纬度坐标gdaldem提供-p参数来标记地理坐标系同时需要配-s设置合理的尺度因子但不同GDAL版本对这个参数的默认行为有差异我一般直接重投影不跟它赌。5.3 河流断流水体空洞让水文分析失效现象水文分析生成的河流网络经过某个大湖或水库区域时断成几截或者河流在某片平坦区域突然消失整个流域边界也是碎的。原因水体在原始DEM里往往没有有效高程值要么是NODATA要么是填充的固定值。这些区域在水文分析时相当于“无数据区”水流经过就断了无法继续向下游累积。解决做水文分析前先做填洼和补洞。填洼用第4章的pysheds代码即可但如果空洞区域过大比如一大片水面占了几千个像元填洼可能不够需要用邻域插值或外部水体边界数据把水面高程补上再把水面区域的流向设置为强制排水。检查方法是看TIF的NODATA占比超过5%就该重视。5.4 内存溢出一上来就merge全国分幅现象用gdal_merge.py直接拼几百个分幅TIF跑了两三个小时最后抛出一个“iomem”或“gdal_merge.py: out of memory”之类的错误前面时间全白费。原因gdal_merge.py在处理大量输入时会在内存里同时维护多个输入文件的数据块窗口。分幅多、窗口尺寸大时内存占用量直线上升再加上输出文件超过4GB时没有BIGTIFF支持写入直接失败。解决先gdalbuildvrt建立虚拟拼接再用gdalwarp带-co BIGTIFFYES -co TILEDYES输出正式TIF。这个组合会让GDAL按块读写内存占用可控。如果单次输出仍然太大可以按纬度带分批warp再合并但一般情况下buildvrtwarp已经足够。5.5 拼接处接缝断层分幅间高程基准或时相不一致现象拼好的DEM在山体阴影图里能看到明显的条带拼接处一边高一边低形成台阶状纹理。放大后看接缝两侧的高程值有系统性的几米到几十米差异。原因分幅数据来自不同批次甚至不同源数据。高程基准、像元中心对齐、时相带来的地表变化都会在接缝处产生不连续。DEM是浮点栅格这种差异不会像影像拼接那样因为羽化而自动消隐。解决先把每个分幅的高程统计信息打印出来看接缝两侧的均值差异。如果差异是固定值用gdal_calc.py做差值修正再拼接如果差异没有规律说明是不同源数据别硬拼回到数据源层面确定处理方案。接缝处坚持用nearest重采样可以减少过渡带的异常值插值但根本解法永远是保证源数据一致。6. 把预处理变成流水线一个DEM质量体检脚本拿去改改就能用处理到这一步相信你已经意识到整个DEM预处理链中最花时间的不是跑命令而是检查每一批分幅数据是否货真价实。第2章我教你用gdalinfo人工看一个文件但到了全国数据几百个分幅的场景人工根本看不过来。我自己现在接到新的DEM数据包第一件事就是跑一个批量体检脚本把每个分幅的关键指标打出来几秒钟就能判断这份数据能不能用。from pathlib import Path import rasterio import numpy as np for tif in sorted(Path(dem).glob(*.tif)): with rasterio.open(tif) as src: data src.read(1) nodata src.nodata if nodata is not None: nodata_mask data nodata else: nodata_mask np.isnan(data) valid data[~nodata_mask] print( tif.name, fsize{src.width}x{src.height}, fres{src.res[0]:.2f}, fnodata{nodata_mask.mean()*100:.2f}%, fmin{valid.min():.0f}, fmax{valid.max():.0f}, fcrs{src.crs}, )这段脚本遍历dem目录下所有TIF逐个读取波段数据统计分辨率、NODATA占比、有效值最小最大和坐标参考。判定规则很简单NODATA占比超过5%的分幅要重点检查高程min低于-500或者max超过9000基本可以判断是异常值或坐标错乱res出现两种以上数值说明分幅间分辨率不统一crs为None的分幅直接标红。脚本里src.read(1)会把整个分幅读入内存对几百MB的大分幅会有IO压力。如果文件特别大可以改成只读取元数据和采样统计但当前DEM分幅一般几十MB直接读没压力。跑完一遍输出后把异常项揪出来先处理再走后续的buildvrt和warp整个流程的返工率能降一半以上。我自己早年处理第一批全国DEM时就没做体检直接buildvrt后跑重投影最后出来的底图在某个省份区域出现了整片空洞排查了很久发现是一个分幅的NoData值和其他的不一样。从那以后任何DEM数据拿过来都先让脚本过一遍再做拼接。体检不花多少时间但能避免你在第5章那些坑里反复横跳。希望帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
