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