简介这份中国喀斯特岩溶空间分布矢量数据集面向地理信息、地质地貌与环境规划领域的研究者与从业者用于分析岩溶地块边界、岩性类型及空间分布规律。资源包共8个文件约1.2MB以SHP矢量数据为核心配套SHX、SBX、SBN空间索引文件DBF属性表PRJ空间参考CPG编码说明及XML元数据构成一套可直接在GIS软件中加载使用的完整数据。属性字段涵盖岩性分类连续与不连续碳酸盐岩、地块面积与周长、岩性文本标签便于评估岩溶发育程度及地表地下水系统规模。目前已有394人学习下载适合用于喀斯特地貌研究、农业与水利规划、旅游开发等场景也可作为GIS空间分析与制图练习的底图数据。1. 喀斯特岩溶空间分布矢量数据集从一张 SHP 到一套可复用的分析底图拿到“中国Karsts喀斯特岩溶空间分布矢量数据集SHP数据”这个标题很多人第一反应是去找下载链接但真正卡住一线工程师的往往不是下载而是拿到 shp 之后怎么用、坐标系对不对、属性表能不能筛、能不能和县域行政区划边界 shp 做叠加。喀斯特岩溶空间分布本质上是把碳酸盐岩出露区、岩溶地貌类型、溶蚀强度等地质信息落到多边形上形成可查询、可裁剪、可统计的矢量图层。它适合做水文地质调查、生态敏感性评价、工程选址避让、碳汇估算前期底图的人。这篇笔记按“数据是什么→怎么读→怎么裁→怎么避坑→怎么进阶”的顺序把一套 SHP 从入库到出图的完整链路讲清楚中间会穿插 dwg转shp、shp转txt、渔网分割shp 这些热搜里高频出现的操作场景。2. 先搞懂喀斯特岩溶 SHP 里到底存了什么2.1 岩溶空间分布的三种几何表达与选型理由喀斯特岩溶空间分布矢量数据集通常不是单一图层而是按比例尺和专题拆成若干 SHP 文件。常见几何类型有三类面状多边形表达碳酸盐岩连片出露区或岩溶地貌分区线状要素表达地下河、溶洞通道、断裂带点状要素表达落水洞、天窗、泉点。做区域统计和叠加分析时面状图层是主力做线性工程避让时线状和点状图层才关键。选型上如果你要做全省岩溶面积占比优先用面状图层因为面积计算直接基于多边形。如果你要做某条铁路的岩溶风险分段线状和点状图层叠加缓冲区更实用。很多公开数据集只给面状线状和点状需要从地质图或水文地质图里单独提取。常见做法是先用面状图层圈定“大范围”再用高精度线点图层做“局部修正”。属性表是第二个要关注的点。一个合格的岩溶 SHP 至少应包含岩性代码、地貌类型、出露面积、数据来源、比例尺这几列。岩性代码常见的有石灰岩、白云岩、白云质灰岩等分类地貌类型可能分峰丛、峰林、溶丘、溶原。没有属性表的 SHP 等于一堆无意义的边界后续没法按类型筛选。提示拿到 SHP 先别急着打开符号系统用ogrinfo或 ArcGIS 的“要素类属性”看一眼字段名和字段类型比在图上瞎点效率高得多。2.2 用 GDAL 在本地读一遍 SHP 的最小命令不管后续用 ArcGIS、QGIS 还是 PostGIS第一步都建议用 GDAL 做一次“体检”。下面这段 bash 命令能一次性输出图层名、几何类型、坐标系、字段定义和要素数量。# 查看 SHP 基本信息不加载图形界面 ogrinfo -so -al china_karst_distribution.shp # 如果目录下有多个 SHP批量看坐标系和要素数 for f in *.shp; do echo $f ogrinfo -so -al $f | grep -E Layer name|Geometry|Feature Count|AUTHORITY done-so表示只输出摘要不遍历要素-al表示所有图层。输出里重点看AUTHORITY[EPSG,xxxx]这决定后续能不能直接和县域行政区划边界 shp 叠加。如果看到EPSG:4326说明是地理坐标系面积计算前必须投影如果看到EPSG:4490或EPSG:4610也是地理坐标系同样要投影。只有看到EPSG:3857、EPSG:4527这类投影坐标系才能直接算面积。要素数量也要留意。一个全国尺度的岩溶面状图层要素数通常在几千到几万之间。如果只有几十个要素大概率是高度概化后的分区图不适合做精细分析。如果超过十万可能是把碎小多边形都保留了需要先做融合。2.3 属性表字段的清洗与标准化原始 SHP 的属性表经常有字段名超长、中文乱码、类型不一致的问题。用 Python 的geopandas可以快速做一轮清洗。import geopandas as gpd gdf gpd.read_file(china_karst_distribution.shp) # 查看字段名和类型 print(gdf.dtypes) print(gdf.columns.tolist()) # 重命名关键字段避免中文乱码导致后续脚本报错 rename_map { 岩性代码: litho_code, 地貌类型: geom_type, 出露面积: area_km2, 数据来源: source } gdf gdf.rename(columns{k: v for k, v in rename_map.items() if k in gdf.columns}) # 统一岩性代码大小写去除首尾空格 if litho_code in gdf.columns: gdf[litho_code] gdf[litho_code].astype(str).str.strip().str.upper() # 检查几何有效性修复自相交 invalid gdf[~gdf.is_valid] print(f无效几何数量: {len(invalid)}) if len(invalid) 0: gdf[geometry] gdf.buffer(0) gdf.to_file(china_karst_clean.shp, encodingUTF-8)这段代码做了四件事重命名字段、统一代码格式、检查几何有效性、用buffer(0)修复自相交。buffer(0)是处理自相交多边形的常用技巧但要注意它可能轻微改变边界对面积精度要求极高的场景要谨慎。encodingUTF-8是为了避免中文属性在跨平台时乱码如果目标平台只认 GBK改成encodingGBK。字段清洗后建议把常用筛选条件写成独立列。比如把“石灰岩”和“白云岩”归为“碳酸盐岩”把“峰丛”“峰林”归为“正向岩溶地貌”这样后续做统计时不用反复写复杂条件。3. 把岩溶 SHP 和行政区划、DEM 叠起来用3.1 与县域行政区划边界 SHP 做相交统计岩溶分布数据单独看意义有限一旦和县域行政区划边界 shp 叠加就能算出每个县的岩溶面积占比。这是生态评估和工程规划里最常用的操作。import geopandas as gpd import pandas as pd karst gpd.read_file(china_karst_clean.shp) counties gpd.read_file(county_boundary.shp) # 统一投影到适合全国面积计算的 Albers 投影 albers EPSG:4527 # 中国 Albers 等面积投影 karst_albers karst.to_crs(albers) counties_albers counties.to_crs(albers) # 相交运算保留县域属性 intersect gpd.overlay(karst_albers, counties_albers, howintersection) # 计算每个县的岩溶面积 intersect[karst_area_km2] intersect.geometry.area / 1e6 counties_albers[county_area_km2] counties_albers.geometry.area / 1e6 # 按县汇总 summary intersect.groupby(county_name)[karst_area_km2].sum().reset_index() summary summary.merge( counties_albers[[county_name, county_area_km2]], oncounty_name, howleft ) summary[karst_ratio] summary[karst_area_km2] / summary[county_area_km2] summary.to_csv(karst_by_county.csv, indexFalse, encodingutf-8-sig)gpd.overlay的howintersection会保留两个图层的所有属性适合做统计。如果县域边界有重叠或缝隙先做一次dissolve或拓扑检查。EPSG:4527是中国 Albers 等面积投影适合全国尺度面积计算如果只做某个省可以用该省对应的 UTM 带号精度更高。karst_ratio这一列是核心产出。通常认为占比超过 30% 的县属于岩溶发育重点区超过 60% 属于强烈发育区。这个阈值不是固定的要根据具体项目的地质背景调整。3.2 用 DEM 提取坡度后与岩溶图层做分区统计岩溶区的地形坡度直接影响溶蚀强度和工程难度。把 DEM 提取的坡度栅格和岩溶面状图层叠加可以算出每个岩溶多边形内的平均坡度。import rasterio from rasterio.mask import mask import numpy as np import geopandas as gpd # 读取 DEM 并计算坡度 with rasterio.open(dem.tif) as src: dem src.read(1) transform src.transform crs src.crs # 简单坡度计算实际项目建议用 richdem 或 GDAL 的 gdaldem dy, dx np.gradient(dem, transform[4], transform[0]) slope np.degrees(np.arctan(np.sqrt(dx**2 dy**2))) # 将坡度写回栅格 with rasterio.open( slope.tif, w, driverGTiff, heightslope.shape[0], widthslope.shape[1], count1, dtypeslope.dtype, crscrs, transformtransform ) as dst: dst.write(slope, 1) # 按岩溶多边形做分区统计 karst gpd.read_file(china_karst_clean.shp).to_crs(crs) with rasterio.open(slope.tif) as src: for idx, row in karst.iterrows(): try: out_image, _ mask(src, [row.geometry], cropTrue, nodatanp.nan) mean_slope np.nanmean(out_image) karst.loc[idx, mean_slope] mean_slope except Exception: karst.loc[idx, mean_slope] np.nan karst.to_file(karst_with_slope.shp, encodingUTF-8)np.gradient算坡度是简化版实际项目建议用gdaldem slope或richdem它们处理了投影单位和边缘效应。mask函数按多边形裁剪栅格cropTrue减少内存占用。如果多边形数量多逐行循环会很慢可以改用rasterstats的zonal_stats底层用rasterio和shapely做了优化。得到mean_slope后可以按坡度分级小于 8 度适合耕作和建设8 到 25 度适合林业大于 25 度属于生态脆弱区。岩溶区如果同时满足“碳酸盐岩出露”和“坡度大于 25 度”就是石漠化敏感区。3.3 渔网分割 SHP 做格网化统计当需要把岩溶分布做成规则格网参与模型计算时渔网分割 shp 是常用手段。下面用geopandas生成 10km×10km 渔网再和岩溶图层相交。import geopandas as gpd from shapely.geometry import box import numpy as np karst gpd.read_file(china_karst_clean.shp).to_crs(EPSG:4527) minx, miny, maxx, maxy karst.total_bounds # 生成 10km 渔网 cell_size 10000 # 单位米 cols list(np.arange(minx, maxx, cell_size)) rows list(np.arange(miny, maxy, cell_size)) cells [box(x, y, x cell_size, y cell_size) for x in cols for y in rows] grid gpd.GeoDataFrame({geometry: cells}, crsEPSG:4527) # 相交并统计每个格网的岩溶面积 intersect gpd.overlay(grid, karst, howintersection) intersect[karst_area_km2] intersect.geometry.area / 1e6 grid_summary intersect.groupby(intersect.index)[karst_area_km2].sum().reset_index() grid grid.reset_index().merge(grid_summary, onindex, howleft).fillna(0) grid.to_file(karst_grid_10km.shp, encodingUTF-8)cell_size根据研究尺度调整全国用 10km 到 50km省级用 1km 到 5km县级用 100m 到 500m。渔网生成后overlay的intersection会保留每个格网和每个岩溶多边形的交集groupby按格网索引汇总。注意grid.index在overlay后可能变化所以先reset_index再合并。格网化后的数据可以直接导入 MaxEnt、InVEST 等模型做空间预测。如果格网内岩溶面积为 0fillna(0)会补零避免模型把缺失值当成有效值。4. 喀斯特 SHP 处理里最容易翻车的五个坑4.1 坐标系没投影就直接算面积现象用gdf.area算出来的面积单位是“平方度”数值小得离谱或者和实际面积差几个数量级。原因SHP 是 EPSG:4326 地理坐标系经纬度直接算面积没有物理意义。解决先to_crs到等面积投影全国用 EPSG:4527省级用对应 UTM 带号县级用高斯克吕格投影。判断方法gdf.crs.is_geographic返回 True 就必须投影。4.2 属性表中文乱码导致筛选失效现象在 QGIS 里看到属性表中文正常但用 Python 读取后字段名变成\u5ca9\u6027筛选条件写不对。原因SHP 的 DBF 文件默认编码可能是 GBK而geopandas默认按 UTF-8 读。解决读取时指定encodingGBK或者先用ogrinfo确认编码再用gpd.read_file(..., encodingGBK)。写入时统一用encodingUTF-8并在文档里注明。4.3 几何自相交导致叠加分析报错现象gpd.overlay报TopologyException或者相交结果里出现面积为零的碎片。原因原始 SHP 的多边形存在自相交、悬挂边或重复节点。解决先gdf.is_valid检查对无效几何用buffer(0)修复。如果buffer(0)后仍有问题用shapely.validation.make_valid做更精细的修复。修复后重新计算面积和修复前对比差异超过 1% 的要人工检查。4.4 渔网分割后格网索引错位现象overlay之后按索引汇总发现某些格网的岩溶面积明显偏大或偏小。原因grid在overlay前没有reset_index导致groupby时索引对不上。解决生成渔网后立即grid grid.reset_index()把原始索引存成一列overlay后用这一列做groupby。另外overlay会丢弃没有交集的格网汇总后要merge回完整渔网并fillna(0)。4.5 把不同比例尺的岩溶图层混用现象把 1:50 万和 1:100 万的岩溶图层合并后边界对不上出现大量狭长碎片。原因不同比例尺的概化程度不同同一地理边界在两张图上的位置有偏差。解决优先使用同一比例尺的数据。如果必须混用先做snap或buffer容差处理容差取数据精度的 1 到 2 倍。更稳妥的做法是以高精度图层为准低精度图层只做属性补充不做几何合并。5. 从 SHP 到 WKT、TXT 和 3D Tiles 的进阶转换5.1 把岩溶多边形导出为 WKT 和 TXT有些平台只接受文本格式的几何比如某些在线空间分析接口或数据库导入工具。把 SHP 转成 WKT 或 TXT 是常见需求。import geopandas as gpd gdf gpd.read_file(china_karst_clean.shp) # 导出为 WKT保留关键属性 gdf[wkt] gdf.geometry.apply(lambda geom: geom.wkt) gdf[[litho_code, geom_type, wkt]].to_csv( karst_wkt.txt, sep\t, indexFalse, encodingutf-8 ) # 如果只需要坐标对导出为 TXT with open(karst_coords.txt, w, encodingutf-8) as f: for idx, row in gdf.iterrows(): coords list(row.geometry.exterior.coords) f.write(ffeature_{idx}\t{row.get(litho_code, NA)}\n) for x, y in coords: f.write(f{x}\t{y}\n)WKT 适合直接粘贴到 PostGIS 的ST_GeomFromText函数里。TXT 坐标对适合导入到没有 GIS 库的轻量级绘图工具。注意exterior.coords只取外环如果多边形有内环孔洞需要额外处理interiors。5.2 用 GDAL 把 SHP 转成 3D Tiles 的预处理shp转3dtiles 是热搜里高频出现的需求但 SHP 是二维矢量转 3D Tiles 前需要先有高度信息。常见做法是给岩溶多边形赋一个高程属性再拉伸成体块。# 第一步给 SHP 添加高程字段可以用 DEM 采样或统一赋值 ogr2ogr -f ESRI Shapefile karst_3d.shp china_karst_clean.shp \ -sql SELECT *, 1000 AS elevation FROM china_karst_clean # 第二步用 GDAL 的 gdal_rasterize 或第三方工具做拉伸 # 实际项目中常用 Cesium ion 或 py3dtiles 做矢量拉伸-sql里的1000 AS elevation是统一赋值实际项目应该用gdalwarp或rasterstats从 DEM 采样每个多边形的平均高程。得到带高程的 SHP 后用py3dtiles或Cesium ion的矢量拉伸功能生成 3D Tiles。注意 3D Tiles 的坐标系通常是 EPSG:4978地心直角坐标系转换时要做基准面变换。5.3 用 QGIS 做快速可视化和出图不是所有场景都需要写代码。QGIS 做岩溶分布出图很快关键是符号化和标注。操作路径参数建议按岩性分类着色图层属性 → 符号化 → 分类字段选litho_code色带选地形色带添加县域边界图层 → 添加图层 → 添加矢量图层设置边界为灰色细线不填充标注岩溶类型图层属性 → 标注字段选geom_type字体大小 8避让开启输出地图项目 → 新建打印布局比例尺 1:100 万添加指北针和图例QGIS 的“分类”符号化会自动读取字段唯一值适合岩性代码这种离散字段。如果字段值太多先用field calculator做分组把低频值归为“其他”。出图时注意图例顺序和颜色对比度岩溶区常用暖色调非岩溶区用冷色调。5.4 一个我常犯的错误忽略数据精度和适用尺度早期做省级岩溶统计时我直接拿全国 1:100 万的岩溶 SHP 去算某个县的岩溶面积占比结果和县里提供的 1:5 万调查数据差了近 20%。后来才明白小比例尺数据在大比例尺应用里只能做趋势参考不能做精确统计。现在我的习惯是先看 SHP 元数据里的比例尺再决定用它做什么级别的分析。全国尺度用 1:100 万到 1:400 万省级用 1:50 万到 1:100 万县级用 1:5 万到 1:10 万。如果手头只有小比例尺数据就在报告里明确标注“数据精度限制”并建议后续用高精度数据复核。希望帮到你。本文还有配套的精品资源点击获取