☰
全国地形地貌地质数据集处理:从DEM到坡度分析的完整实践
2026/10/3 5:04:36 网站建设 项目流程

简介:面向GIS专业人员、科研工作者及地理爱好者的全国地形地貌地质数据集,整合了全国地质图、地质矿产分布、各省耕地面积、水系、地貌、地形、森林分布、土地利用、土壤类型与植被分布等十类专题数据,可用于地质勘探、灾害风险评估、资源规划、环境保护、粮食安全分析以及城市与交通建设等研究场景。压缩包共609个文件,约772.86MB,以ADF栅格数据、TIF高程影像、SHP矢量要素为主,辅以DAT、DBF属性表及PRJ投影信息等,目录结构清晰,可直接在ArcGIS、ENVI中加载并开展空间分析、地图制图与模型构建。已有4129人学习下载。其中全国地质图与矿产地层信息可辅助地震带判别和找矿部署;DEM地形数据可计算坡度、坡向及高程剖面;土壤、植被和土地利用图层则可支撑农业生产区划、碳汇估算、水土保持与生态保护研究。整套资源兼顾自然地理与人文地理要素,既是区域规划、学科教学和论文写作的实用底图,也是工程选址、灾害防治与可持续发展研究的综合数据基础。

1. 全国地形地貌地质数据集:先搞清楚包里有什么再动手

做地形分析最怕的不是没有数据,而是数据来了不知道底细。我在一个山区公路选线项目里拿到一套全国地形地貌地质数据集,初看目录很完整,有DEM栅格、等高线矢量、地貌分类图层,还有一整文件夹的地质报告扫描件。真正用起来才发现,这些数据来自不同年代、不同坐标系、不同制图规范,直接丢进ArcGIS里叠图,地形和地貌能差出几百米。这份数据集的价值在于它把分散在全国各地的地形、地貌、地质图件和配套文档收拢到了一个资源包里,省去了满世界找资料的功夫,但代价是你要在动手前先把坐标系、图层关系和文档的可信度捋清楚。适合谁用呢,做区域地质调查、工程选址、自然资源普查的从业者,还有做课程设计的师生。接下来的内容,我会按“拆目录、统坐标、做分析、躲坑、验结果”的路子把这包数据讲透。

2. 资源包目录拆解:栅格、矢量、文档三类数据怎么配合使用

2.1 三种数据形态:DEM、图件图层、地质报告之间是什么关系

拿到资源包第一件事,不是急着加载图层,而是先把目录结构整理清楚。常见做法是把数据分成三个子目录。第一个是栅格数据,里面是数字高程模型,也就是DEM,分辨率常见的有30米和90米两档,沿海平原和西部山区都有覆盖。第二个是矢量数据,包含等高线、断层线、地层界线、地貌类型面状图层这些,通常是从地形图和地质图上矢量化出来的成果。第三个是文档资料,这个最容易被忽略,里面是区域地质调查报告、图幅说明书、地层柱状图、构造纲要图,格式以PDF和图片为主。

这三类数据的配合逻辑是一个递进关系。DEM是基础底图,提供地表起伏的连续表达;地貌类型图层是对地表形态的分类解释,比如喀斯特地貌、黄土地貌、流水地貌这些;地质图层进一步往下走,解释的是地下的物质组成和构造骨架。举个例子,你在DEM上看到一个陡峭的坡面,地貌图层会告诉你这是构造侵蚀剥蚀地貌还是滑坡堆积地貌,地质图层则告诉你这个位置的岩性是花岗岩还是页岩。三者叠起来才能回答“这个地形为什么是这样”以及“在这里搞建设风险大不大”这两个问题。单独拿任何一层出来,信息量都是不完整的。

2.2 文档资料的正确打开顺序:先读说明书的哪个部分

这套资源包里的大量文档资料,分量最重的当属区域地质调查报告。一份正规报告的结构是固定的:地理位置、区域地层、岩浆岩、构造、矿产地、结束语。我一般建议拿到报告先翻“地层”这一章,因为地层是后续所有分析的底座。报告里会附一张地层简表,按界、系、统、组四级把工作区内的岩层列出来,每一套地层旁边标注了岩性描述和分布范围。比如你看到“二叠系下统栖霞组(P1q):深灰色厚层灰岩”,立刻就能知道这片区域的地质背景是海相碳酸盐岩沉积,那么降雨条件下容易发育喀斯特现象,这就是一个从阅读到应用的直接推导。

