
简介这份中国喀斯特岩溶空间分布矢量数据集面向GIS、地理学与地质科研人员及学生用于岩溶地貌区域划分、溶蚀作用模拟与地质灾害监测等空间分析场景。资源包共8个文件约1.2MB以SHP矢量格式存储包含shp图形数据、dbf属性表、prj坐标系统定义、shx索引及cpg、sbn、sbx、shp.xml等配套文件可在主流GIS软件中直接读取。数据以面状Polygon记录岩溶地块边界属性字段含rock_type岩性分类连续与不连续碳酸盐岩、Shape_Area与Shape_Len面积周长、RTypeLabel岩性文本标签便于快速筛选与统计。已有263人学习下载适合需要一手空间数据开展岩溶地貌研究、制图与建模的读者也可为土地资源管理和生态保护提供参考依据。1. 喀斯特岩溶空间分布矢量数据集从一张 SHP 到一套可复现的岩溶分析底图做西南地区水文或工程地质的同行多半遇到过这种局面手头有一份 DEM、一份降雨栅格想叠加岩溶发育程度做分区评价结果发现岩溶边界只能靠文字描述在图上手描。中国喀斯特岩溶空间分布矢量数据集 SHP 数据解决的正是这个“底图从哪来”的问题。它把碳酸盐岩出露与岩溶发育的空间范围落成面状矢量字段里通常带岩性类型或发育等级能直接进 ArcGIS、QGIS 或 PostGIS 做叠加、裁剪、统计。适合做区域地质调查、岩溶塌陷易发性评价、隧道选线避让、地下水脆弱性制图的从业者。这篇笔记按“数据长什么样 → 怎么加载和检查 → 怎么和业务数据叠加 → 坑在哪 → 怎么进阶验证”的顺序讲新手能照着跑通熟手能直接看参数和边界条件。2. 拿到 SHP 先别急着画图坐标系、字段与拓扑的三步体检喀斯特岩溶空间分布这类数据最怕的是“看着对、算出来错”。坐标系没对齐、字段类型不对、面要素自相交都会让后面的叠加分析结果变成玄学。我一般拿到任何一份岩溶 SHP先做三件事确认坐标系与投影、读字段结构和属性分布、查几何有效性。这三步做完才决定要不要重投影、要不要修拓扑。2.1 用 GeoPandas 读元数据确认 CRS 和字段类型import geopandas as gpd # 读取喀斯特岩溶空间分布 SHP注意 encoding 在中文属性下常需指定 gdf gpd.read_file(karst_distribution.shp, encodingutf-8) print(要素数量:, len(gdf)) print(坐标系:, gdf.crs) print(几何类型:, gdf.geom_type.unique()) print(字段与类型:) print(gdf.dtypes) print(属性前 5 行:) print(gdf.drop(columnsgeometry).head())这段代码的逻辑是先看数据规模再看坐标系是不是地理坐标EPSG:4326还是投影坐标如 CGCS2000 高斯投影。参数上encoding在属性含中文时建议显式指定否则容易读出乱码gdf.crs若返回 None说明缺少.prj文件必须找原始说明补上不能猜。字段类型里要重点看岩性编码字段是整型还是字符串后面做分类统计时处理方式不同。2.2 检查几何有效性与面积分布识别碎面和自相交# 几何有效性检查 invalid gdf[~gdf.is_valid] print(无效几何数量:, len(invalid)) # 面积分布判断是否存在异常碎面 gdf_proj gdf.to_crs(epsg4547) # CGCS2000 / 3-degree Gauss-Kruger zone 39 gdf_proj[area_km2] gdf_proj.area / 1e6 print(gdf_proj[area_km2].describe()) # 按岩性字段统计面积 if rock_type in gdf_proj.columns: print(gdf_proj.groupby(rock_type)[area_km2].sum().sort_values(ascendingFalse))逻辑说明先判断无效几何再投影到等面积或等距投影算面积避免用经纬度直接算面积导致数量级错误。参数上EPSG:4547 适用于中国中东部 3 度带西南地区要根据经度选对应带号选错会让面积偏差百分之几到十几。按岩性汇总面积能快速判断数据是否合理比如纯碳酸盐岩面积远大于总喀斯特区面积就说明字段或范围有问题。2.3 拓扑修复与字段规范化为后续叠加做准备from shapely.validation import make_valid # 修复无效几何 gdf[geometry] gdf[geometry].apply( lambda geom: make_valid(geom) if not geom.is_valid else geom ) # 统一岩性字段命名便于后续脚本复用 rename_map {岩性: rock_type, 类型: rock_type, 发育程度: karst_level} gdf gdf.rename(columns{k: v for k, v in rename_map.items() if k in gdf.columns}) # 导出修复后的数据保留原始字段 gdf.to_file(karst_clean.shp, encodingutf-8)这里用make_valid处理自相交和环方向错误比buffer(0)更稳不会把细长面压没。字段重命名是为了后续脚本不依赖中文列名减少编码翻车。导出时保留.shp格式注意 Shapefile 单文件字段名上限 10 个字符长字段名会被截断必要时改用 GeoPackage。提示Shapefile 对字段名长度和中文支持有限若属性字段多或含长中文名建议同时导出一份 GeoPackage 作为工作副本。3. 把岩溶 SHP 接进业务分析叠加、裁剪与分区统计的完整链路数据体检通过后真正产生价值的是把它和业务图层叠加。常见场景有三类和行政区划叠加算各县岩溶面积占比、和 DEM 或坡度叠加做发育程度分区、和工程线路叠加做避让分析。这一章按“叠加 → 裁剪 → 统计 → 出图”的链路走每步给可抄的代码和参数说明。3.1 与行政区划叠加用空间连接算县域岩溶面积import geopandas as gpd karst gpd.read_file(karst_clean.shp).to_crs(epsg4547) county gpd.read_file(county_boundary.shp).to_crs(epsg4547) # 空间连接每个岩溶面落到所属县 joined gpd.sjoin(karst, county, howinner, predicateintersects) # 按县汇总岩溶面积 joined[area_km2] joined.area / 1e6 result joined.groupby(county_name)[area_km2].sum().reset_index() result[ratio] result[area_km2] / county.set_index(county_name).area * 100 print(result.sort_values(area_km2, ascendingFalse).head(10))逻辑是先统一投影再做空间连接predicateintersects比within更宽容能处理边界压线的情况。参数上howinner只保留有岩溶的县若要保留全部县用left。面积汇总前要确保两个图层投影一致否则连接结果会错位。占比计算用县总面积做分母注意县边界若有重叠或空洞分母会偏大。3.2 按坡度分级裁剪生成岩溶发育程度分区import rasterio from rasterio.mask import mask import geopandas as gpd import numpy as np karst gpd.read_file(karst_clean.shp).to_crs(epsg4547) with rasterio.open(slope.tif) as src: slope src.read(1) transform src.transform # 用岩溶面裁剪坡度栅格 out_image, out_transform mask(src, karst.geometry, cropTrue) out_image out_image[0] # 按坡度分级统计岩溶面积 bins [0, 8, 15, 25, 35, 90] labels [平缓, 较缓, 中等, 较陡, 陡峭] classes np.digitize(out_image, bins) - 1 for i, label in enumerate(labels): count np.sum(classes i) print(f{label}: {count} 像元)这段代码用rasterio.mask按岩溶边界裁剪坡度栅格cropTrue会缩小输出范围节省内存。参数上坡度分级阈值按项目规范调整西南岩溶区常用 8°、15°、25°、35° 作为分界。np.digitize返回的索引从 1 开始减 1 后对应标签。注意裁剪后像元值可能含 NoData统计前要过滤。3.3 导出分区结果与制图字段方便 ArcGIS 直接出图# 将分级结果写回矢量按主导坡度等级赋属性 from shapely.geometry import shape import rasterio from rasterio.mask import mask # 简化示例按岩溶面与坡度等级求交后赋等级 karst[slope_level] 中等 # 实际应按空间统计结果赋值 karst.to_file(karst_slope_class.shp, encodingutf-8) # 导出 WKT 便于入库或跨平台使用 karst[wkt] karst.geometry.apply(lambda g: g.wkt) karst[[rock_type, slope_level, wkt]].to_csv(karst_wkt.csv, indexFalse)导出时保留rock_type和slope_level两个关键字段方便在 ArcGIS 里按符号系统直接渲染。WKT 导出适合入库或传给前端做轻量展示注意 WKT 字符串较长CSV 读取时留意字段截断。若后续要做三维展示可在此基础上转 3dtiles但岩溶面数据转三维通常需要拉伸或贴地处理不是直接转换。注意叠加分析前务必确认所有图层的投影一致经纬度与投影坐标混用是岩溶面积统计最常见的翻车点。4. 喀斯特 SHP 使用中的避坑与排查五个真实踩坑记录这一章按“现象 → 原因 → 解决”写五条我实际遇到过的坑覆盖坐标系、字段、拓扑、性能和格式转换。4.1 面积算出来偏大或偏小一个数量级现象用 GeoPandas 直接算面积结果和 ArcGIS 差很多。原因数据是地理坐标系EPSG:4326gdf.area按度算不是平方米。解决先to_crs到投影坐标系再算面积西南地区按经度选 3 度带或 6 度带不确定就用等面积投影如 Albers。4.2 中文属性读出来是乱码现象read_file后岩性字段显示为问号或方块。原因Shapefile 的.dbf编码未指定常见为 GBK 或 UTF-8。解决读取时加encodinggbk或encodingutf-8试导出时统一用 UTF-8长期方案是转 GeoPackage。4.3 空间连接后要素数量暴增现象sjoin后行数远大于原始岩溶面数量。原因一个岩溶面跨多个县或县边界有重叠导致一对多连接。解决先检查县边界是否有重叠必要时用dissolve合并统计时按唯一 ID 去重或改用predicatewithin只保留完全落入的要素。4.4 裁剪栅格时内存溢出现象用大范围岩溶面裁剪高分辨率 DEM 时进程被 kill。原因mask默认读全图再裁内存不够。解决先用cropTrue缩小范围或分块读取也可先用岩溶面做clip再统计避免整图加载。4.5 转 GeoPackage 后字段丢失现象SHP 转 GPKG 后部分字段为空。原因Shapefile 字段名截断或类型不兼容转换时映射失败。解决转换前重命名字段为短英文名检查字段类型转换后用gpd.read_file回读验证字段完整性。5. 进阶验证与复用用岩溶 SHP 做一套可复现的易发性底图数据用顺之后真正拉开差距的是验证和复用。我一般会做两件事一是用已知岩溶塌陷点做空间验证看岩溶分布与灾害点的重合率二是把整套处理流程脚本化换一个区域只需改路径和投影参数。验证时用sjoin把灾害点落到岩溶面上统计落入比例若比例明显偏低要回头检查岩溶边界是否偏保守或灾害点坐标是否有偏移。脚本化时把投影、字段映射、分级阈值抽成配置字典避免每次手改代码。验证项方法合格参考坐标系一致性对比 CRS 与项目规范全部图层同投影几何有效性is_valid检查无效几何占比 1%面积合理性与文献或统计年鉴对比偏差在合理范围灾害点重合率空间连接统计根据区域经验判断字段完整性回读检查关键字段无缺失最后说个习惯我每次拿到新的岩溶 SHP都会先裁一小块跑通全流程再放大到全域。这样翻车成本低也容易定位是数据问题还是脚本问题。希望帮到你。本文还有配套的精品资源点击获取