简介：这份中国喀斯特岩溶空间分布矢量数据集面向地理信息、地质地貌与环境规划领域的研究者及从业者，用于分析岩溶地块边界、岩性差异与空间分布格局。资源包共8个文件，约1.2MB，以SHP矢量数据为核心，配套SHX、SBX、SBN空间索引，DBF属性表，PRJ空间参考，CPG编码说明及XML元数据，可直接在ArcGIS、QGIS等软件中加载与检索。属性字段包含岩性分类（连续与不连续碳酸盐岩）、地块面积、周长及岩性文本标签，便于评估岩溶发育程度与地表地下水系统规模。已有394人学习下载，适合用于喀斯特地貌研究、农业规划、水利工程与旅游开发等场景，帮助读者快速获取覆盖多省份的岩溶分布底图与属性信息。

1. 喀斯特岩溶空间分布矢量数据集：从一张 SHP 说起

做西南地区水文、地质灾害或生态评估的同行，大概率都遇到过同一个尴尬：手上有 DEM、有降雨栅格、有县域行政区划边界 shp，唯独缺一层“哪里是喀斯特”的底图。碳酸盐岩分布区和非可溶岩区的产流机制、土壤侵蚀模数、边坡稳定性完全不是一回事，缺了这层约束，后面所有模型都在裸奔。中国 Karsts 喀斯特岩溶空间分布矢量数据集 SHP 数据，解决的正是这件事——它把全国可溶岩分布范围落成面状矢量，能直接叠加到 ArcGIS 或 QGIS 里参与空间分析。这篇笔记面向三类人：要拿它做区划统计的规划人员、要拿它做模型掩膜的科研人员，以及需要把 shp 转成其他格式（shp 转 txt、shp 格式矢量数据导出为 wkt、shp 转 3dtiles）做二次分发的工程人员。下面按“数据长什么样 → 怎么用 → 坑在哪”的顺序讲透。

2. 先搞懂喀斯特岩溶矢量数据的属性结构与坐标基准

拿到一个 shp，第一件事不是急着拖进软件，而是先确认它到底装了什么。喀斯特岩溶空间分布这类专题数据，和常见的行政区划边界 shp 不一样，它的价值一半在几何、一半在属性表。属性表里通常记录了岩性分组、可溶岩类型、分布面积等字段，这些字段决定了你后面能不能按“纯碳酸盐岩 / 碳酸盐岩夹碎屑岩 / 碎屑岩夹碳酸盐岩”分级统计。

2.1 面状要素、属性字段与岩性分级逻辑

喀斯特发育的本质是碳酸盐岩的可溶性，所以数据集的几何形态是面（Polygon），每个面代表一片连续的可溶岩出露区。属性表里常见的字段包括岩性代号、岩石组合类型、面积（一般以平方米或平方公里计）。分级逻辑通常遵循“可溶岩占比”这条主线：连续型碳酸盐岩、互层型、夹层型，可溶性依次递减。你在做区划时如果只按面积统计而不区分岩性等级，结论会失真——同样是喀斯特区，纯灰岩区和白云岩夹层的产流响应差得很远。

判断字段含义有个笨办法但很管用：把属性表导出成 CSV，按岩性字段做频数统计，看类别数是否和文献里的分级数对得上。对不上就说明字段被裁剪或重编码过，需要回头找原始说明。

2.2 坐标系与投影：为什么必须先看 prj 文件

SHP 是 ESRI 定义的一套多文件格式，.shp存几何、.dbf存属性、.shx存索引、.prj存坐标系。很多人只拷贝.shp一个文件，结果打开后坐标是乱的，或者干脆提示“未知空间参考”。喀斯特数据集常见的基准是 CGCS2000 或 WGS84 地理坐标系，也可能是 Albers 等积投影。等积投影对面积统计友好，地理坐标系对叠加在线底图友好，两者不能混用。

提示：如果.prj缺失，不要凭经验硬猜。先用ogrinfo读一遍几何范围，再结合已知的经纬度范围反推坐标系，比直接定义投影安全得多。

2.3 用 GDAL 快速体检一个喀斯特 SHP 文件

下面这段命令不依赖任何桌面软件，纯命令行就能把数据的基本情况摸清楚，适合批量处理前做质检。