图幅说明书的使用顺序也讲究。1:20万图幅说明书通常先是“概况”,讲交通位置、地形特征、水系发育情况,然后是“地质特征”,最后是“矿产”。注意一点,说明书里写的坐标基准是老图幅编号对应的经纬度范围,不是CGCS2000的经纬度,后面用的时候要小心。文档资料里还经常有航卫片解译记录和野外路线观察记录,这些原始记录的价值很高,里面记录了露头位置、产状测量数据、岩石标本编号,这些信息在矢量化图层里往往没有完整转出来。读文档时,我习惯顺手做一张表格,把每个地层代号、岩性关键词、产状数据、出露位置摘出来,后面做地层对比或者填图验证时直接查表,比翻几百页PDF快得多。

2.3 图例与符号系统:看不懂图例就谈不上用数据

地质图、地貌图、地形图各有各的图例系统,这套资源包里如果有标准图例文件那是运气好,没有的话就要靠自己去对照。地质图图例是色块加代号,用颜色区分时代,用代号表示岩石地层单位,比如用浅蓝色表示三叠系,用绿色表示侏罗系。地貌图的图例是形态符号系统,用不同线条和网点表示冲积平原、洪积扇、剥蚀台地等类型。地形图的图例最直观,主要是等高线、高程点、陡崖符号、植被符号。

在加载矢量图层之前,我强烈建议先单独打开图例文件或者说明书里的图例页,把每一类图斑对应什么含义弄清楚。踩过的坑是直接把地层界线图层和现代交通图层一起叠加,结果显示地层被道路切割得支离破碎,后来才发现两条线属性的语义完全不同,一条是地质界线,一条是道路中线,本来就该分开管理。正确的做法是把图例信息整理成一个属性对照表,在GIS里给每个图层配上独立的符号系统,地质图层用地质图规范配色,地貌图层用地貌分类配色,DEM做灰度渐变或山体阴影,这样视觉上信息才不会打架。

3. 坐标系与投影统一:数据叠不合的根源与批量转换方案

3.1 新老坐标基准混用:先查每一层的坐标系再开始干活

从这套资源包里随便挑几个图层查看属性,坐标系大概率是几种基准混着的。老的数据用的是北京54坐标系,中间年代的数据用西安80坐标系,新一些的用CGCS2000,部分从公开网站下载的DEM是WGS84地理坐标系。直接叠图,平面偏移量在北京54和CGCS2000之间可以达到百米量级,在西部地区甚至更大。所以第一步是在ArcGIS或QGIS里逐个图层查看属性表的坐标系信息,我一般是写一个小循环把目录下所有矢量文件的空间参考信息打出来,快速摸清底数。

在QGIS里有个方便的办法,用Python控制台可以批量打印图层的CRS信息,代码如下:

from qgis.core import QgsVectorLayer, QgsRasterLayer import os data_dir = r"D:\geodata\vector" for fname in os.listdir(data_dir): if fname.endswith(".shp"): layer = QgsVectorLayer(os.path.join(data_dir, fname), fname, "ogr") if layer.isValid(): crs = layer.crs() print(fname, " -> ", crs.authid(), crs.description())

这段代码遍历数据目录里的全部shapefile文件,逐个加载并输出坐标系编号和描述。打印结果里如果同时出现EPSG:4214(北京54地理坐标)、EPSG:4610(西安80地理坐标)、EPSG:4490(CGCS2000地理坐标)这些编号,说明数据源混用了多个基准,后续必须统一转换才能做叠加分析。注意,这里的EPSG编号是地理坐标系的,如果数据源里存储的是投影坐标系,输出会变成EPSG:2437这些高斯-克吕格投影编号,对应的中央经线不同也要处理。

3.2 矢量数据批量转换:用GDAL做坐标基准统一

确定好每层数据的坐标基准之后,接下来要做的是统一转换。我的习惯是全部转到CGCS2000坐标系,这也是当前国家标准要求的数据基准。ArcGIS里有“Project”工具能做转换,但如果图层数量多,用GDAL命令行批量操作效率更高。

gdalwarp -s_srs EPSG:4610 -t_srs EPSG:4490 -r bilinear input_dem.tif output_dem_cgcs2000.tif ogr2ogr -s_srs EPSG:4214 -t_srs EPSG:4490 -f "ESRI Shapefile" input_shp_54.shp output_shp_2000.shp

