行政區(qū)shp數(shù)據(jù)處理全流程:從字段清洗到空間掛接)
簡(jiǎn)介這份資源面向地理信息、城市規(guī)劃與人口統(tǒng)計(jì)等方向的學(xué)習(xí)者和從業(yè)者提供2017年中國(guó)縣級(jí)行政區(qū)劃的矢量邊界數(shù)據(jù)采用1:100萬(wàn)比例尺適合大范圍空間分析與制圖實(shí)踐。壓縮包共7個(gè)文件、約46.33MB以shp主文件承載幾何圖形dbf保存縣名與代碼等屬性prj定義CGS_WGS_1984坐標(biāo)系sbn、sbx、shx等索引文件加速檢索xml記錄元數(shù)據(jù)配套結(jié)構(gòu)完整可直接在主流GIS軟件中加載使用。已有615人學(xué)習(xí)下載說(shuō)明其在教學(xué)與項(xiàng)目練習(xí)中具有一定參考價(jià)值。讀者可借助該數(shù)據(jù)完成行政區(qū)劃底圖制作、人口密度與地形氣候等專題疊加分析并理解Shapefile各組成文件的分工與坐標(biāo)系統(tǒng)一要點(diǎn)為后續(xù)空間統(tǒng)計(jì)與規(guī)劃決策打下數(shù)據(jù)基礎(chǔ)。1. 拿到一份 2017 縣級(jí)行政區(qū) shp先別急著往地圖里拖一份標(biāo)注為 2017 年、覆蓋縣級(jí)行政區(qū)劃的 shp 文件落到手里第一件事往往不是打開(kāi) QGIS 看效果而是判斷它到底能不能用、該怎么用?;A(chǔ)地理數(shù)據(jù)庫(kù)里的縣級(jí)行政區(qū)圖層本質(zhì)是一份帶屬性表的矢量面數(shù)據(jù)承載的是行政區(qū)名稱、代碼、幾何邊界三樣?xùn)|西。它最常見(jiàn)的用途是給統(tǒng)計(jì)數(shù)據(jù)做空間掛接、給業(yè)務(wù)系統(tǒng)做底圖、給分析任務(wù)做區(qū)域裁剪。但 2017 這個(gè)年份很關(guān)鍵區(qū)劃調(diào)整年年有撤縣設(shè)區(qū)、合并、更名都會(huì)讓邊界和代碼對(duì)不上今天的口徑。所以真正要問(wèn)的是三件事——坐標(biāo)系是什么、屬性字段能不能對(duì)上你的業(yè)務(wù)主鍵、幾何有沒(méi)有拓?fù)鋯?wèn)題。這篇筆記按我實(shí)際處理這類數(shù)據(jù)的順序走一遍從字段結(jié)構(gòu)、坐標(biāo)轉(zhuǎn)換、入庫(kù)掛接到避坑和驗(yàn)證新手能照著復(fù)現(xiàn)熟手能直接跳到參數(shù)和邊界那幾節(jié)。2. 拆開(kāi) shp 看結(jié)構(gòu)字段、編碼與坐標(biāo)系的三個(gè)確認(rèn)點(diǎn)shp 從來(lái)不是單個(gè)文件而是一組同名文件的集合。拿到壓縮包先解壓用命令行把目錄列清楚比在圖形界面里點(diǎn)開(kāi)更可靠因?yàn)槟隳芰⒖贪l(fā)現(xiàn)缺了 .prj 或 .dbf 這種致命問(wèn)題。2.1 先確認(rèn)文件組是否完整一個(gè)能正常讀取的 shapefile 至少要有 .shp幾何、.shx索引、.dbf屬性三個(gè)文件.prj投影定義和 .cpg編碼聲明屬于強(qiáng)烈建議存在的補(bǔ)充文件。缺 .prj 意味著坐標(biāo)系未知缺 .cpg 意味著中文字段可能亂碼。# 列出解壓目錄下所有同名文件組確認(rèn)完整性 ls -lh ./county_2017/ | awk {print $9, $5} # 期望看到成對(duì)出現(xiàn)的 .shp .shx .dbf .prj .cpg這段命令只是把文件名和大小打出來(lái)重點(diǎn)看兩件事同名文件是否成套、.dbf 大小是否合理。如果 .dbf 只有幾 KB 而 .shp 有幾十 MB說(shuō)明屬性表幾乎是空的這種數(shù)據(jù)只能當(dāng)?shù)讏D用做不了屬性掛接。參數(shù)上沒(méi)什么可調(diào)的-lh只是讓人眼可讀awk取第 9 列文件名和第 5 列大小不同系統(tǒng)ls列序略有差異看不準(zhǔn)就直接ls -l人工核對(duì)。2.2 用 Python 讀出字段清單和坐標(biāo)系圖形界面能看但批量判斷和留痕要靠代碼。用 geopandas 讀進(jìn)來(lái)先打印 CRS 和字段類型這一步?jīng)Q定了后面所有轉(zhuǎn)換動(dòng)作。import geopandas as gpd # 讀取時(shí)顯式指定編碼中文屬性表常見(jiàn) GBK/GB18030 gdf gpd.read_file(./county_2017/county_2017.shp, encodinggb18030) print(要素?cái)?shù)量:, len(gdf)) print(坐標(biāo)系:, gdf.crs) print(字段與類型:) print(gdf.dtypes) print(前兩行屬性:) print(gdf.drop(columnsgeometry).head(2))邏輯說(shuō)明encoding參數(shù)是這類國(guó)內(nèi)基礎(chǔ)地理數(shù)據(jù)最容易翻車的地方默認(rèn) UTF-8 讀 GBK 編碼的 .dbf 會(huì)拋 UnicodeDecodeError 或讀出亂碼字段名。先試gb18030它兼容 GBK 且覆蓋更全。gdf.crs如果打印出None說(shuō)明 .prj 缺失或損壞必須手動(dòng)補(bǔ)坐標(biāo)系不能靠猜。字段類型里要特別關(guān)注行政區(qū)代碼字段——如果它被讀成 int64前導(dǎo)零會(huì)丟比如110101變成110101還好但010101這種就會(huì)變成10101掛接時(shí)全對(duì)不上。2.3 行政區(qū)代碼字段的三種常見(jiàn)命名與清洗不同來(lái)源的縣級(jí)圖層代碼字段可能叫ADCODE、CODE、XZQDM、AREA_CODE名稱字段可能叫NAME、XZQMC、COUNTY。不要假設(shè)先看實(shí)際列名再寫映射。清洗的核心是把代碼統(tǒng)一成 6 位字符串并去掉可能存在的空格。# 按實(shí)際列名調(diào)整這里假設(shè)代碼列叫 ADCODE名稱列叫 NAME code_col ADCODE name_col NAME # 統(tǒng)一轉(zhuǎn) 6 位字符串補(bǔ)前導(dǎo)零去空格 gdf[code_col] gdf[code_col].astype(str).str.strip().str.zfill(6) gdf[name_col] gdf[name_col].astype(str).str.strip() # 檢查代碼長(zhǎng)度分布異常長(zhǎng)度說(shuō)明數(shù)據(jù)有問(wèn)題 print(gdf[code_col].str.len().value_counts()) # 正常應(yīng)幾乎全是 6出現(xiàn) 5 或 7 要逐條排查參數(shù)說(shuō)明zfill(6)是補(bǔ)零的關(guān)鍵astype(str)必須在它之前否則對(duì)數(shù)值型補(bǔ)零無(wú)效。str.strip()處理的是從 Excel 或數(shù)據(jù)庫(kù)導(dǎo)出時(shí)混入的首尾空格這類空格肉眼看不見(jiàn)但會(huì)讓merge靜默失敗——這是血淚經(jīng)驗(yàn)里最常見(jiàn)的一類。長(zhǎng)度分布打印出來(lái)如果出現(xiàn)非 6 位不要急著刪先定位是哪些行往往是某個(gè)特殊區(qū)劃或數(shù)據(jù)錄入錯(cuò)誤。3. 坐標(biāo)系與投影從地理坐標(biāo)到可計(jì)算面積的轉(zhuǎn)換路徑縣級(jí)行政區(qū)數(shù)據(jù)最常見(jiàn)的坐標(biāo)系是 CGCS2000 地理坐標(biāo)單位度或某種高斯-克呂格投影單位米。前者適合存儲(chǔ)和展示后者適合算面積、算距離。選錯(cuò)坐標(biāo)系不會(huì)報(bào)錯(cuò)但會(huì)讓面積統(tǒng)計(jì)差出幾十倍這是最隱蔽的坑。3.1 判斷當(dāng)前是地理坐標(biāo)還是投影坐標(biāo)看gdf.crs的輸出如果unit是degree就是地理坐標(biāo)如果是metre就是投影坐標(biāo)。也可以直接看幾何坐標(biāo)數(shù)值范圍經(jīng)度在 73 到 135 之間、緯度在 3 到 54 之間基本可判定為地理坐標(biāo)。# 打印 CRS 詳細(xì)信息和幾何范圍 print(gdf.crs.to_string() if gdf.crs else CRS 缺失) print(幾何邊界:, gdf.total_bounds) # 地理坐標(biāo)下 total_bounds 大致是 [73, 3, 135, 54] 量級(jí)如果total_bounds出現(xiàn)幾百萬(wàn)量級(jí)的數(shù)值說(shuō)明已經(jīng)是投影坐標(biāo)單位大概率是米。這一步的判斷直接決定下一步要不要做投影轉(zhuǎn)換。注意total_bounds返回的是[minx, miny, maxx, maxy]別把經(jīng)緯順序搞反。3.2 面積計(jì)算必須切到等積投影在地理坐標(biāo)下直接調(diào).area得到的是平方度沒(méi)有任何物理意義。正確做法是先投影到適合中國(guó)區(qū)域的等積坐標(biāo)系再算面積。常見(jiàn)選擇是 Albers 等積投影參數(shù)按中國(guó)范圍設(shè)定。# 定義中國(guó)區(qū)域 Albers 等積投影 albers_crs projaea lat_125 lat_247 lat_00 lon_0105 x_00 y_00 datumWGS84 unitsm no_defs gdf_proj gdf.to_crs(albers_crs) gdf_proj[area_km2] gdf_proj.geometry.area / 1e6 print(gdf_proj[[NAME, area_km2]].sort_values(area_km2, ascendingFalse).head(5))參數(shù)說(shuō)明lat_1和lat_2是標(biāo)準(zhǔn)緯線取 25 和 47 覆蓋中國(guó)大部分緯度帶能壓低面積變形lon_0105是中央經(jīng)線取中國(guó)中部經(jīng)度。這套參數(shù)是業(yè)內(nèi)常用配置不是唯一解但對(duì)全國(guó)尺度統(tǒng)計(jì)足夠。/1e6是把平方米轉(zhuǎn)平方公里。如果只做單省分析可以把標(biāo)準(zhǔn)緯線和中央經(jīng)線收窄到該省范圍面積精度會(huì)更好。注意to_crs會(huì)重算幾何原始gdf不受影響別在轉(zhuǎn)換后的對(duì)象上再調(diào)原始字段做混用。3.3 投影轉(zhuǎn)換后的幾何有效性檢查投影轉(zhuǎn)換本身可能放大原本就存在的幾何問(wèn)題比如自相交、重復(fù)點(diǎn)。轉(zhuǎn)換后必須做一次有效性檢查否則后續(xù)空間疊加會(huì)報(bào)拓?fù)溴e(cuò)誤。# 檢查無(wú)效幾何并嘗試修復(fù) invalid gdf_proj[~gdf_proj.geometry.is_valid] print(無(wú)效幾何數(shù)量:, len(invalid)) if len(invalid) 0: gdf_proj[geometry] gdf_proj.geometry.buffer(0) print(修復(fù)后無(wú)效數(shù)量:, (~gdf_proj.geometry.is_valid).sum())buffer(0)是修復(fù)自相交的常用手段原理是把幾何做一次零距離緩沖多數(shù)自相交會(huì)被規(guī)整。但它不是萬(wàn)能藥對(duì)嚴(yán)重破損的幾何可能產(chǎn)生空幾何或碎片修復(fù)后要再查一次is_valid和is_empty。如果無(wú)效數(shù)量占比很高說(shuō)明源數(shù)據(jù)質(zhì)量差應(yīng)該回到源頭核對(duì)而不是硬修。4. 入庫(kù)與掛接把縣級(jí)行政區(qū)接到業(yè)務(wù)統(tǒng)計(jì)表上數(shù)據(jù)清洗完真正的落地動(dòng)作是入庫(kù)和掛接。這一步的目標(biāo)是讓縣級(jí)行政區(qū)圖層能和你的統(tǒng)計(jì)表按行政區(qū)代碼一一對(duì)上對(duì)不上的要能定位原因而不是靠肉眼找。4.1 導(dǎo)出為 GeoPackage 或 PostGIS 的取舍shp 有字段名長(zhǎng)度限制10 字符和編碼問(wèn)題長(zhǎng)期使用建議轉(zhuǎn)成 GeoPackage 或直接入 PostGIS。GeoPackage 單文件、無(wú)字段名限制、UTF-8 原生適合中小規(guī)模PostGIS 適合多用戶并發(fā)和空間查詢。# 導(dǎo)出 GeoPackage字段名不受 10 字符限制 gdf_proj.to_file(./county_2017.gpkg, layercounty, driverGPKG) # 如需入 PostGIS用 sqlalchemy 引擎 # from sqlalchemy import create_engine # engine create_engine(postgresql://user:passhost:5432/gisdb) # gdf_proj.to_postgis(county_2017, engine, if_existsreplace, indexFalse)參數(shù)說(shuō)明driverGPKG必須顯式指定否則 geopandas 可能按擴(kuò)展名推斷失敗。if_existsreplace會(huì)覆蓋同名表生產(chǎn)環(huán)境慎用建議先append再校驗(yàn)。入 PostGIS 前確認(rèn)目標(biāo)庫(kù)已啟用 PostGIS 擴(kuò)展否則建表會(huì)失敗。字段名在 GeoPackage 里可以保留中文但跨系統(tǒng)交換時(shí)仍建議用英文列名加中文別名減少編碼依賴。4.2 按行政區(qū)代碼做左連接并統(tǒng)計(jì)匹配率掛接的核心是左連接以業(yè)務(wù)統(tǒng)計(jì)表為主表縣級(jí)圖層為右表按代碼連接統(tǒng)計(jì)匹配率。匹配率低于預(yù)期時(shí)問(wèn)題一定出在代碼格式或區(qū)劃變更上。import pandas as pd # 模擬業(yè)務(wù)統(tǒng)計(jì)表代碼為 6 位字符串 stats pd.DataFrame({ ADCODE: [110101, 130102, 310101, 999999], value: [120, 340, 210, 50] }) merged stats.merge( gdf_proj[[ADCODE, NAME, geometry]], onADCODE, howleft ) matched merged[NAME].notna().sum() print(f匹配率: {matched}/{len(merged)}) print(未匹配代碼:, merged.loc[merged[NAME].isna(), ADCODE].tolist())邏輯說(shuō)明howleft保證業(yè)務(wù)表所有行保留未匹配的NAME為 NaN便于定位。匹配率是判斷數(shù)據(jù)可用性的硬指標(biāo)低于 95% 就要停下來(lái)查。未匹配代碼常見(jiàn)原因有三類業(yè)務(wù)表代碼帶前導(dǎo)零丟失、圖層是 2017 口徑而業(yè)務(wù)表用了新代碼、代碼字段混入了空格或全角字符。逐類排查比盲目重試有效。4.3 區(qū)劃變更導(dǎo)致的代碼對(duì)不上怎么處理2017 年的縣級(jí)圖層和今天的業(yè)務(wù)數(shù)據(jù)之間必然存在撤縣設(shè)區(qū)、合并更名的情況。處理方式不是改圖層而是建一張代碼映射表把舊代碼映射到新代碼掛接前先過(guò)一遍映射。# 維護(hù)一張變更映射表舊代碼 - 新代碼 code_map pd.DataFrame({ old_code: [130102, 310101], new_code: [130102, 310101] # 示例占位實(shí)際按真實(shí)變更填寫 }) stats stats.merge(code_map, left_onADCODE, right_onold_code, howleft) stats[ADCODE] stats[new_code].fillna(stats[ADCODE]) stats stats.drop(columns[old_code, new_code])這張映射表要作為獨(dú)立資產(chǎn)維護(hù)每次區(qū)劃調(diào)整后更新而不是散落在腳本里。參數(shù)上注意fillna的順序先用新代碼填充沒(méi)有映射的保留原代碼這樣不會(huì)把未變更的記錄誤改。映射表本身要有來(lái)源記錄否則半年后沒(méi)人說(shuō)得清某條映射的依據(jù)。5. 避坑與排查縣級(jí)行政區(qū) shp 處理中最容易翻車的五件事5.1 中文屬性亂碼字段名讀出來(lái)是問(wèn)號(hào)現(xiàn)象gdf.columns打印出一堆\uXXXX或問(wèn)號(hào)NAME列的值也是亂碼。原因.dbf 的編碼聲明缺失或與 .cpg 不一致geopandas 按默認(rèn) UTF-8 解碼 GBK 內(nèi)容。解決讀取時(shí)顯式傳encodinggb18030如果仍亂碼用dbfread單獨(dú)讀 .dbf 確認(rèn)原始編碼再?zèng)Q定是否用iconv轉(zhuǎn)碼后重寫 .cpg。不要反復(fù)用不同編碼試讀后直接覆蓋原文件先備份。5.2 面積算出來(lái)大得離譜或小得離譜現(xiàn)象某縣面積算出幾十萬(wàn)平方公里。原因在地理坐標(biāo)下直接調(diào).area得到的是平方度或者投影參數(shù)的標(biāo)準(zhǔn)緯線設(shè)錯(cuò)。解決先確認(rèn)gdf.crs的unit地理坐標(biāo)必須先to_crs到等積投影再算面積。投影參數(shù)按分析范圍設(shè)定全國(guó)尺度用 Albers 標(biāo)準(zhǔn)緯線 25/47單省分析收窄參數(shù)。算完拿幾個(gè)已知面積的縣做交叉驗(yàn)證偏差超過(guò) 5% 就回頭查參數(shù)。5.3 merge 之后匹配率莫名很低現(xiàn)象代碼肉眼看著一樣merge后大量 NaN。原因一側(cè)是 int 一側(cè)是 str或一側(cè)帶前導(dǎo)零一側(cè)沒(méi)有或混入了不可見(jiàn)空格。解決連接前兩側(cè)都做astype(str).str.strip().str.zfill(6)并在 merge 前各自打印str.len().value_counts()確認(rèn)長(zhǎng)度一致。這類問(wèn)題不會(huì)報(bào)錯(cuò)只能靠匹配率統(tǒng)計(jì)發(fā)現(xiàn)所以每次掛接都要打印匹配率別省這一步。5.4 投影轉(zhuǎn)換后幾何報(bào)拓?fù)溴e(cuò)誤現(xiàn)象to_crs后做空間疊加拋TopologyException。原因源數(shù)據(jù)本身有自相交或重復(fù)點(diǎn)投影放大了問(wèn)題。解決轉(zhuǎn)換后立即做is_valid檢查對(duì)無(wú)效幾何用buffer(0)修復(fù)修復(fù)后再查is_empty。如果無(wú)效比例超過(guò) 1%不要硬修回到源數(shù)據(jù)核對(duì)硬修出來(lái)的邊界可能已經(jīng)偏離實(shí)際。5.5 用 2017 數(shù)據(jù)直接對(duì)今天的業(yè)務(wù)表現(xiàn)象掛接后某些區(qū)縣永遠(yuǎn)匹配不上或匹配上了但邊界明顯不對(duì)。原因2017 至今的區(qū)劃調(diào)整沒(méi)有處理。解決建代碼映射表掛接前先映射對(duì)邊界敏感的分析明確標(biāo)注數(shù)據(jù)年份不要默認(rèn)它代表現(xiàn)狀。年份不是元數(shù)據(jù)里的一個(gè)裝飾字段它決定了這份數(shù)據(jù)能回答什么問(wèn)題、不能回答什么問(wèn)題。6. 驗(yàn)證與進(jìn)階用抽樣核對(duì)和版本留痕把數(shù)據(jù)用穩(wěn)數(shù)據(jù)能不能信不靠感覺(jué)靠可重復(fù)的驗(yàn)證。我一般會(huì)做兩層核對(duì)一層是抽樣人工核對(duì)挑幾個(gè)熟悉的縣把算出的面積和公開(kāi)的行政區(qū)面積對(duì)比另一層是內(nèi)部一致性核對(duì)檢查代碼唯一性、幾何是否重疊、名稱是否有重復(fù)。# 內(nèi)部一致性三項(xiàng)檢查 print(代碼重復(fù):, gdf_proj[ADCODE].duplicated().sum()) print(名稱為空:, gdf_proj[NAME].isna().sum()) print(幾何為空:, gdf_proj.geometry.is_empty.sum()) # 抽樣核對(duì)取幾個(gè)縣打印面積人工比對(duì) sample gdf_proj.sample(5, random_state42)[[NAME, area_km2]] print(sample)duplicated().sum()應(yīng)為 0出現(xiàn)重復(fù)說(shuō)明同一行政區(qū)有多條記錄掛接時(shí)會(huì)產(chǎn)生笛卡爾積式的行數(shù)膨脹。is_empty檢查能抓出投影修復(fù)后產(chǎn)生的空幾何。抽樣用固定random_state保證可復(fù)現(xiàn)方便下次對(duì)比。進(jìn)階用法上我習(xí)慣把整個(gè)處理流程固化成一條可重跑的腳本鏈讀取 → 編碼清洗 → 代碼標(biāo)準(zhǔn)化 → 投影轉(zhuǎn)換 → 有效性修復(fù) → 導(dǎo)出 → 匹配率報(bào)告。每一步的輸出都落盤中間產(chǎn)物帶日期后綴這樣任何一次結(jié)果異常都能回溯到具體環(huán)節(jié)。參數(shù)不寫死在代碼里而是集中到一個(gè)配置字典投影參數(shù)、編碼、代碼字段名都從配置讀換一批數(shù)據(jù)只改配置不改邏輯。檢查項(xiàng)合格標(biāo)準(zhǔn)不合格時(shí)的動(dòng)作代碼唯一性重復(fù)數(shù)為 0定位重復(fù)記錄確認(rèn)是否多版本混入代碼長(zhǎng)度全部為 6 位逐條排查前導(dǎo)零或錄入錯(cuò)誤幾何有效性無(wú)效占比低于 1%低于閾值修復(fù)超過(guò)則回源核對(duì)掛接匹配率高于 95%查代碼格式與區(qū)劃變更映射面積抽樣偏差與參考值偏差低于 5%復(fù)查投影參數(shù)與單位換算這張表我每次處理新批次數(shù)據(jù)都會(huì)過(guò)一遍五項(xiàng)全綠才敢往業(yè)務(wù)里推。最后說(shuō)個(gè)習(xí)慣任何一份帶年份的行政區(qū)數(shù)據(jù)我都會(huì)在文件名和圖層名里保留年份并在屬性表加一列data_year不是為了好看是為了半年后有人問(wèn)“這個(gè)邊界怎么和現(xiàn)在不一樣”時(shí)能一句話說(shuō)清是數(shù)據(jù)口徑問(wèn)題而不是處理錯(cuò)誤。希望幫到你。本文還有配套的精品資源點(diǎn)擊獲取