# 查看图层概要：几何类型、要素数、坐标系、字段定义 ogrinfo -so -al china_karst.shp # 输出属性表前若干行，确认岩性字段的实际取值 ogrinfo -al -where "1=1" china_karst.shp | head -n 40 # 把坐标系单独打印出来，确认是否为等积投影 gdalsrsinfo china_karst.prj

-so表示只输出摘要不读几何，速度快；-al表示遍历所有图层；gdalsrsinfo会打印 EPSG 代码和投影参数。如果ogrinfo报“Unable to open datasource”，八成是.shx或.dbf没跟着一起放。参数上，-where后面接的是 OGR SQL 语法，不是标准 SQL，字段名带中文时要用双引号包住。

3. 把喀斯特 SHP 接进实际分析流程的三种做法

数据体检通过后，就进入真正的使用环节。不同岗位的诉求差别很大：做区划的要按行政单元统计面积，做模型的要把它转成栅格掩膜，做可视化的要把它转成 3D 瓦片。下面三种做法覆盖绝大多数场景。

3.1 按县域行政区划边界做面积统计

这是最典型的“喀斯特 + 行政区划”叠加需求。思路是：把喀斯特面与县域面做相交（Intersect），再按县代码分组求和。用 GeoPandas 几行就能跑完，比在桌面软件里点来点去可复现得多。

import geopandas as gpd # 读取喀斯特分布与县域边界，统一到同一投影再做叠加 karst = gpd.read_file("china_karst.shp").to_crs(epsg=4526) county = gpd.read_file("county_boundary.shp").to_crs(epsg=4526) # 相交后按县代码聚合面积，面积单位取决于投影，4526 下为平方米 inter = gpd.overlay(karst, county, how="intersection") inter["area_km2"] = inter.geometry.area / 1e6 stat = inter.groupby("county_code")["area_km2"].sum().reset_index() stat.to_csv("karst_area_by_county.csv", index=False)

关键点在to_crs：两个图层必须投影一致才能 overlay，否则几何运算结果会错得离谱。epsg=4526是 CGCS2000 高斯投影的一个分带，实际用哪个分带取决于你的研究区经度，跨带时建议改用 Albers 等积投影。how="intersection"只保留相交部分，如果你还想保留没有喀斯特分布的县，改成how="union"再补零。

3.2 转成栅格掩膜喂给水文或生态模型

模型通常吃栅格，所以要把矢量面栅格化。用 GDAL 的gdal_rasterize最省事，指定一个参考栅格就能保证像元对齐。

# 以 DEM 为模板，把喀斯特面烧录成 1/0 掩膜 gdal_rasterize -a "karst_flag" -init 0 -burn 1 \ -te 100.0 25.0 110.0 35.0 -tr 0.001 0.001 \ -ot Byte -l china_karst china_karst.shp karst_mask.tif

-a指定用哪个属性字段的值，-init 0把背景填 0，-burn 1把面内烧成 1。-te和-tr定义范围和分辨率，如果已有 DEM，更稳妥的做法是用-te和-tr抄 DEM 的参数，避免像元错位。-l指定图层名，SHP 单图层时也不能省，否则部分版本会报错。

3.3 导出为 WKT 或 TXT 做跨平台分发

有些平台不吃 SHP，只认文本坐标。把面要素导出成 WKT 或逐点 TXT 是常见做法，注意面要素导出 WKT 时是POLYGON((...))结构，环的闭合和顺序不能乱。

import geopandas as gpd gdf = gpd.read_file("china_karst.shp") # 导出 WKT，保留属性字段，方便对方按字段过滤 gdf["wkt"] = gdf.geometry.apply(lambda g: g.wkt) gdf[["karst_type", "wkt"]].to_csv("karst_wkt.csv", index=False) # 导出逐点 TXT，适合只吃坐标串的老系统 with open("karst_points.txt", "w", encoding="utf-8") as f: for geom in gdf.geometry: for x, y in geom.exterior.coords: f.write(f"{x:.6f},{y:.6f}\n")

WKT 导出保留了拓扑信息，适合对方重建面；TXT 逐点导出丢了环的边界，只适合点云式消费。坐标精度.6f对应约 0.1 米，做区域尺度分析足够，再高就是浪费存储。

4. 避坑与排查：喀斯特 SHP 处理中最容易翻车的五件事

这一节全是血泪经验，每条都按“现象 → 原因 → 解决”写，遇到问题直接对号入座。