第一条命令处理栅格数据,-s_srs指定输入坐标系为西安80地理坐标,-t_srs指定输出为CGCS2000地理坐标,-r bilinear是重采样方法,用双线性内插可以保持地形表面的光滑性,如果你是做坡度分析,双线性比最近邻要好。第二条命令处理矢量数据,把北京54坐标系的数据转换到CGCS2000,-f "ESRI Shapefile"指定输出格式。需要注意一点,GDAL的转换是七参数还是三参数取决于你的环境,默认情况下用的是一个简化模型,在局部地区可能有两到三米的残余误差。如果你在工程选址这种需要厘米级精度的场景,建议向测绘部门索取所在区域的精确转换参数,在GDAL命令里通过-ct参数指定一个完整的坐标转换管线。

3.3 投影分带问题:跨带数据拼接时的计算规则

除了基准不同,投影分带也是一个容易翻车的点。高斯-克吕格投影分为3度带和6度带两种分带方式,不同图幅可能落在不同的投影带内。资源包里的1:20万地质图按经纬度分幅,每幅图跨越的经度范围可能是2度或3度,相邻两幅图可能出现在不同投影带。直接合并成一个图层后,在投影带边界处会出现明显的错位或断裂,看起来像数据出错了,其实是投影带切换导致。

我一般用两种方案应对。第一种是做分析时侧重一个图幅范围,不跨带;第二种是统一转成Albers等积投影或Lambert等角投影,这类投影适合全国或大区域分析。做全国尺度的地貌统计时,我习惯把矢量图层统一转成Albers投影,因为Albers是等积投影,面积量算不扭曲,统计各个地貌类型的面积占比时才靠谱。投影转换命令和前面格式类似,只是目标EPSG编号换成Albers投影的编号,国内常用的是EPSG:102008或自定义中央经线参数。转换后用ArcGIS的“Project”工具或QGIS的“Reproject Layer”再做一次验证,确保没有产生飞点,也就是远离源数据区域的不合理坐标点。

4. 坡度、坡向与地貌分析的完整流程:参数与每一步验证

4.1 从DEM到坡度图:先填洼地再算坡度,顺序不能反

把数据坐标统一之后,开始正式做地形分析。最常见的第一步是计算坡度。从我拿到这套数据集的经验来看,DEM数据往往是经过初步处理的,但局部区域仍然可能存在数据空洞和洼地。计算坡度之前如果不填洼地,洼地边缘会产生虚假的陡坡,后续的坡度和坡向统计就会失真。

在ArcGIS里填洼地的工具叫“Fill”,在QGIS里对应的是“Fill sinks (wang & liu)”这个算法。填洼的阈值参数一般设成DEM分辨率的2到3倍。30米分辨率的DEM,阈值设在60到90米之间比较稳妥。阈值设大了会过度平缓地形,把真实的洼地也填掉;设小了又填不干净,结果里还有零零散散的小坑。

填完洼地之后开始计算坡度,这里有一个关键选择:用度(degree)还是用百分比(percent)。我一般推荐用度,因为后续做坡度分级时阈值更直观。ArcGIS的“Slope”工具里把输出测量单位选成DEGREE,QGIS的“Slope”工具则在参数里选择“Degrees”。

import rasterio import numpy as np with rasterio.open(r"D:\geodata\dem_filled.tif") as src: dem = src.read(1) transform = src.transform # 计算x和y方向的梯度 x_grad, y_grad = np.gradient(dem, transform[0], -transform[4]) # 坡度(度) slope = np.degrees(np.arctan(np.sqrt(x_grad**2 + y_grad**2))) slope = np.where(slope < 0, 0, slope)

这段Python代码用numpy的gradient函数计算DEM的梯度,然后通过反正切公式把梯度转换成坡度值。transform[0]是像素宽度,transform[4]是像素高度,注意这个值是负数,所以要用负号取绝对值。代码最后把计算出的负值强制归零,防止数值误差产生非法角度。这个脚本适合在没有ArcGIS环境的时候快速算坡度,算出来的结果和GIS工具基本一致,但要注意边界处的梯度值因为没有邻居像素,会略偏低,分析时可以不统计边界一圈。

4.2 坡度分级:不同行业标准下的分级阈值怎么选

有了连续的坡度栅格,接下来要做的通常是分级。不同应用场景的分级标准差异很大,这也是一个容易照搬经验翻车的地方。我做工程选址时用的是《土地利用现状分类》的坡度分级标准,把坡度分成小于2度、2到6度、6到15度、15到25度、大于25度五档。做地质灾害易发性分析时,分级就粗糙一些,小于10度、10到30度、大于30度三档就够了,因为地质灾害的坡度敏感区间是宽泛的,分太细反而得不到清晰规律。

在工具操作层面,ArcGIS的“Reclassify”工具可以完成重新分级,QGIS里对应“Raster calculator”做条件赋值。我更习惯在Python里直接完成,因为可以同时输出各级别的面积占比,一步到位:

