簡介全國12052個大型礦產(chǎn)點(diǎn)位矢量SHP數(shù)據(jù)是一套面向GIS專業(yè)人士、礦產(chǎn)規(guī)劃與地質(zhì)研究人員的空間數(shù)據(jù)集可直接用于礦產(chǎn)分布展示、空間查詢與綜合分析。壓縮包共8個文件涵蓋shp、dbf、prj、sbn、sbx、cpg、shp.xml、shx等類型其中shp保存點(diǎn)位幾何信息、dbf記錄礦產(chǎn)屬性字段、prj定義投影坐標(biāo)、shp.xml提供元數(shù)據(jù)整體僅881KB結(jié)構(gòu)緊湊、便于部署可直接在ArcGIS、QGIS等平臺中打開使用。目前已有131人學(xué)習(xí)適合地質(zhì)勘探、環(huán)境影響評估、土地資源管理及城市規(guī)劃等應(yīng)用場景。依托這份數(shù)據(jù)用戶可獲取礦產(chǎn)地名稱、經(jīng)緯度、利用現(xiàn)狀、地質(zhì)工作程度、礦床成因類型、規(guī)模及礦種等關(guān)鍵指標(biāo)免去自行采集與矢量化工作可快速開展區(qū)域礦產(chǎn)格局分析、疊加其他地理要素輔助決策為科研報(bào)告和規(guī)劃方案提供可靠的數(shù)據(jù)底圖。1. 全國12052個大型礦產(chǎn)點(diǎn)位SHP數(shù)據(jù)一份能直接上工程的點(diǎn)位底圖做地質(zhì)或GIS項(xiàng)目最耗時間往往不是算法而是底圖數(shù)據(jù)的整理。這份“全國12052個大型礦產(chǎn)點(diǎn)位矢量SHP數(shù)據(jù)”把全國范圍內(nèi)的大型礦產(chǎn)地按點(diǎn)要素存放每個點(diǎn)帶礦種、規(guī)模、經(jīng)緯度、行政區(qū)等屬性拿到后就能在ArcGIS、QGIS里直接掛接也方便轉(zhuǎn)成KML、GeoJSON或3DTiles。它能解決的問題很直接你不需要再翻幾十個網(wǎng)頁拼數(shù)據(jù)也不用擔(dān)心點(diǎn)位缺失或坐標(biāo)漂移。適合正在做資源評價(jià)、礦權(quán)分析、生態(tài)紅線疊置或野外查證規(guī)劃的地質(zhì)、GIS從業(yè)者。本文我會順著“先體檢、再轉(zhuǎn)換、后分析”的順序把這套數(shù)據(jù)的用法和坑一次說清。2. 拿到SHP后先做數(shù)據(jù)體檢坐標(biāo)系、字段與幾何合法性SHP數(shù)據(jù)不像Excel打開不報(bào)錯不代表能用。我遇到過好幾個項(xiàng)目數(shù)據(jù)加載后點(diǎn)跑到了國外、屬性全是亂碼、空間查詢結(jié)果為空最后追根溯源都是因?yàn)榈谝徊經(jīng)]做體檢。全國12052個大型礦產(chǎn)點(diǎn)位數(shù)據(jù)量大、來源不一建議先花十分鐘做以下三件事查點(diǎn)數(shù)、驗(yàn)證坐標(biāo)系、檢查幾何問題。2.1 用QGIS讀取SHP并快速統(tǒng)計(jì)點(diǎn)位屬性打開QGIS直接把.shp文件拖進(jìn)圖層區(qū)右鍵圖層選擇“打開屬性表”能看到最前面的幾個要素和屬性。但肉眼檢查12052條記錄效率太低建議用Python控制臺做快速統(tǒng)計(jì)。QGIS內(nèi)置的Python控制臺可以直接寫腳本# 在QGIS菜單“處理-處理控制臺”或“擴(kuò)展-微型Python控制臺”中運(yùn)行 layer iface.activeLayer() print(f要素?cái)?shù)量: {layer.featureCount()}) # 輸出所有字段名和類型確認(rèn)有沒有“規(guī)?！薄暗V種”這樣的分類字段 for f in layer.fields(): print(f.name(), f.typeName()) # 取前3條要素看屬性值是否完整 for feat in layer.getFeatures(): print(feat.attributes()) if feat.id() 2: break這段代碼通過iface.activeLayer()獲取當(dāng)前選中的圖層featureCount()返回12152還是12052一目了然。字段名能幫你判斷數(shù)據(jù)里有沒有“煤礦”“鐵礦”“規(guī)模等級”“坐標(biāo)X”“坐標(biāo)Y”這類列。為什么我要先做這一步因?yàn)楹竺婧Y選、符號化都要依賴字段名如果字段名是拼音縮寫先用紙筆記下來省得分析時找不到字段。參數(shù)說明getFeatures()默認(rèn)遍歷所有要素這里只取前三條避免打印刷屏。如果輸出的要素?cái)?shù)量不等于標(biāo)題里的12052不要急著懷疑數(shù)據(jù)可能下拉框選了子集或隱藏了要素關(guān)掉過濾器再試一次。2.2 坐標(biāo)系識別與投影轉(zhuǎn)換參數(shù)SHP的坐標(biāo)系信息存在同名的.prj文件里。用命令行工具gdalinfo可以快速看到坐標(biāo)系統(tǒng)定義# 查看SHP的完整元數(shù)據(jù)重點(diǎn)看Coordinate System段 gdalinfo 全國大型礦產(chǎn)點(diǎn)位.shp # 只輸出摘要包含要素邊界和投影 ogrinfo -al -so 全國大型礦產(chǎn)點(diǎn)位.shpgdalinfo會打印PROJCS或GEOGCS信息例如GEOGCS[WGS 84, DATUM[WGS_1984]]這說明它是經(jīng)緯度坐標(biāo)EPSG:4326。如果顯示的是PROJCS[CGCS2000 / 3-degree Gauss-Kruger zoned ...]則代表已經(jīng)做過投影坐標(biāo)值可能是米而不是度。注意很多地質(zhì)點(diǎn)位的原始坐標(biāo)是1980西安坐標(biāo)系或CGCS2000如果.prj缺失或?qū)戝e加載后點(diǎn)位會偏移幾百米到幾公里。不確定時先不要動源文件在QGIS圖層面板上右鍵“圖層屬性-源”查看“參考系”。如果顯示“未知”或“用戶定義”別繼續(xù)做分析否則后面所有空間操作都是錯的。我一般的做法是調(diào)出點(diǎn)位的經(jīng)緯度字段在在線地圖上比對幾個已知礦點(diǎn)確認(rèn)坐標(biāo)落在真實(shí)位置上。QGIS中可以直接用“縮放至圖層”看點(diǎn)位輪廓是不是與中國地圖重合。如果確認(rèn)需要將WGS84經(jīng)緯度轉(zhuǎn)成CGCS2000投影坐標(biāo)可以右鍵圖層選擇“導(dǎo)出-另存為”在CRS選擇框里填EPSG:4490CGCS2000地理坐標(biāo)或EPSG:4500系列投影帶。關(guān)鍵是在“打開要素”前勾選“開啟CRS變換”否則只是給坐標(biāo)打上錯誤的標(biāo)簽。2.3 幾何校驗(yàn)自相交、空幾何、重復(fù)點(diǎn)點(diǎn)位數(shù)據(jù)的幾何錯誤通常不顯眼但會影響空間連接和緩沖區(qū)計(jì)算。先用Python的geopands做一個快速體檢import geopandas as gpd gdf gpd.read_file(全國大型礦產(chǎn)點(diǎn)位.shp, encodingutf-8) # 檢查無效幾何和空幾何 print(無效幾何數(shù)量:, (~gdf.is_valid).sum()) print(空幾何數(shù)量:, gdf.is_empty.sum()) # 假設(shè)有經(jīng)度lon、緯度lat字段按坐標(biāo)和礦種名稱一起判斷重復(fù) duplicate gdf[gdf.duplicated(subset[lon, lat])] print(完全重復(fù)坐標(biāo)數(shù)量:, len(duplicate))is_valid在GeoPandas里基于OGR規(guī)則判斷拓?fù)浜戏ㄐ渣c(diǎn)要素一般是有效的但MultiPoint類型偶爾會出現(xiàn)坐標(biāo)重復(fù)。這里假設(shè)數(shù)據(jù)里有l(wèi)on和lat字段如果實(shí)際字段名不同把subset里的列名改掉即可。重復(fù)點(diǎn)可能是同礦多井位也可能是錄入錯誤需要結(jié)合業(yè)務(wù)判斷如果同一個經(jīng)緯度出現(xiàn)多條記錄且礦種相同大概率是重復(fù)如果礦種不同可能是伴生礦不應(yīng)刪除。如果檢查出大量錯誤若不是關(guān)鍵步驟就直接過濾掉不要在臟數(shù)據(jù)上浪費(fèi)時間。需要修復(fù)時可以用熱詞里提到的shapechk工具打開SHP后執(zhí)行Check它會標(biāo)出問題要素Repair會重建幾何。但強(qiáng)烈建議修復(fù)前備份因?yàn)樾迯?fù)算法可能改動點(diǎn)的位置。3. 把礦產(chǎn)點(diǎn)位SHP轉(zhuǎn)成業(yè)務(wù)數(shù)據(jù)KML、GeoJSON與Excel點(diǎn)轉(zhuǎn)SHP拿到SHP后最常見的需求是“把它轉(zhuǎn)成我的工作流能用的格式”。手機(jī)上用地圖軟件看就轉(zhuǎn)KML發(fā)網(wǎng)頁就轉(zhuǎn)GeoJSON原始數(shù)據(jù)是表格就先用Excel轉(zhuǎn)SHP。數(shù)據(jù)格式轉(zhuǎn)換看似簡單坐標(biāo)系和編碼問題卻每天都在坑人——范圍不對、中文亂碼、點(diǎn)位偏移都從這里來。3.1 SHP轉(zhuǎn)KML本地轉(zhuǎn)換與參數(shù)陷阱KML默認(rèn)是經(jīng)緯度坐標(biāo)系EPSG:4326如果你的SHP本身就是WGS84經(jīng)緯度直接轉(zhuǎn)換就行如果它是投影坐標(biāo)務(wù)必要指定源坐標(biāo)系否則轉(zhuǎn)換出的KML點(diǎn)位偏移會非常大。用GDAL自帶的ogr2ogr可以一行完成# 假設(shè)源SHP是CGCS2000經(jīng)緯度轉(zhuǎn)成KML ogr2ogr -f KML 全國大型礦產(chǎn)點(diǎn)位.kml 全國大型礦產(chǎn)點(diǎn)位.shp \ -s_srs EPSG:4490 -t_srs EPSG:4326 \ -dsco NameField礦種參數(shù)-s_srs指源坐標(biāo)系-t_srs指目標(biāo)坐標(biāo)系-dsco NameField礦種表示把“礦種”字段作為KML要素名稱這樣在Google Earth里點(diǎn)開符號直接顯示“鐵礦”“銅礦”等。如果源坐標(biāo)系不對比如實(shí)際是WGS84但你寫了4490每個點(diǎn)會偏移幾十米到幾百米肉眼在衛(wèi)星影像上看得特別清楚。注意KML對UTF-8支持較好如果SHP是GBK編碼建議先轉(zhuǎn)編碼再轉(zhuǎn)格式。轉(zhuǎn)換后用Earth打開KML檢查12052個點(diǎn)是否全部出現(xiàn)。如果出現(xiàn)點(diǎn)特別多或特別少用ogrinfo查看KML要素?cái)?shù)量ogrinfo -al -so 全國大型礦產(chǎn)點(diǎn)位.kml如果要素?cái)?shù)量對不上多半是源SHP里存在無效幾何轉(zhuǎn)換時被跳過。這時回到第2.3節(jié)做幾何修復(fù)別在KML階段硬調(diào)。3.2 Excel表格轉(zhuǎn)SHP從經(jīng)緯度生成點(diǎn)圖層經(jīng)常有人拿來一份礦點(diǎn)調(diào)查表里面只有“礦點(diǎn)名稱”“經(jīng)度”“緯度”“規(guī)?！睕]有SHP。這時候需要把Excel轉(zhuǎn)成點(diǎn)SHP。最可靠的做法是用Python的Pandas加GeoPandasimport pandas as pd import geopandas as gpd from shapely.geometry import Point # 讀取Excel注意字段名不要帶空格 df pd.read_excel(礦點(diǎn)調(diào)查表.xlsx) print(df.head()) # 構(gòu)造幾何列Point(經(jīng)度, 緯度)千萬不能寫成Point(緯度, 經(jīng)度) geometry [Point(x, y) for x, y in zip(df[經(jīng)度], df[緯度])] gdf gpd.GeoDataFrame(df, geometrygeometry) # 設(shè)置坐標(biāo)系為WGS84經(jīng)緯度 gdf.set_crs(epsg4326, inplaceTrue) # 導(dǎo)出SHP編碼使用UTF-8防止屬性表中文亂碼 gdf.to_file(礦點(diǎn)轉(zhuǎn)Shp.shp, encodingutf-8)參數(shù)說明Point(x, y)中x是經(jīng)度y是緯度。如果你習(xí)慣把坐標(biāo)寫成“緯度,經(jīng)度”這里就要改zip(df[緯度], df[經(jīng)度])。很多人在這里翻車所有點(diǎn)都落進(jìn)海里或堆在赤道上就是因?yàn)閄Y寫反了。set_crs(epsg4326)很關(guān)鍵如果不設(shè)置導(dǎo)出SHP會沒有.prj文件其他軟件加載時會問你坐標(biāo)系。另存SHP時字段名會被限制為10個字符比如“礦產(chǎn)資源類型”會變成“礦產(chǎn)資源類”最好提前把字段改成簡短英文導(dǎo)出后再映射回來。3.3 用屬性篩選和空間連接裁切出目標(biāo)區(qū)域點(diǎn)位實(shí)際業(yè)務(wù)里很少用全國全量數(shù)據(jù)。比如做塔里木河流域小流域分析只需要落在流域邊界內(nèi)的礦點(diǎn)。這時候不要手動一條條選用空間連接最安全。先準(zhǔn)備好兩個數(shù)據(jù)全國礦產(chǎn)點(diǎn)位SHP和一個流域邊界SHP用GeoPandas做import geopandas as gpd # 讀取全國礦點(diǎn)數(shù)據(jù)務(wù)必指定編碼 mineral gpd.read_file(全國大型礦產(chǎn)點(diǎn)位.shp, encodingutf-8) # 讀取研究區(qū)邊界并投影到和礦點(diǎn)一致的坐標(biāo) region gpd.read_file(塔里木河流域邊界.shp) print(礦點(diǎn)坐標(biāo)系:, mineral.crs) print(區(qū)域坐標(biāo)系:, region.crs) # 如果不同統(tǒng)一區(qū)域坐標(biāo)系到礦點(diǎn)坐標(biāo)系 if mineral.crs ! region.crs: region region.to_crs(mineral.crs) # 空間連接提取位于區(qū)域內(nèi)或與區(qū)域相交的點(diǎn) joined gpd.sjoin(mineral, region, predicatewithin) joined.to_file(塔里木礦點(diǎn).shp, encodingutf-8) print(區(qū)域內(nèi)礦點(diǎn)數(shù)量:, len(joined))gpd.sjoin默認(rèn)是左連接會保留左側(cè)表所有要素。使用predicatewithin表示點(diǎn)在邊界內(nèi)如果區(qū)域邊界本身有縫隙部分點(diǎn)被遺漏可以改用intersects。crs不一致時程序會報(bào)錯只有做了to_crs統(tǒng)一才能得到正確結(jié)果。輸出SHP時留意坐標(biāo)系是否保留to_file會沿用當(dāng)前GeoDataFrame的crs如果之前沒有設(shè)置導(dǎo)出的文件沒有.prj后續(xù)工具依然會坐標(biāo)錯亂。4. 從點(diǎn)到面熱力圖、緩沖區(qū)與距離分析的工程應(yīng)用12052個點(diǎn)如果只是疊加到地圖上看密密麻麻根本看不出規(guī)律。工程上要把點(diǎn)變成“面”或“距離”信息才有決策價(jià)值。常見的三種做法點(diǎn)位密度熱力圖、固定半徑緩沖區(qū)、計(jì)算與斷層或河流的距離。這些分析都能沉淀成新的矢量結(jié)果供環(huán)評、規(guī)劃或風(fēng)險(xiǎn)評估使用。4.1 制作點(diǎn)位密度熱力圖核密度的參數(shù)選擇熱力圖核密度的目的是發(fā)現(xiàn)礦點(diǎn)聚集區(qū)。QGIS的“熱力圖Kernel Density Estimation”工具很容易用但有兩個參數(shù)要好好選半徑和權(quán)重。如果源數(shù)據(jù)是地理坐標(biāo)度半徑單位也是度1度大約111公里這會導(dǎo)致結(jié)果粗得沒法看。正確的做法是先把數(shù)據(jù)投影到以米為單位的坐標(biāo)系例如CGCS2000投影帶或Web墨卡托再計(jì)算。用QGIS的處理工具箱可以選擇圖層“全國大型礦產(chǎn)點(diǎn)位.shp”勾選“熱力圖”設(shè)置半徑比如50000米、像素大小1000米、衰減系數(shù)0.1。如果想在Python里重復(fù)調(diào)參可以使用QGIS的Processing接口import processing from qgis.core import QgsVectorLayer layer QgsVectorLayer(全國大型礦產(chǎn)點(diǎn)位.shp, points, ogr) # 注意圖層必須已經(jīng)是投影坐標(biāo)系單位米否則半徑單位無效 params { INPUT: layer, RADIUS: 50000, # 50公里 RADIUS_UNITS: 2, # 2代表米 實(shí)際QGIS版本里可接受0像素1米此處需按版本調(diào)整 DECAY: 0.1, # 指數(shù)衰減越小衰減越快 OUTPUT: 熱度.tif } result processing.run(qgis:heatmapkerneldensityestimation, params) print(result[OUTPUT])RADIUS選取沒有唯一正確值半徑偏小熱點(diǎn)碎成一片半徑偏大看不出局部聚集。我一般先運(yùn)行50公里和20公里兩組結(jié)果再把點(diǎn)位疊加到熱力圖上看哪個更符合礦集區(qū)形態(tài)。DECAY控制距離衰減速度0.1是常用值意味著距離中心越遠(yuǎn)影響力下降越快。如果你跑出來的熱力圖全部黑一塊大概率是圖層坐標(biāo)仍然為經(jīng)緯度半徑被解釋為度必須重新投影。4.2 緩沖區(qū)分析評估礦點(diǎn)對周邊環(huán)境的壓力緩沖區(qū)分析最簡單也最容易出錯錯就錯在單位上。GeoPandas中buffer的寬度單位與數(shù)據(jù)坐標(biāo)系一致如果數(shù)據(jù)是WGS84經(jīng)緯度buffer(5000)代表0.005度約500米而不是5公里。所以必須先把數(shù)據(jù)投影到米制坐標(biāo)。import geopandas as gpd mineral gpd.read_file(全國大型礦產(chǎn)點(diǎn)位.shp, encodingutf-8) # 統(tǒng)一投影到Web墨卡托EPSG:3857單位是米 mineral_proj mineral.to_crs(EPSG:3857) # 生成5km緩沖區(qū) buffer_gdf mineral_proj.copy() buffer_gdf[geometry] mineral_proj.geometry.buffer(5000) # 導(dǎo)出緩沖區(qū)圖層并保留礦點(diǎn)屬性 buffer_gdf.to_file(礦點(diǎn)5km緩沖區(qū).shp, encodingutf-8)copy()避免修改原文件。buffer(5000)的5000是米因?yàn)镋PSG:3857下坐標(biāo)單位是米。注意EPSG:3857在高緯度地區(qū)有面積變形如果做全國尺度的科學(xué)分析建議改用分省或分帶的CGCS2000投影坐標(biāo)而不是全局Web墨卡托。我習(xí)慣在每個省級項(xiàng)目中使用CGCS2000 / 3-degree Gauss-Kruger zone能最大限度減少緩沖面積誤差。生成緩沖區(qū)后可以疊加自然保護(hù)區(qū)和生態(tài)紅線用gpd.sjoin判斷重疊情況篩出威脅范圍。注意緩沖區(qū)是面要素而原SHP是點(diǎn)要素兩者可以用空間連接但字段會重復(fù)導(dǎo)出前先刪除不需要的列。4.3 計(jì)算與斷層或河流的最近距離礦點(diǎn)選址安全評估里有一個剛需每個礦點(diǎn)到最近斷層的距離。GeoPandas沒有直接的“最近距離”函數(shù)可以用逐點(diǎn)計(jì)算import geopandas as gpd faults gpd.read_file(斷層線.shp) pts gpd.read_file(全國大型礦產(chǎn)點(diǎn)位.shp, encodingutf-8) # 確保兩個圖層投影到同一個米制坐標(biāo)系否則算出的距離單位是度 faults faults.to_crs(EPSG:3857) pts pts.to_crs(EPSG:3857) # 逐個點(diǎn)計(jì)算到斷層線的最小距離單位為米 distances pts.geometry.apply(lambda p: faults.geometry.distance(p).min()) pts[斷層距離_km] (distances / 1000).round(3) pts.to_crs(epsg4326).to_file(礦點(diǎn)斷層距離.shp, encodingutf-8)faults.geometry.distance(p)會計(jì)算該點(diǎn)上每個斷層段落的距離min取最小值。如果斷層數(shù)據(jù)里包含多個要素這個操作會遍歷所有斷層線12052個點(diǎn)處理起來也很快。注意最后to_crs(epsg4326)是為了讓輸出回到常用坐標(biāo)系方便疊加底圖。距離單位是米除以1000得到公里。我遇到過有人直接輸出“距離”字段卻忘了換算單位結(jié)果讀圖的人誤把500米看成500公里所以在字段名里注明單位“_km”是必要的。5. 避坑指南12052個點(diǎn)常見的5個坑與排查方法光有理論步驟還不夠?qū)嶋H操作中翻車往往集中在幾個細(xì)節(jié)。我在多個項(xiàng)目里反復(fù)踩過這些坑現(xiàn)在把它們匯總成現(xiàn)象、原因、解決三步你遇到類似情況可以直接對號入座。5.1 現(xiàn)象點(diǎn)全部跑到海里、坐標(biāo)超出中國范圍原因坐標(biāo)系被混淆。最常見的是把CGCS2000或西安80數(shù)據(jù)當(dāng)成WGS84加載或者把GCJ-02火星坐標(biāo)錯標(biāo)為WGS84點(diǎn)位整體偏移幾十米到幾百米。另一種情況是投影坐標(biāo)的中央經(jīng)線設(shè)置錯誤導(dǎo)致同一坐標(biāo)在不同投影帶下位移幾百公里。解決立即停止使用這個圖層。先查看.prj文件或QGIS源信息里的CRS然后用一個已知礦點(diǎn)的經(jīng)緯度在在線地圖上核對。如果偏移固定可以檢查原始Excel里的坐標(biāo)來源如果是GPS測的通常為WGS84如果是國土部門提供的多半是CGCS2000。確定正確坐標(biāo)系后右鍵圖層“導(dǎo)出-另存為”正確設(shè)置CRS再做任何后續(xù)分析。5.2 現(xiàn)象屬性表中文全部亂碼或變成問號和菱形原因SHP的屬性表DBF編碼不匹配。國內(nèi)很多SHP基于GBK或GB2312編碼而QGIS和GeoPandas默認(rèn)讀取UTF-8所以中文顯示亂碼。反之如果SHP是UTF-8某些老舊GIS軟件按GBK讀也會亂。解決在QGIS中打開SHP時點(diǎn)擊“處理-選項(xiàng)”在“數(shù)據(jù)源編碼”里手動選擇GBK或UTF-8直到屬性表恢復(fù)正常。命令行里可以用ogr2ogr把編碼轉(zhuǎn)成UTF-8ogr2ogr -f ESRI Shapefile 編碼轉(zhuǎn)換后.shp 原始亂碼.shp \ -lco ENCODINGUTF-8需要注意的是-lco ENCODINGUTF-8只寫入編碼聲明并不會真正改變DBF內(nèi)部字節(jié)所以源文件如果本身就是GBK轉(zhuǎn)換后屬性內(nèi)容不變只是聲明變了QGIS再用UTF-8讀就正常。如果亂碼已經(jīng)嵌進(jìn)幾何字段只能從源數(shù)據(jù)重新編碼轉(zhuǎn)換不要指望圖層修復(fù)功能。5.3 現(xiàn)象圖層能打開但符號化失敗或要素?cái)?shù)量不對原因幾何類型不是標(biāo)準(zhǔn)點(diǎn)可能是MultiPoint或帶Z值也可能存在空幾何或重復(fù)地理坐標(biāo)。QGIS在某些情況下會對無效幾何自動過濾導(dǎo)致顯示數(shù)量減少。12052個數(shù)據(jù)里混入幾十個異常點(diǎn)多見于人工錄入坐標(biāo)不全。解決先用第2.3節(jié)的腳本檢查is_empty、is_valid、重復(fù)坐標(biāo)。如果存在空幾何用以下代碼過濾并導(dǎo)出干凈圖層import geopandas as gpd gdf gpd.read_file(全國大型礦產(chǎn)點(diǎn)位.shp, encodingutf-8) clean gdf[~gdf.geometry.is_empty gdf.geometry.is_valid] clean.to_file(清潔礦點(diǎn).shp, encodingutf-8)如果必須在原圖層基礎(chǔ)上修復(fù)幾何可以使用QGIS的“修復(fù)幾何”算法Processing工具它會根據(jù)拓?fù)湟?guī)則重建點(diǎn)要素。修復(fù)前先備份修復(fù)后重新統(tǒng)計(jì)要素?cái)?shù)量確保仍是12052個點(diǎn)。5.4 現(xiàn)象緩沖區(qū)或空間連接結(jié)果為空兩點(diǎn)明明離得很近卻匹配不上原因兩個圖層的坐標(biāo)系不一致。比如礦點(diǎn)是CGCS2000河流邊界是WGS84在QGIS視圖上因?yàn)椤皩?shí)時CRS變換”顯示重疊但做空間分析時如果未統(tǒng)一CRS系統(tǒng)會強(qiáng)行計(jì)算相交結(jié)果為空。解決在分析前用下面的強(qiáng)制統(tǒng)一坐標(biāo)系region region.to_crs(mineral.crs)或者反過來把所有數(shù)據(jù)都投影到同一EPSG代碼。其次檢查邊界幾何是否存在空洞或分段如果流域邊界是線而不是面需要先用“線轉(zhuǎn)面”工具閉合。還有一點(diǎn)容易被忽略空間連接條件within要求點(diǎn)嚴(yán)格在內(nèi)部如果點(diǎn)正好落在邊界線上換成intersects能多匹配一部分。5.5 現(xiàn)象字段被截?cái)唷?shù)值變成科學(xué)計(jì)數(shù)法原因SHP的DBF格式是上世紀(jì)老標(biāo)準(zhǔn)字段名最多10個字符數(shù)值字段精度有限。當(dāng)你的Excel字段名超過10個字符轉(zhuǎn)換為SHP時會被偷偷截?cái)啾热纭暗V床規(guī)模級別”變成“礦床規(guī)模級”。當(dāng)坐標(biāo)或規(guī)模數(shù)值很大DBF可能以科學(xué)計(jì)數(shù)法存儲讀出來變成“4.56789E08”。解決轉(zhuǎn)換前先把字段名改成不超過10字符的英文字母或拼音縮寫并建立字段映射表。導(dǎo)出后立即用QGIS屬性表檢查數(shù)值字段類型如果字段類型不是Double而變成整數(shù)需要先在Excel里把小數(shù)位保留夠。對12052個點(diǎn)來說坐標(biāo)字段尤其重要建議在源Excel中把經(jīng)緯度格式設(shè)置成普通數(shù)值不要用文本格式否則轉(zhuǎn)SHP后運(yùn)算精度會丟失。6. 進(jìn)階用Python批量處理這12052個礦點(diǎn)從屬性提取到分級出圖當(dāng)我需要快速摸清這批礦點(diǎn)的規(guī)律時會直接用一段短代碼把屬性統(tǒng)計(jì)和分級符號圖一次跑出來。這個套路比在ArcGIS里一步步點(diǎn)菜單高效得多而且結(jié)果可復(fù)現(xiàn)。import geopandas as gpd import matplotlib.pyplot as plt import warnings warnings.filterwarnings(ignore) # 讀入SHP編碼按原始數(shù)據(jù)情況調(diào)整 df gpd.read_file(全國大型礦產(chǎn)點(diǎn)位.shp, encodingutf-8) # 1. 統(tǒng)計(jì)礦種分布輸出前10類 top_types df[礦種].value_counts().head(10) print(top_types) top_types.to_csv(礦種統(tǒng)計(jì).csv) # 2. 按礦種或規(guī)模字段分級繪制點(diǎn)圖 fig, ax plt.subplots(figsize(12, 10)) df.plot(axax, column礦種, categoricalTrue, legendTrue, markersize4) ax.set_axis_off() plt.tight_layout() plt.savefig(礦種分級圖.png, dpi300) # 3. 提取關(guān)鍵屬性輸出成精簡版表格 df[[礦名稱, 礦種, 規(guī)模, lon, lat]].to_excel(礦點(diǎn)簡表.xlsx, indexFalse)value_counts()統(tǒng)計(jì)各礦種數(shù)量to_csv方便后續(xù)處理。df.plot(axax, column礦種, categoricalTrue)按礦種字段分類符號化圖例自動生成。如果字段“礦種”有大量缺失會被稱為NaN并被當(dāng)作一類可以先df df.dropna(subset[礦種])過濾。導(dǎo)出Excel前建議把經(jīng)緯度保留到6位小數(shù)避免失真。驗(yàn)證成果我有一招隨機(jī)抽出20個點(diǎn)位用to_excel輸出再對照在線地圖或野外檢查點(diǎn)看坐標(biāo)有沒有系統(tǒng)偏移。這步雖簡單但能發(fā)現(xiàn)坐標(biāo)系錯誤。有一回我跳過驗(yàn)證直接做緩沖結(jié)果成果發(fā)出去后被人指出所有點(diǎn)偏了3公里因?yàn)樵碨HP的坐標(biāo)實(shí)際上不是WGS84而是Xian80。從那以后我拿到任何SHP第一件事永遠(yuǎn)是核對坐標(biāo)系和抽查點(diǎn)位這個習(xí)慣救了我很多次。希望這套流程和避坑清單也能幫你在處理全國12052個大型礦產(chǎn)點(diǎn)位時少走彎路。本文還有配套的精品資源點(diǎn)擊獲取