现象一：打开后要素全挤在一个点上，或者跑到非洲西海岸去了。原因：.prj缺失或坐标系定义错误，软件按默认 WGS84 解释，而数据实际是投影坐标。解决：先用gdalsrsinfo看.prj，缺失时根据坐标数值范围反推——数值在 6 到 7 位数通常是投影坐标，在 -180 到 180 之间是地理坐标，确认后用gpd.read_file(...).set_crs(epsg=xxxx)补上再转。

现象二：按县统计面积，结果比县域总面积还大。原因：喀斯特面与县域面投影不一致，overlay 时几何被错误拉伸；或者喀斯特数据本身有重叠面。解决：overlay 前强制to_crs统一投影，并用gdf.geometry.is_valid检查无效几何，无效的用buffer(0)修复后再算。

现象三：栅格化后掩膜边缘锯齿严重，和矢量边界对不上。原因：gdal_rasterize没指定-te和-tr，默认按数据范围生成，像元格与 DEM 错位。解决：始终用参考栅格的-te -tr参数，必要时加-at开启抗锯齿，但注意-at会让边缘像元出现中间值，做二值掩膜时反而要慎用。

现象四：属性表里岩性字段全是乱码。原因：.dbf的字符编码不是 UTF-8，常见于老数据。解决：用gpd.read_file("xxx.shp", encoding="gbk")试读，或先用ogrinfo看字段值，确认编码后再读。写回时统一用 UTF-8，避免下游再踩。

现象五：shp 转 3dtiles 后模型悬空或贴地错误。原因：SHP 是二维面，没有高程，转 3D 时高度字段缺失或单位不一致。解决：先给面要素补一个高程字段（可采样 DEM 得到），再走“面拉伸成体”的流程，拉伸高度单位要和坐标系水平单位一致，投影坐标下是米，地理坐标下是度，混用必翻车。

5. 进阶：用渔网分割做喀斯特分布密度制图

前面讲的都是“把喀斯特数据当输入”，这一节换个角度，把它当被统计对象，做一张能直接进报告的密度图。核心思路是渔网分割 shp：生成规则格网，统计每个格网内喀斯特面积占比，再按占比分级设色。这个方法比按行政区统计更客观，不受边界形状影响，适合做全国或流域尺度的连续制图。

import geopandas as gpd import numpy as np from shapely.geometry import box karst = gpd.read_file("china_karst.shp").to_crs(epsg=4526) minx, miny, maxx, maxy = karst.total_bounds cell = 50000 # 50 公里格网，按研究区尺度调整 # 生成渔网格网 cols = np.arange(minx, maxx, cell) rows = np.arange(miny, maxy, cell) cells = [box(x, y, x + cell, y + cell) for x in cols for y in rows] grid = gpd.GeoDataFrame(geometry=cells, crs=karst.crs) # 相交后按格网聚合，计算喀斯特面积占比 inter = gpd.overlay(grid, karst, how="intersection") inter["k_area"] = inter.geometry.area inter["cell_id"] = inter.geometry.apply( lambda g: (int((g.bounds[0] - minx) // cell), int((g.bounds[1] - miny) // cell)) ) agg = inter.groupby("cell_id")["k_area"].sum().reset_index() agg["ratio"] = agg["k_area"] / (cell * cell) # 回写格网几何，导出制图 grid["cell_id"] = list(zip( ((grid.geometry.bounds.minx - minx) // cell).astype(int), ((grid.geometry.bounds.miny - miny) // cell).astype(int) )) result = grid.merge(agg, on="cell_id", how="left").fillna({"ratio": 0}) result.to_file("karst_density_grid.shp", encoding="utf-8")

格网大小cell是最关键的参数：太小则大量格网为空、图面碎；太大则抹平空间差异。经验值是研究区短边的 1/50 到 1/100，全国尺度取 50 公里起步比较稳。cell_id用格网左下角坐标反算，保证 overlay 前后能对上，这一步如果直接用索引对齐，相交后要素数变化会导致错位。分级设色时建议用自然断点法而不是等间距，喀斯特占比的分布本身是偏态的，等间距会把大部分格网压进最低一级。

最后说个我自己的习惯：每次拿到新的喀斯特 SHP，先跑一遍ogrinfo -so -al，把要素数、坐标系、字段名记在一张便签上贴在项目 README 里。这个动作花不到一分钟，但能省掉后面无数次“为什么叠加不上”的返工。数据这东西，后悔药从来都是提前吃的。希望帮到你。

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