import rasterio import numpy as np with rasterio.open(r"D:\geodata\slope_deg.tif") as src: slope = src.read(1) profile = src.profile classes = np.zeros_like(slope, dtype=np.uint8) classes[(slope >= 0) & (slope < 2)] = 1 classes[(slope >= 2) & (slope < 6)] = 2 classes[(slope >= 6) & (slope < 15)] = 3 classes[(slope >= 15) & (slope < 25)] = 4 classes[slope >= 25] = 5 # 统计各分级面积占比 total = classes.size for i in range(1, 6): pct = (classes == i).sum() / total * 100 print(f"分级{i}: {pct:.2f}%")

代码逻辑很简单,先生成一个和坡度栅格一样大小的空数组,然后按阈值条件给每个像元赋分级代号,最后统计每个级别的占比并打印。重点说一下为什么用uint8类型存储分级结果,因为分级号只有1到5,用8位无符号整形就够了,文件体积比浮点型小四分之三,后续处理速度也快。如果你要输出的分级图用于制图,记得在输出TIFF文件时把profile里的dtype改成uint8,再把nodata设成0,这样在GIS里显示不会有黑边。

4.3 地貌类型图层与坡度叠加:寻找高坡度和脆弱地貌的耦合区域

坡度图算出来只是第一步,真正的价值在于和地貌类型图层叠加分析。我的一个常用操作是把坡度分级和地貌类型做交叉统计,找出“高坡度+易滑地貌”的组合区域,这类区域往往是地质灾害的高发区。具体做法是对地貌图层做“Zonal histogram”操作,以地貌类型为分区对象,统计每个地貌分区内坡度分级的占比。

QGIS里可以这样操作:先用地貌面图层把坡度栅格裁剪出来,然后通过“Zonal statistics”按地貌类型区域统计平均坡度和最大坡度。ArcGIS里对应的工具是“Zonal Statistics as Table”。这个操作得到的表格是后续风险评估的直接输入数据。注意一个细节,做交叉统计前一定要确认地貌图层的几何没有自相交或重叠,否则统计结果会重复计算。全部检查完成后,把统计表格导成CSV,用透视表拉出各地貌类型的坡度分布直方图,马上就能看出哪些地貌在陡坡区间占比高。

5. 避坑指南:坐标混用、接边断裂与OCR误识别实录

5.1 西安80和CGCS2000混用导致叠图偏移几百米

现象是地形图和地质图叠加之后,等高线和地层界线之间出现一个稳定的平面偏移,偏移量在城区的建筑物上特别明显,同一栋楼在两张图上位置差了大约三四百米,放到大比例尺图上一眼就能发现。原因是我前期没有逐个图层检查坐标系,地形图是WGS84导出的数据,地质图是西安80坐标系,两种基准在同一地区的椭球面差异被投影放大后形成了这个偏移。解决的办法是把所有图层重新投影到统一的CGCS2000坐标系,对矢量数据用ogr2ogr做转换,对栅格数据用gdalwarp做转换,转换后用一幅图的明显地物点做精度校验,比如河流交叉点或独立山头的高程点,确认误差在10米以内才继续后续分析。

5.2 跨图幅拼接接边处要素断裂

现象是拼接两幅相邻地质图后,原本应该连续的地层界线在接边处出现错开或断开,有的地层多边形边界在接缝处差了几十米到上百米,搭不起来。原因有两方面,一方面是两幅图矢量化时的精度不一致,老图幅是手工数字化,误差大;另一方面是投影分带不同,接边区域被强行变换到同一个投影下,产生了形变。我的解决方法是先检查接边处各要素的坐标残差,如果残差在50米以内,用ArcGIS的“Snap”工具做微调,把断点吸附到相邻图幅的对应端点上;如果残差超过200米,多半是坐标系没统一对,回到第3章的转换流程把投影重做一遍。处理完后用“Integrate”工具整合公共边,保证多边形拓扑闭合,再做一次拓扑检查,确保没有悬挂弧段。

5.3 扫描件OCR后地层代号误识别

现象是文档资料里的区域地质调查报告扫描件经过OCR识别之后,地层代号出现大量错误,最典型的是把“C2”识别成“C2”后面的字母被丢掉,把“P1q”识别成“P1g”,“灰岩”被识别成“灰告”,甚至“大理岩”被识别成“大理石”。原因很直接,老式印刷体在扫描后对比度差,OCR引擎对地质代号里的数字和上下标支持不好,识别置信度低的字段被强行输出。我的处理办法是不直接信任OCR结果,而是把识别出的文本对照图幅说明书中的地层简表做人工校对,重点检查地层代号和岩性关键词。我一般把地层代号列表提取出来,然后和GB标准的地层代号规范表做模糊匹配,匹配不上的手动改正。这是一个费时但绕不过去的环节,因为地层代号如果错了,后续所有基于地层信息的分析结果都会跟着错。

