ARTICLE DETAIL

资讯详情

深耕郑州网站建设与运营推广的一线实战洞察。

大型矿产点位SHP数据全解析:清洗、投影与空间分析实战

大型矿产点位SHP数据全解析:清洗、投影与空间分析实战 简介这份全国12052个大型矿产点位矢量SHP数据是一套面向GIS从业者、地质勘探人员及资源规划决策者的基础地理信息数据包可用于矿产资源分布分析、开发利用现状研判及空间可视化展示。包体共含8个文件以.shp主文件存储空间几何信息.dbf保存矿产地名称、经纬度、利用现状、地质工作程度、矿床成因类型、规模及矿种等属性.prj提供投影坐标.shx、.sbn、.sbx分别用于空间与属性索引加速.cpg和.xml则补充编码及元数据信息整体约881KB轻量易用。目前已有131人学习下载。数据覆盖全国范围点位信息完整可直接导入ArcGIS、QGIS等平台进行点位渲染、缓冲区分析与叠加制图也可为矿产储量估算、环境影响评价、土地利用规划等场景提供基础支撑适合作为教学实验、科研建模与行业应用的前期数据素材。1. 全国12052个大型矿产点位SHP数据一份能直接做空间分析的“家底清单”做区域找矿潜力评价、矿山生态修复范围划定或者给矿业权设置方案配图时最怕的不是分析方法不够而是缺一份能直接落到地图上的基础数据。全国12052个大型矿产点位矢量SHP数据就是把大型矿产地目录从纸质表格或Excel台账转成了带坐标、带属性、可直接叠加分析的点图层一个点对应一处大型矿产地属性表里记录矿种、规模、所在行政区等信息。这份数据解决的核心问题是把矿产家底从“数字”变成“SHP矢量数据”让你能在GIS里做缓冲区分析、叠加分析、核密度估算和各类空间统计。适合矿产地质工程师、GIS数据分析师、做国土空间规划的人也适合刚接触SHP文件、想拿一份真实数据练手的新手。本文按实际落地路径从数据结构讲到清洗、转换和避坑最后给一套最小可行的分析流程。2. 先翻SHP的“底裤”投影、字段与编码三大基础2.1 拿到SHP先别急着加图层先看这三样SHP文件从来不是一个文件而是一组文件的集合。点图层至少包含主文件.shp存几何形状、索引文件.shx、属性表.dbf再加上可选的投影文件.prj、编码文件.cpg。很多人拿到手直接拖进ArcGIS或QGIS发现属性表是乱码、点全部跑海里回头排查才发现是.prj缺失或.cpg没同步。我拿到这份全国12052个大型矿产点位SHP后的第一个动作是打开文件夹看文件清单。有.prj说明坐标系统有显式声明没有.prj就得靠坐标值反推。常见的反推方式很简单看X和Y的数量级。X是6位比如500000上下、Y是7位比如3000000到4000000这是高斯-克吕格投影的投影坐标带号决定中央经线X是4位、Y是3位经度100-130纬度20-50这是经纬度坐标大概率是WGS84或CGCS2000。检查文件清单用一条命令就够了ls -lh 大型矿产点位.*逻辑很简单如果.shp有几百MB但.dbf只有几十KB说明属性字段很少甚至可能只有FID反过来.dbf很大但.shp很小说明几何信息偏少这在点图层里属于正常现象。参数上注意.dbf单文件理论上限是2GB字段名最长10个字节中文编码下1个汉字占2字节所以字段名超过5个汉字就会被截断这是后面乱码和字段丢失的根源之一。2.2 用Python把属性表翻个底朝天五个必查字段SHP能不能用取决于属性表里有没有可用的关键字段。很多矿产地SHP数据为了兼容老GIS软件字段名会存成拼音缩写或数字编号比如KC、GM、XZQ这时就得先搞清楚每个字段对应什么含义。我习惯用GeoPandas快速读一遍结构import geopandas as gpd from pathlib import Path data_path Path(大型矿产点位.shp) gdf gpd.read_file(data_path, encodinggbk) # dbf常见编码是GBK print(点数:, gdf.shape[0]) print(集合类型:, gdf.geom_type.unique()) print(坐标系:, gdf.crs) print(字段列表:, gdf.columns.tolist()) print(gdf.head(30).to_string())这段代码做四件事确认点数是不是12052、确认几何类型是Point还是MultiPoint、确认坐标系统、预览前30行属性。逻辑上说如果点数对不上说明你拿到的有可能是经过裁剪或合并的版本如果出现MultiPoint后续转GeoJSON或做空间连接时都要先用.explode()拆成单点。参数上重点提两个encodinggbk不是万能的老数据还有用GB2312、GB18030的新版数据可能是UTF-8。一次读不对就依次换编码试。再看具体字段统计# 按矿种字段统计摸清数据里有什么、缺什么 if 矿种 in gdf.columns: stats gdf.groupby(矿种).size().sort_values(ascendingFalse) print(矿种数量前20:) print(stats.head(20)) else: print(属性表里没有矿种字段检查拼音缩写字段) # 检查关键字段的空值比例 key_fields [矿种, 规模, 省, 市, 东经, 北纬] for col in key_fields: if col in gdf.columns: null_cnt gdf[col].isna().sum() print(f{col}: 空值 {null_cnt} 条, 占 {null_cnt / len(gdf):.2%})空值检查这一步不能省。一份12052个点的SHP哪怕只有0.5%的点缺矿种字段在按矿种做分组统计时就会少掉60个点图上直接出现空洞。这套流程在任何SHP矢量数据上都通用不只是矿产地数据。2.3 坐标系陷阱WGS84、CGCS2000与投影坐标哪个是对的点位数据的坐标系是整个分析工作的地基但大多数翻车都发生在这里。现实中存在三种情况一是SHP自带的.prj声明的坐标系和点坐标实际用的坐标系不一致二是点坐标实际是西安80或北京54却被标成了WGS84三是坐标本身没问题但前面的处理步骤里无意间做了一次投影变换。判断坐标系是否匹配有一个土办法抽取几个点落在全国范围内看经纬度是否合理。中国陆地经度约73°E到135°E纬度约18°N到53°N。如果点全部在非洲或大西洋不用怀疑坐标被错读了。常见情形是原始数据用高斯投影坐标X约8位带号、Y约7位但加载时被按经纬度解读直接偏出十万八千里。如果确认当前坐标系和实际不匹配用.to_crs()做转换# 假设数据实际是CGCS2000经纬度但被标成了WGS84 gdf_corrected gdf.set_crs(EPSG:4490) # CGCS2000 地理坐标系 # 再转成WGS84经纬度方便叠加在线底图 gdf_wgs84 gdf_corrected.to_crs(EPSG:4326) print(gdf_wgs84.crs) # 若要做距离或面积量算必须转成投影坐标系 gdf_proj gdf_wgs84.to_crs(EPSG:32649) # WGS84 UTM 49N覆盖东经111-120 print(gdf_proj.geometry.area.describe())这里有个行业约定要说明两个地理坐标系之间的转换在大比例尺制图上可能看不出差异CGCS2000与WGS84在多数地区偏差不到1米。但转到投影坐标系后做缓冲区或面积统计不转换直接用经纬度计算的结果是完全错误的长度单位会被当度来处理。EPSG:32649只覆盖东经111°到120°跨带数据要按点数分区域分别转换这属于容易被忽视的细节。3. 把点位数据整理成能直接分析的样子去重、关联与Excel点转SHP3.1 全量体检重复点、缺失值和无效坐标的清理SHP矢量数据从上游流转到手里通常已经过了好几手。常见问题包括目录式台账里同一矿产地因为多次勘探被录了两遍坐标在Excel转SHP时被保留成文本导致读取异常还有部分点落在国界外的争议区域。直接用原始数据做分析结果会和真实情况产生系统性偏差。我的习惯是先做一次全量体检再动手分析。体检脚本如下import geopandas as gpd gdf gpd.read_file(大型矿产点位.shp, encodinggbk) # 1. 剔除关键字段为空的点 before len(gdf) gdf gdf.dropna(subset[矿种, 规模]) print(f剔除空属性后: {before} - {len(gdf)}) # 2. 按空间位置去重经纬度完全相同且矿种相同视为重复 gdf[lon] gdf.geometry.x gdf[lat] gdf.geometry.y dedup gdf.drop_duplicates(subset[lon, lat, 矿种]) print(f空间矿种去重后: {len(gdf)} - {len(dedup)}) # 3. 剔除非法坐标经度不在73-135纬度不在3-54 valid dedup[ (dedup[lon].between(73, 135)) (dedup[lat].between(3, 54)) ] print(f剔除越界坐标后: {len(dedup)} - {len(valid)})逻辑说明第一步剔除空值是防止分组统计时出现“NaN矿种”的脏分组第二步去重按坐标矿种的组合键因为同一矿种在完全相同的坐标上出现两次从矿产地角度就是重复登记第三步过滤越界坐标虽然简单粗暴但能快速发现坐标系错读的问题。参数说明dedup[lon].between(73, 135)之类的范围是按中国国界给的如果你的分析场景包含南海诸岛纬度范围要放宽到北纬0度如果你同时拿到的还有海域矿产地经度范围要放宽到140度。这里的重点是“边界条件取决于业务范围”不是写死一套规则。3.2 空间关联把点位落到行政区划图上的两种做法做矿产地分析几乎逃不开“点落到哪个省、哪个市”这个问题。常见的行业做法是空间连接Spatial Join在ArcGIS里叫“空间连接”工具在QGIS里叫“连接属性按位置”。我用Python时习惯这样写import geopandas as gpd # 读取全国行政区划面图层注意选择包含南海诸岛的版本 admin gpd.read_file(全国行政区划.shp, encodingutf-8) admin admin[[省, 市, geometry]].copy() # 空间连接点落在哪个面就赋哪个面的属性 joined gpd.sjoin(valid, admin, howleft, predicatewithin) print(未匹配到行政区划的点:, joined[省].isna().sum()) # 如果有点落空原因通常是坐标系不一致或点在边界争议区 unmatched joined[joined[省].isna()] if len(unmatched) 0: unmatched.to_file(unmatched_points.shp, encodingutf-8)这段代码背后的逻辑要讲清楚我没有用ArcGIS默认的“相交intersects”而用了“within包含于”是因为intersects会把落在边界线上的点同时匹配给两侧的行政区导致一个点计到两个省里统计总量时会虚增。within更严格但缺点是在边界上的点会被当成未匹配需要单独导出检查。关于行政区划底图做矿产地分析时建议用带海域、带完整国界线的那种省级面图层而不是简化版。很多做GIS的人去搜“中国行政区划shp怎么下载”会看到省界线文件省1和省2两种省1是带南海诸岛、界限为国界级的完整版省2是简化版。点数据叠加分析一律用省1否则南海诸岛的矿产地会被直接丢掉。3.3 Excel台账转SHPArcGIS手动流程与Python脚本很多矿产地数据最初是以Excel表格形式存在的比如一列经度、一列纬度外加矿种、规模、类型等属性。热词里“arxgis excel点转shp”搜得很多这确实是最常用的操作之一。ArcGIS的流程是文件→添加数据→添加XY数据选X字段经度、Y字段纬度坐标系选WGS84或CGCS2000然后右键图层→数据→导出为SHP。这里最大的坑是Excel表头里有合并单元格、首行空行或坐标列中间混入文本ArcGIS会直接把整列读成字符串导出的点全是一个位置。用Python处理同样场景更可控import pandas as pd import geopandas as gpd df pd.read_excel(大型矿产地台账.xlsx, sheet_name矿产地) # 去掉全空行、去掉坐标非数值行 df df.dropna(subset[经度, 纬度]) df[经度] pd.to_numeric(df[经度], errorscoerce) df[纬度] pd.to_numeric(df[纬度], errorscoerce) df df.dropna(subset[经度, 纬度]) # 从XY列构造点几何 gdf gpd.GeoDataFrame( df, geometrygpd.points_from_xy(df[经度], df[纬度]), crsEPSG:4326 ) # 输出前把字段名长度控制在5个中文字符内 gdf gdf.rename(columns{矿山名称: 矿名, 所在省份: 省}) gdf.to_file(矿产地_from_excel.shp, encodingutf-8)这里有三个参数细节一是errorscoerce会把文本型数字转成NaN而不是报错配合后面的dropna能自动过滤脏行二是points_from_xy的参数顺序是先经度后纬度很多人写反成先纬度后经度结果点全跑到海里三是导出SHP时字段名有长度限制中文超过5个字会被截断或乱码所以要提前改短。如果你最终要交给MapGIS用户还要注意SHP转MapGIS线文件时属性字段会丢失后面避坑章再细说。4. 按使用场景导出格式shp转kml、转json、转txt与3dtiles4.1 ArcGIS/QGIS里shp转kml参数和投影的讲究KML是Google Earth和各类Android地图App通用的交换格式很多现场勘查人员习惯在手机上打开KML看点位。ArcGIS里对应的工具是“图层转KML”Conversion Tools→KML→Layer To KMLQGIS里则是右键图层→导出→另存为格式选“Keyhole Markup Language [KML]”。参数上最容易踩坑的是坐标系统。KML内部强制使用WGS84经纬度EPSG:4326如果你的SHP是CGCS2000或高斯投影坐标工具会自动做投影变换但不一定做得对。我在ArcGIS里一般先把图层用“投影”工具显式转成WGS84再转KML而不是让转换工具隐式处理。输出设置里还有一个“Elevation高度”参数矿产地点位没有高度值要设为0或勾选“Clamp to ground”贴地否则在手机上看点位会浮在空中或扎到地下。在行业应用里“arcgis shp转kml”大多是给野外核查用的所以不需要带属性全字段转KML之前先在属性表里用“删除字段”把没用的列清掉KML文件体积会小很多在手机上打开更快。4.2 命令行方案shp转GeoJSON、shp转txt一把梭在线转换网站能办的事本地命令行都能办而且可控性高得多。GeoJSON是目前WebGIS和前端可视化的事实标准把SHP转成JSON格式后可以直接用于Leaflet、Mapbox GL、Cesium等前端框架。GDAL里的ogr2ogr是行业标准转换工具# SHP转GeoJSON同时重投影到WGS84 ogr2ogr -f GeoJSON 大型矿产点位.geojson 大型矿产点位.shp \ -t_srs EPSG:4326 -lco COORDINATE_PRECISION6 # SHP转带坐标的CSV/TXT制表符分隔 ogr2ogr -f CSV 大型矿产点位.txt 大型矿产点位.shp \ -lco GEOMETRYAS_XY -lco SEPARATORTAB运行逻辑说明第一个命令把字段和几何都带进GeoJSONCOORDINATE_PRECISION6表示坐标保留6位小数约0.1米的精度能显著缩小文件体积第二个命令把几何转成X和Y两列配合原本的属性字段输出到TXT适合给数据统计软件或不在GIS环境里的同事。这里要说一个“json转shp网站”最常见的坑在线转换服务通常支持换格式但属性编码经常丢失中文全部乱码而且大文件传上去不是超时就是被清空。12052个点的量级不大但属性字段多时在线转换一样会翻车。本地用ogr2ogr不存在这些问题。4.3 shp转3dtiles什么时候值得做什么时候别碰“shp转3dtiles”最近搜得很多因为Cesium在Web端做三维矿山可视化很火。但要明确一点3dtiles是为倾斜摄影、BIM、大规模建筑白模设计的把一万多个点转成3dtiles视觉效果上并不会比GeoJSON图标更好反而要多维护一套切片缓存。我的建议是分场景。如果只是做全国点位分布展示直接把SHP转成GeoJSON丢给前端就够了如果是做矿山三维场景中的标识点那也需要结合倾斜摄影模型一起切片而不是单独转。真正值得转3dtiles的情况是数据量达到百万级点位浏览器扛不住DOM节点才需要做点云切片或3dtiles的批量加载。如果确实要转常见做法是先用上面提到的ogr2ogr把SHP转成GeoJSON再用Cesium ion或开源工具链做格式转换。重点说参数3dtiles切片时要设置几何误差GeometricError点数据一般设1到10米LOD层数不必多两到三层就够。这个方向工具链还不太成熟如果你不是专门做Web三维的不必在这个格式上花太多时间。5. 避坑手册乱码、错位、重复点与格式转换的五个翻车现场5.1 属性表全是乱码现象在ArcGIS或QGIS里打开SHP属性表的中文内容变成“锟斤拷”“绗┞”一类符号矿种、规模字段完全不可读。原因SHP的属性表.dbf文件是用GBK或GB2312编码保存的但GIS软件按UTF-8读取或者反过来。解决在QGIS的“图层属性→数据源→数据源编码”里手动改成“GBK”或“GB18030”在ArcGIS里用“表选项→更改数据源编码”调整。如果拿到的是别人转好给你的新SHP用gpd.read_file(..., encodingutf-8)试试如果还是乱码再换gbk。5.2 点位整体飞到海里或境外现象把SHP叠加到在线底图上点位全部落在非洲西海岸或太平洋中间。原因数据本身是高斯投影坐标X带带号、Y七位加载时被按WGS84经纬度解读反之也有。解决先看.prj文件内容用ogrinfo -so 大型矿产点位.shp 大型矿产点位读元数据没有.prj就按前面2.3节的方式反推坐标范围再用set_crs或assign修复。这属于血泪经验我在早期做这类数据时直接把图层一叠加就开始做缓冲区后来发现全部偏了几公里返工了整整两天。5.3 12052个点里有重复点去重键选错现象按矿种统计出来的数量比预期高或者同一个矿种在图上叠了好几个点。原因原始台账按开采区块登记同一座矿山多个矿权重复落点或者Excel转SHP时坐标精度不同导致同一点被记成相邻两点。解决不要只用坐标去重要按“坐标矿种行政区名称”组合键去重。如果属性表里有“矿山编号”或“开采许可证号”优先按这个业务键去重几何位置去重会误删合法的共伴生矿产点。5.4 dwg转shp后整体偏移MAPGIS线文件属性丢光现象把矿区的CAD设计图转成SHP后和遥感影像叠不上整个图幅偏移几十米到上百米。原因CAD图纸用的是建筑坐标系或西安80坐标系转SHP时直接原样搬过来没有做向CGCS2000/WGS84的换算。解决至少找三个均匀分布的控制点如矿区拐点、已知坐标的井口在ArcGIS里用“空间校正”做仿射变换或者让数据提供方直接给出坐标转换七参数。另外SHP转MapGIS线文件时CAD转SHP生成的辅助线、标注线会混进线图层转MapGIS之前先在SHP里把非矿权边界线剔除干净属性字段也会因为MapGIS对dBase字段类型兼容有限而丢失导出时每个字段都要单独检查。这不是技术玄学是真实会发生的损失。5.5 渔网分割SHP后点位数量对不上现象用ArcGIS“渔网”工具做格网统计生成渔网后与点位做空间连接发现统计总和比原始点数少。原因渔网网格被当成面要素与点做连接时默认使用“相交”落在线上的点被同时计给两个网格或者渔网的坐标系和点位坐标系不一致部分点在投影变换后落在网格外。解决做渔网之前先把点位和渔网都转成同一投影坐标系空间连接时用“包含CompletelyContains”而不是“相交Intersects”作为匹配方式。另外一个反直觉的坑是渔网的“模板范围”会受原始范围边界影响生成后要裁剪掉没有意义的海上网格再统计。6. 进阶用核密度和缓冲区叠加把这份数据用出价值拿到清洗后的12052个矿产地SHP只做改格式和配图太可惜了至少应该跑一次核密度分析把“全国大型矿产分布的热点区域”直观呈现出来。在QGIS里用“热图Heatmap”工具核半径Kernel radius先按100公里设置看整体格局再按30公里看局部聚集渲染模式选对数或分位数能让热点层次更清楚。如果用PythonGeoPandas配合scipy.stats.gaussian_kde就能跑import geopandas as gpd from scipy.stats import gaussian_kde import numpy as np gdf gpd.read_file(大型矿产点位_清洗后.shp, encodingutf-8) coords np.vstack([gdf.geometry.x, gdf.geometry.y]) kde gaussian_kde(coords, bw_method0.05) # 在数据范围内生成格网点计算核密度值 xmin, ymin, xmax, ymax gdf.total_bounds x np.linspace(xmin, xmax, 500) y np.linspace(ymin, ymax, 500) X, Y np.meshgrid(x, y) density kde(np.vstack([X.ravel(), Y.ravel()]))这个做法的关键在于bw_method参数它控制核密度的带宽值越小热点越碎值越大格局越宏观建议先跑0.05再看效果迭代。更贴合行业应用的一步是缓冲区和生态敏感区叠加分析。比如评估大型矿产开发对自然保护区的空间影响范围做法是对每个矿产地做10公里缓冲区再与生态红线或自然保护区面图层叠加计算重叠面积。用GeoPandas一条链子就串完buf gpd.GeoDataFrame( gdf[[矿种]], geometrygdf.geometry.buffer(10000), # 10公里单位取决于投影坐标系 crsgdf.crs ) # 叠加求交并计算面积 overlay gpd.overlay(buf, eco_areas, howintersection) overlay[面积] overlay.geometry.area / 1e6 # 转平方公里 impact overlay.groupby([矿种]).面积.sum().sort_values(ascendingFalse) print(impact.head(10))注意这里的单位陷阱buffer(10000)在EPSG:4326经纬度下表示10000度出来的多边形完全错误。必须先转为投影坐标系再缓冲。这套方法同样适用于流域分析类课题——搜“arcswat做小流域分析要用什么矢量数据”的人最终也会需要把矿产地、排污口这类点源数据作为污染源输入进去做法和上面的缓冲区叠加完全一致。核心习惯是所有空间分析之前统一坐标系所有面积量算之前换投影所有统计结果出来之后先回到原始SHP里抽几个点手工验证一遍。我从入行到现在凡是省掉验证步骤的分析基本都出过偏差。希望帮到你。本文还有配套的精品资源点击获取
返回列表