拓冰网站定制 
简介面向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 rD:\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:4490CGCS2000地理坐标这些编号说明数据源混用了多个基准后续必须统一转换才能做叠加分析。注意这里的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”工具里把输出测量单位选成DEGREEQGIS的“Slope”工具则在参数里选择“Degrees”。import rasterio import numpy as np with rasterio.open(rD:\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(rD:\geodata\slope_deg.tif) as src: slope src.read(1) profile src.profile classes np.zeros_like(slope, dtypenp.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(rD:\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 # 中央经线 带号 * 33度带或带号 * 6 - 36度带 df[central_meridian] df[band] * 3 df.to_csv(rD:\geodata\points_fixed.csv, indexFalse)这段代码把Y坐标拆成带号和去掉500公里假偏移后的真实横坐标再根据带号推算中央经线。逻辑上需要注意带号的计算前提是坐标数据是高斯投影坐标如果数据源是经纬度这个脚本就不适用。处理后的点数据再结合带中央经线做投影定义加载到GIS里就能落到正确位置。6. 数据自查与出图交付几个能救命的小技巧在把分析结果交付出去之前我习惯做三遍快速自查每一遍都有对应的工具和技巧。第一遍查DEM方法是在GIS里做山体阴影渲染把光照方向设在西北方向如果发现DEM表面有横向条带或者规则的棋盘格状纹理说明原始DEM有采集噪声需要用低通滤波做平滑。第二遍等高线叠加DEM检查等高线高程与DEM栅格值是否一致如果同一位置等高线标着500米而DEM栅格读出480多米说明等高线的基准面和DEM不一致这个数据不能用来做精细剖面分析。第三遍查地貌图斑的拓扑关系用“Check geometries”工具检查面图层是否有重叠或缝隙有缝隙的话栅格化后统计面积会漏。出图的时候有一个技巧我想多说一句。地貌类型图直接用单色填充出图颜色层次感差不直观。我的做法是把DEM山体阴影栅格设置为图层背景透明度和灰度调整到位再把地貌面图层叠加上去面的填充透明度设为30%左右这样既能看清地形起伏又能识别地貌类型边界。ArcGIS里在图层属性的Display选项卡里调透明度QGIS里在图层属性的Transparency选项卡里调。最终出图时图例里必须同时包含地形、地貌、地质三套符号的说明否则看图的人分不清线划的语义。还有一个小习惯在那次被坐标系坑过之后我每次拿到新的地形地质数据集第一件事就是在QGIS里跑一遍CRS遍历脚本把所有图层的坐标系信息导出成一张表格随项目文档一起归档。后续任何人接手这个项目先看这张表立刻就能知道哪些图层可以直接叠加哪些需要先转换。坐标转换做完之后再花十分钟做精度验证验证方法是在两幅图上各取三个明显同名地物点量算转换后的残差。从那以后这套流程我再也没跳过已经成了肌肉记忆。希望帮到你能少走这些弯路。本文还有配套的精品资源点击获取 全国地形地貌地质数据集处理：从DEM到坡度分析的完整实践
简介面向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 rD:\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:4490CGCS2000地理坐标这些编号说明数据源混用了多个基准后续必须统一转换才能做叠加分析。注意这里的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”工具里把输出测量单位选成DEGREEQGIS的“Slope”工具则在参数里选择“Degrees”。import rasterio import numpy as np with rasterio.open(rD:\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(rD:\geodata\slope_deg.tif) as src: slope src.read(1) profile src.profile classes np.zeros_like(slope, dtypenp.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(rD:\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 # 中央经线 带号 * 33度带或带号 * 6 - 36度带 df[central_meridian] df[band] * 3 df.to_csv(rD:\geodata\points_fixed.csv, indexFalse)这段代码把Y坐标拆成带号和去掉500公里假偏移后的真实横坐标再根据带号推算中央经线。逻辑上需要注意带号的计算前提是坐标数据是高斯投影坐标如果数据源是经纬度这个脚本就不适用。处理后的点数据再结合带中央经线做投影定义加载到GIS里就能落到正确位置。6. 数据自查与出图交付几个能救命的小技巧在把分析结果交付出去之前我习惯做三遍快速自查每一遍都有对应的工具和技巧。第一遍查DEM方法是在GIS里做山体阴影渲染把光照方向设在西北方向如果发现DEM表面有横向条带或者规则的棋盘格状纹理说明原始DEM有采集噪声需要用低通滤波做平滑。第二遍等高线叠加DEM检查等高线高程与DEM栅格值是否一致如果同一位置等高线标着500米而DEM栅格读出480多米说明等高线的基准面和DEM不一致这个数据不能用来做精细剖面分析。第三遍查地貌图斑的拓扑关系用“Check geometries”工具检查面图层是否有重叠或缝隙有缝隙的话栅格化后统计面积会漏。出图的时候有一个技巧我想多说一句。地貌类型图直接用单色填充出图颜色层次感差不直观。我的做法是把DEM山体阴影栅格设置为图层背景透明度和灰度调整到位再把地貌面图层叠加上去面的填充透明度设为30%左右这样既能看清地形起伏又能识别地貌类型边界。ArcGIS里在图层属性的Display选项卡里调透明度QGIS里在图层属性的Transparency选项卡里调。最终出图时图例里必须同时包含地形、地貌、地质三套符号的说明否则看图的人分不清线划的语义。还有一个小习惯在那次被坐标系坑过之后我每次拿到新的地形地质数据集第一件事就是在QGIS里跑一遍CRS遍历脚本把所有图层的坐标系信息导出成一张表格随项目文档一起归档。后续任何人接手这个项目先看这张表立刻就能知道哪些图层可以直接叠加哪些需要先转换。坐标转换做完之后再花十分钟做精度验证验证方法是在两幅图上各取三个明显同名地物点量算转换后的残差。从那以后这套流程我再也没跳过已经成了肌肉记忆。希望帮到你能少走这些弯路。本文还有配套的精品资源点击获取