简介中国Karsts喀斯特岩溶空间分布矢量数据集面向GIS地理信息研究人员、规划人员与环境科学家提供了一份覆盖多个省份和地区的岩溶地块面状数据。数据采用Polygon矢量格式每个地块记录岩性分类连续/不连续碳酸盐岩、面积、周长及岩性文本标签可支撑喀斯特地貌发育程度分析、地表与地下水系统规模评估以及水利工程、旅游开发等应用场景。资源包为ZIP压缩格式共8个文件整体大小约1.2MB。其中SHP存储空间几何DBF承载属性表PRJ定义坐标参考XML提供元数据SBN/SBX建立空间索引以提升检索效率CPG处理字符编码SHX关联几何与属性记录结构完整可直接导入ArcGIS、QGIS等常见GIS平台使用。已有388人学习下载适合地理信息系统应用、地质环境研究及区域规划等领域的中高级用户作为基础数据支撑。1. 数据到手先别急着出图喀斯特SHP要过的第一道关做好西南岩溶区划、地下水评价的同行大概率卡在同一个环节搜到一套“中国Karsts喀斯特岩溶空间分布矢量数据集SHP数据”下载、解压、拖进GIS结果属性表乱码、坐标系是经纬度、图斑边缘还叠着不少重复区。这套数据本身不复杂把它理解成“以点、线、面三种几何结构写入ESRI Shapefile的岩溶空间分布描述”配合属性表里的类型代码和面积字段用来做裁剪、统计和出图。常见的问题反而在拿到之后先读字段和投影再决定转不转KML或GeoJSON然后做裁剪和渔网统计最后一关是交付前检查。下面按这个流程走一遍适合做区划、制图和资源评价的GIS从业者参考。2. 数据包里装的是什么先读图层再谈统计这个标题下的SHP数据交付形态通常是一个文件夹包含.shp、.dbf、.shx、.prj四个必要文件有时多一个.cpg。.shp记录几何坐标.dbf记录属性表.shx是空间索引.prj写坐标系描述.cpg写字符编码。四个文件里只要少一个就会出现“缺字段”“无法定位”“打开全黑”这类现象。SHP格式能沿用这么多年结构简单是原因之一但简单也意味着容错低少一个伴随文件就会出问题。2.1 点、线、面三种图层结构先分清手里到底是哪一种几何。喀斯特岩溶空间分布数据集并不是只有一张面图层常见做法是把三类要素分开存放面图层岩溶发育分区、岩溶盆地、峰丛洼地分布区属性表里一般有类型代码和区域名称线图层地下暗河、断层、岩溶管道走向属性表里一般有等级或类型字段点图层溶洞、泉点、天坑、竖井坐标点属性表里一般有名称、出露高程。判断一个SHP文件是什么几何类型最直接的办法是用GDAL/OGR看概览一行命令ogrinfo -al -so karst.shp输出里的Feature Count是要素个数Geometry字段会标出Polygon、LineString或Point。这一步不要跳过我见过有人把地下河线图层当面图层统计面积结果面积列全是0整张统计表作废重来。想进一步知道属性表列名用Python读字段定义from osgeo import ogr ds ogr.Open(karst.shp, 0) layer ds.GetLayer() defn layer.GetLayerDefn() for i in range(defn.GetFieldCount()): fld defn.GetFieldDefn(i) print(fld.GetName(), fld.GetTypeName())参数0表示只读打开不会改动原始文件。字段类型常见有OFTString、OFTInteger、OFTReal分别对应文本、整型、浮点数。喀斯特数据集的字段名各地版本差别很大有type、code、name、area也有LSD、KSFL这类拼音缩写有的字段干脆全空。先确认“哪一列是类型、哪一列是面积”后面过滤和统计才不糊涂。如果想把每个类型的图斑数、面积总量一次性摸清接着跑一段统计from collections import Counter type_counter Counter() area_sum {} for feat in layer: t feat.GetField(type) type_counter[t] 1 area_sum[t] area_sum.get(t, 0.0) feat.GetField(area_km2) for t, cnt in type_counter.items(): print(t, cnt, round(area_sum[t], 2))这段代码不做几何运算只读.dbf属性字段速度快。为什么要有这一步很多公开下载的矢量数据集图层名和字段说明并不严格对应发布者“字段填错了”的情况不少用计数结果和说明文件交叉验证能提前发现字段挂错、单位写错这类黑匣子问题。提示ogrinfo命令找不到时安装gdal-bin或者直接在QGIS数据源管理器里读属性表效果一样。2.2 先解决投影面积和叠加全看它SHP数据的坐标系常见就两类。一类是WGS84经纬度.prj里写GEOGCS[WGS 84]另一类是投影坐标系多写CGCS2000或西安80的高斯-克吕格带中央经线和带号。“先查投影”是处理一切外部SHP的第一原则因为后续所有面积计算、裁剪、叠加全部依赖坐标系正不正确。判断命令gdalsrsinfo karst.shp输出会直接显示坐标系名称。如果是经纬度坐标系坐标单位是度不能直接算面积也不能直接跟投影坐标系的省界做叠加。需要先把图层重投影到等积投影上再继续干活。以西南喀斯特区为例我通常转到阿尔伯斯等积圆锥投影import os from osgeo import ogr, osr src_ds ogr.GetDriverByName(ESRI Shapefile).Open(karst_wgs84.shp, 0) src_lyr src_ds.GetLayer() src_srs src_lyr.GetSpatialRef() target_srs osr.SpatialReference() target_srs.ImportFromEPSG(4555) # 示例CGCS2000/Albers供西南区域参考 target_srs.SetProjParm(central_meridian, 105.0) out_path karst_albers.shp driver ogr.GetDriverByName(ESRI Shapefile) if os.path.exists(out_path): driver.DeleteDataSource(out_path) out_ds driver.CreateDataSource(out_path) out_lyr out_ds.CreateLayer(karst_albers, target_srs, ogr.wkbMultiPolygon) out_lyr.CreateField(ogr.FieldDefn(type, ogr.OFTString)) out_lyr.CreateField(ogr.FieldDefn(area_km2, ogr.OFTReal)) for src_ftr in src_lyr: geom src_ftr.GetGeometryRef().Clone() geom.TransformTo(target_srs) out_ftr ogr.Feature(out_lyr.GetLayerDefn()) out_ftr.SetGeometry(geom) out_ftr.SetField(type, src_ftr.GetField(type)) out_ftr.SetField(area_km2, geom.GetArea() / 1_000_000) out_lyr.CreateFeature(out_ftr) out_ds None src_ds None参数解释ImportFromEPSG(4555)是我在西南地区常用的CGCS2000阿尔伯斯等积投影示例不要在所有地区照抄。central_meridian是中央经线西南喀斯特集中区习惯用105华南地区可以按项目需要改成108或111。投影坐标系的几何面积单位是平方米代码里除以1e6转成平方公里这是面积字段最容易错的地方。改用EPSG:6933之类的世界等积投影也行关键是“等积”二字面积统计类的活必须用等积投影不能用墨卡托或经纬度直接算。如果打开图层时软件提示“未知投影”多半是.prj缺失。不要用记事本手写一个.prj去猜先找原始说明文件里的坐标描述找不到就只能用控制点做配准或根据已知坐标反推投影。这是SHP使用里最典型的玄学现场——图看着没问题但量算全是错的。3. 把SHP转成下游能用的格式KML、GeoJSON和CSVDWG、Excel、SHP、KML、GeoJSON这些格式在喀斯特项目中经常混合出现。地质队给的CAD图先转成SHP环保部门的Excel点位表也要生成SHP最后交到在线平台又得从SHP转KML或GeoJSON。这是整个流程里耗时最多的一段。3.1 用ogr2ogr完成shp转KML和转GeoJSON两步最常见的转换命令是ogr2ogr -f KML karst.kml karst_albers.shp ogr2ogr -f GeoJSON -t_srs EPSG:4326 karst.geojson karst_albers.shp第一条把SHP转成KMLKML面向谷歌地球和手机地图坐标必须WGS84ogr2ogr会按目标格式自动重投影。第二条转GeoJSON时加了-t_srs EPSG:4326强制输出WGS84经纬度如果不加GeoJSON会保留源坐标系的数值放到Leaflet这类前端地图里就偏位。-f指定输出格式-t_srs指定目标坐标系这两个参数做格式转换时优先确认。如果手头只有dwg格式的地质图就要先在GIS里把dwg转成shp再走上面命令。CAD里的多段线转过来可能是线不是面要检查Geometry类型需要合成面时用buffer或polygonize处理后才能继续做面积统计。这个环节出错的概率不高但一旦出错后面所有裁剪结果都跟着错。还要提一句“shp转txt/csv”的需求。SHP不是所有人都能直接打开把属性表导成CSV给不了解GIS的同事做汇总统计很常见ogr2ogr -f CSV karst.csv karst.shp -lco GEOMETRYAS_XY-lco GEOMETRYAS_XY表示把几何坐标以X、Y两列的形式写到CSV里适合只核对属性、不关心拓扑的场景。如果不需要坐标列去掉-lco直接输出纯属性表。注意CSV导出的中文在Excel打开时可能出现乱码多数情况下需要加-lco ENCODINGUTF-8或事后用文本编辑器转成带BOM的UTF-8。批量转换时写个循环更快for f in *.shp; do ogr2ogr -f GeoJSON -t_srs EPSG:4326 ${f%.shp}.geojson $f done这段循环把当前目录下所有SHP全部转成GeoJSON${f%.shp}是去掉后缀的字符串处理输出文件名会自动变成xx.geojson。用在批量交付前整理数据很合适。3.2 用Excel点表生成SHP经纬度清点和编码很多岩溶点数据、泉点调查表都以Excel形式出现对应“excel点转shp”的常见需求。步骤如下确认经纬度列没有混入文本删除缺坐标的行用GeoPandas生成点要素最终写SHP。import pandas as pd import geopandas as gpd from shapely.geometry import Point df pd.read_excel(karst_points.xlsx, sheet_name泉点) df df.dropna(subset[经度, 纬度]) df[geometry] df.apply( lambda row: Point(row[经度], row[纬度]), axis1 ) gdf gpd.GeoDataFrame(df, crsEPSG:4326) gdf.to_file(karst_points.shp, encodingutf-8)参数说明dropna(subset[...])会删除经纬度任一为空的整行避免产生坐标(0,0)之类的野点。crsEPSG:4326声明源坐标是WGS84经纬度Excel里的十进制度数坐标必须对应这个坐标系。.to_file(encodingutf-8)决定.dbf属性表的编码。Excel里最常翻车的是经纬度写成“25°30′30″”这种度分秒文本还有把度分秒拆成三列的情况。上面这段代码只支持十进制度数遇到文本格式要先在Excel里用公式转十进制或写个正则预处理。另一个注意点是中文属性列名utf-8的.dbf在ArcGIS里打开是乱码要么改英文列名要么把encoding改成GBK。我一般统一写成英文列名加GBK输出gdf gdf.rename(columns{泉点名称: name, 出露高程: elev}) gdf.to_file(karst_points.shp, encodingGBK)写成GBK后ArcGIS和QGIS都能正确显示中文字段内容。缺点是属性表列名必须是单字节字符所以上面用了name、elev。如果不想写代码QGIS里“新建SHP文件”也可以做同样的事新建点图层选WGS84坐标系添加name和elev字段然后用“粘贴要素”把Excel坐标粘贴进去。这种方式适合几十个点的量数据量大还是GeoPandas快。新建SHP文件时几何类型选错会直接影响后面能不能用点数据就选Point岩溶分区这类区域数据选Polygon不要图省事全选Point。4. SHP使用避坑五个现象、原因和现场解法4.1 打开SHP以后图层是空的现象下载的SHP添加进ArcGIS后缩放到图层还是一张空白画布属性表里却能看到几百条记录。原因.shx索引文件缺失或.shp文件没解压完整坐标部分写坏了。解决先到文件夹确认四个伴随文件都在再查要素数量Feature Count大于0但图形不显示大概率是几何无效。用GDAL做一次修复ogr2ogr -makevalid -skipfailures karst_fixed.shp karst_broken.shp-makevalid重建无效几何-skipfailures跳过修不动的坏要素避免整批转换中断。如果这个命令也报“corrupted”之类的错原始数据源问题太大别在损坏文件上继续耗回头重新下载或换一个版本。SHP文件是老格式单文件体积超过2GB、要素数量极多时也容易出问题这类情况不叫损坏是格式瓶颈。遇到超大矢量数据转成GeoPackage再处理更稳ogr2ogr -f GPKG karst_fixed.gpkg karst.shpGeoPackage把几何、属性、坐标系封在一个文件里规避了SHP多文件缺漏的问题也适合在移动端和FME流程里替代SHP。4.2 属性表汉字乱码字段还缺一半现象.dbf里的中文全变问号字段只有两三个原始表里的中文列名也没了。原因SHP的.dbf编码本身就缺乏标准说明同一个文件在不同软件里解读不一样。QGIS按UTF-8读ArcGIS按系统区域设置读GBK和UTF-8互相换着解析自然乱码。字段缺一半则可能是早期DBF字段名10字节限制导致截断常见于年代较早的岩溶数据。解决在QGIS的数据源管理器里手动切换UTF-8或GBK重新加载字段截断的只能根据说明文件恢复或重做属性挂接。非中文字段名是规避这系列问题最省事的手段。另外持shapechk这类修复工具可以重建.shp索引但要先确认真正问题是缺失文件还是文件损坏否则修完还是打不开。4.3 手改.prj导致整个图层坐标系错乱现象为了让两个图层“对齐”有人把.prj用记事本改成WGS84结果图还是偏甚至比原来偏得更厉害。原因.prj只描述坐标系不负责把坐标数值转换到另一个坐标系。你改了坐标系的名坐标值还是原来的投影坐标等于把米当度用。解决要改坐标系就用重投影工具不要改描述文件。前面2.2节的TransformTo就是正确做法。如果只是想让两个图层临时叠加ArcGIS里可以开“on the fly”投影QGIS里设置项目CRS不影响底层数据。4.4 用经纬度图层直接算面积数字离谱现象属性表里算面积得到0.000387这种数明显不对。原因WGS84经纬度下几何单位是度面积单位是平方度平方公里无法直接换算因为球面面积随纬度变化。解决先转等积投影再算面积。代码里geom.GetArea()得到的单位是投影坐标单位平方米再除以1e6就是平方公里。这一步做完再核对一下区域总面积和公开口径出入不大才算过。如果坚持在经纬度图层里看面积可以用QGIS的area($geometry)表达式但得到的结果是平方度不能直接当面积用这是公式层面的限制不是软件bug。4.5 多边形重叠导致面积重复统计现象喀斯特岩溶分布图斑中同一区域内不同年份的图斑边界部分重叠Dissolve之前按类型求和会重复计算。原因数据由多个单位多次矢量化或图斑之间本身存在包含关系。解决先对同类型图斑做Dissolve合并相邻多边形再重新计算面积。GeoPandas的写法karst gpd.read_file(karst_albers.shp) dissolved karst.dissolve(bytype, aggfuncsum) dissolved[area_km2] dissolved.geometry.area / 1_000_000 dissolved.to_file(karst_dissolved.shp, encodingGBK)dissolve(bytype)按类型字段合并aggfuncsum对其他数值字段求和最后重算面积。如果不先剔除重叠后面渔网统计和区块汇总结果都会偏高这就是“统计越算越大”的根源。5. 让岩溶分布数据下场裁剪、渔网统计和属性清理数据转成能用的格式之后要做的事情就具体了。最常接到的需求只要某一流域或行政区范围内的岩溶区按固定网格统计岩溶面积占比把Excel里的岩溶点属性挂到已有SHP上。5.1 用行政区划或流域边界裁剪岩溶SHP先统一坐标系再做裁剪。常见错误是直接把WGS84经纬度的岩溶数据扔到CGCS2000投影的流域边界里做空间操作结果图偏出几十公里。下面的代码前提是两个GeoDataFrame已经采用同一投影import geopandas as gpd karst gpd.read_file(karst_albers.shp) boundary gpd.read_file(流域边界.shp) # 建议与karst统一到同一投影坐标系 target boundary[boundary[name] 南盘江流域] if name in boundary.columns else boundary result gpd.overlay(karst, target, howintersection) result[area_km2] result.geometry.area / 1_000_000 result.to_file(karst_clip_result.shp, encodingGBK)gpd.overlay(howintersection)只保留落在目标范围之内的图形等价于GIS里的Clip。与ArcGIS裁剪相比GeoPandas会把边界处相交的多边形切开避免生成重叠图斑。裁剪后原始属性字段会保留项目方常关心的“分类”字段通常还在但源面积字段的数值还是裁剪前的一定要用geometry.area重新计算面积。流域边界这类基础shp可以是南盘江、淮河、塔里木河流域这类水文分区数据只要边界和岩溶数据处于同一坐标系即可不需要额外转换。如果只裁剪不重投影矢量结果能出图但量算就会翻车。建议在执行gpd.overlay前后各做一次geometry.area对比验证总量变化是否在预期范围内。裁剪出来的结果如果图斑边界出现锯齿状缺口多半是边界数据精度不够或两个数据源的容差设置不一致。5.2 渔网分割按规则网格统计岩溶面积另一个高频需求是“把岩溶分布切到10公里网格上统计每个网格内岩溶面积”。这就是常说的渔网分割。栅格化的TIF可以统计但矢量成果便于二次编辑。生成渔网的常见做法import geopandas as gpd from shapely.geometry import box karst gpd.read_file(karst_dissolved.shp) minx, miny, maxx, maxy karst.total_bounds cell 10_000 # 单位依赖坐标系投影坐标系下为米 cols range(int(minx), int(maxx), cell) rows range(int(miny), int(maxy), cell) grid_polys [] for y in rows: for x in cols: grid_polys.append(box(x, y, x cell, y cell)) grid gpd.GeoDataFrame({geometry: grid_polys}, crskarst.crs) grid[grid_id] range(len(grid)) intersection gpd.overlay(grid, karst, howintersection) stats intersection.dissolve(bygrid_id)[geometry].area / 1_000_000代码按岩溶区外接矩形生成方格box(x, y, xcell, ycell)构造正方形最后做相交并汇总面积。cell变量完全按需调整网格越小格子数量越多overlay计算越慢。数据量大时我建议先按区域分块处理再合并避免一次生成几十万个多边形把内存占满。这里有个容易忽略的点若图斑之间仍存在重叠渔网统计前必须先行处理掉否则每个网格内叠加了多层重叠图斑面积虚增。所以说第4章的Dissolve不只是面积统计步骤也是渔网统计的前置步骤。属性挂接这类清理工作也常在这一步做。老版本的岩溶SHP经常缺少“岩溶类型中文名”只有编码数字。可以用Excel表按编码列merge回填import pandas as pd code_dict pd.read_excel(类型编码表.xlsx) karst karst.merge(code_dict, left_oncode, right_oncode, howleft)merge前必须确认编码列在两个文件里的类型一致比如都是文本“0102”一边是数字102一边是文本“0102”合并结果就会产生大量NaN。这是属性挂接最容易踩的坑比空间操作本身更隐蔽。6. 数据质检与交付按这张检查表少走弯路6.1 每批数据都跑一遍的检查项交付前我习惯用一张检查表把数据过一遍。表不长但每条都对应一个翻车现场。检查项怎么看通过标准几何类型ogrinfo -al -so点/线/面和需求一致坐标系gdalsrsinfo或.prj已转等积投影有明确单位字段完整性打开属性表逐列看分类、面积字段非空率不低于95%图斑重叠Dissolve前后数量对比无重复计数面积总量与口径一致编码读取QGIS换GBK/UTF-8各开一次中文不出现乱码检查完成后再跑一遍ogr2ogr -makevalid生成最终交付文件。不要直接在原始下载目录里改来改去交付副本单独放一个文件夹避免后续误操作污染原始数据。最后说一个自己吃过亏的细节。有一年处理喀斯特数据集拿到就拖进ArcGIS图上面色块正常“看起来完全没毛病”直接裁剪出图。结果交付前一天核对面积发现所有数值比测绘口径少了近三分之一。原因就是.prj把投影坐标写成了经纬度图形不缩放根本看不出坐标单位的问题。那次以后立了个规矩任何外部SHP第一件事永远是跑ogrinfo -al -so和gdalsrsinfo两行命令花不了一分钟却能把后续所有统计和叠加的危险提前拆掉。希望帮到你。如果手头也碰到“图没问题、数字对不上”的岩溶数据重投影之后再算一次面积通常就是它。本文还有配套的精品资源点击获取