5.4 断层线是示意性的,不是精确定位边界

现象是拿着地质图上的断层线直接去现场验证,结果在图上标注的位置没有找到断裂迹象,而在图面外偏移几十米的地方发现了断层破碎带。原因是在中比例尺区域地质图上,断层线表达的是断层的平面投影位置,精度受制于图面比例尺和原始调查精度,不能等同于实测剖面定位。解决的办法是把断层数据用作分析背景而不是定位依据,做工程选址时把断层两侧按规范要求设置避让距离,一般要求是断层两侧各避让100到500米,具体取决于工程等级。我处理断层图层时,会先在属性表里加一个“精度等级”字段,把标注为实测的断层和推测断层分开管理,推测断层用虚线表达,分析时降低权重。

5.5 文字报告里提取的坐标点漏带号

现象是文档资料里的坐标点导入GIS后全部跑到了视野之外,或者出现在非洲或海上的位置。原因是报告文本里记录的坐标是高斯坐标但没写带号,比如Y坐标是“38512345”,前两位“38”是带号,后面是500公里加偏移后的数值,导入时如果把整个数值当成了普通横坐标,位置肯定不对。解决方法是写配置文件时强制把带号从Y坐标里拆出来,X坐标保持完整,然后在GIS里设置好正确的中央经线。我常用的一段Python代码来处理坐标点:

import pandas as pd df = pd.read_csv(r"D:\geodata\points_raw.csv") # Y坐标格式:带号(2位) + 500km偏移后的值 df["band"] = df["y_raw"] // 1_000_000 df["y_corrected"] = df["y_raw"] % 1_000_000 - 500_000 # 中央经线 = 带号 * 3(3度带)或带号 * 6 - 3(6度带) df["central_meridian"] = df["band"] * 3 df.to_csv(r"D:\geodata\points_fixed.csv", index=False)

这段代码把Y坐标拆成带号和去掉500公里假偏移后的真实横坐标,再根据带号推算中央经线。逻辑上需要注意,带号的计算前提是坐标数据是高斯投影坐标,如果数据源是经纬度,这个脚本就不适用。处理后的点数据再结合带中央经线做投影定义,加载到GIS里就能落到正确位置。

6. 数据自查与出图交付:几个能救命的小技巧

在把分析结果交付出去之前,我习惯做三遍快速自查,每一遍都有对应的工具和技巧。第一遍查DEM,方法是在GIS里做山体阴影渲染,把光照方向设在西北方向,如果发现DEM表面有横向条带或者规则的棋盘格状纹理,说明原始DEM有采集噪声,需要用低通滤波做平滑。第二遍等高线叠加DEM,检查等高线高程与DEM栅格值是否一致,如果同一位置等高线标着500米而DEM栅格读出480多米,说明等高线的基准面和DEM不一致,这个数据不能用来做精细剖面分析。第三遍查地貌图斑的拓扑关系,用“Check geometries”工具检查面图层是否有重叠或缝隙,有缝隙的话栅格化后统计面积会漏。

出图的时候有一个技巧我想多说一句。地貌类型图直接用单色填充出图,颜色层次感差,不直观。我的做法是把DEM山体阴影栅格设置为图层背景,透明度和灰度调整到位,再把地貌面图层叠加上去,面的填充透明度设为30%左右,这样既能看清地形起伏,又能识别地貌类型边界。ArcGIS里在图层属性的Display选项卡里调透明度,QGIS里在图层属性的Transparency选项卡里调。最终出图时,图例里必须同时包含地形、地貌、地质三套符号的说明,否则看图的人分不清线划的语义。

还有一个小习惯,在那次被坐标系坑过之后,我每次拿到新的地形地质数据集,第一件事就是在QGIS里跑一遍CRS遍历脚本,把所有图层的坐标系信息导出成一张表格,随项目文档一起归档。后续任何人接手这个项目,先看这张表,立刻就能知道哪些图层可以直接叠加,哪些需要先转换。坐标转换做完之后再花十分钟做精度验证,验证方法是在两幅图上各取三个明显同名地物点,量算转换后的残差。从那以后,这套流程我再也没跳过,已经成了肌肉记忆。希望帮到你,能少走这些弯路。

本文还有配套的精品资源,点击获取

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询