密度數(shù)據(jù)集全解析：從Shapefile到柵格處理實戰(zhàn))
簡介珠江三角洲城市群區(qū)域開發(fā)密度數(shù)據(jù)集面向城市規(guī)劃、經(jīng)濟地理與區(qū)域可持續(xù)發(fā)展研究者匯集1998年、2006年、2012年三個時點的地理空間信息可用于揭示近二十年間城市空間擴展、經(jīng)濟集聚與開發(fā)強度的時空演變規(guī)律。壓縮包共44個文件、約10.83MB以Shapefile矢量格式.shp、.shx、.dbf、.prj等和GeoTIFF柵格數(shù)據(jù).tif、.tfw、.ovr為主分層存儲研究區(qū)域邊界、經(jīng)濟密度、發(fā)展緊湊度、發(fā)展強度以及1公里網(wǎng)格市域GDP在ArcGIS、QGIS等常用GIS平臺中可直接打開和橫向比較。已有74人瀏覽學習。借助其中包含的城市行政邊界、河流海岸線、GDP空間分布、建成區(qū)緊湊度、建筑覆蓋率與人口密度等指標研究者能逐項對比不同年份的空間結(jié)構(gòu)變化識別經(jīng)濟熱點與城鄉(xiāng)差異亦可延伸至城市可持續(xù)性評估、交通規(guī)劃、環(huán)境保護及社會經(jīng)濟不平等分析為區(qū)域發(fā)展政策制定提供量化支撐。1. 這份珠三角密度數(shù)據(jù)集為什么值得拆開看拿到一個名為“珠江三角洲城市群區(qū)域開發(fā)密度數(shù)據(jù)集199820062012.rar”的壓縮包第一反應(yīng)不是解壓而是先問它能回答什么問題我給城市地理項目做前期數(shù)據(jù)梳理時常遇到這種混合了矢量邊界和柵格GDP的“研究型數(shù)據(jù)集”。它不像標準測繪產(chǎn)品有統(tǒng)一規(guī)范往往帶有一套自己的命名邏輯比如目錄里的2c_EconomicDensity、2b_DevelopmentCompactness、2a_DevelopmentIntensity以及后面的3_NAGDP_1km。這套數(shù)據(jù)覆蓋廣州、深圳、佛山、東莞、中山、珠海、江門、肇慶、惠州九市跨越三個時間斷面用兩種數(shù)據(jù)模型矢量面與1km柵格表達同一區(qū)域的發(fā)展演變。對規(guī)劃師和計量研究者來說最直接的用途是重現(xiàn)“2006是拐點”之類的空間經(jīng)濟命題對GIS工程師來說它則是一份練習時空數(shù)據(jù)清洗、投影轉(zhuǎn)換和柵格-矢量聯(lián)動的標準樣本。下面從解壓這一步開始按實際分析鏈路把每個文件的作用、讀取方式和常見坑逐層拆開盡量做到拿這份數(shù)據(jù)的人能直接照著操作。2. 先把.rar拆開文件結(jié)構(gòu)、坐標系統(tǒng)與柵格/矢量格式核對拿到壓縮包先不要急著在ArcGIS里拖進圖層。數(shù)據(jù)集使用的.rar壓縮格式在Windows上有WinRAR、7-Zip等常見工具在Linux服務(wù)器上則依賴unrar或7z命令行。解壓前先確認目標目錄不含中文路徑否則后續(xù)用Python處理時容易在路徑解析上踩編碼坑。2.1 解壓不是雙擊了事Linux命令與完整性校驗在Linux下常見的做法是先安裝unrar然后執(zhí)行解壓命令sudo apt install unrar # Debian/Ubuntu unrar x 珠江三角洲城市群區(qū)域開發(fā)密度數(shù)據(jù)集199820062012.rar ./prd_dataset/參數(shù)說明x表示保留壓縮包內(nèi)的目錄結(jié)構(gòu)./prd_dataset/是解壓目標目錄。解壓后建議用unrar t 壓縮包名.rar做一次完整性測試因為后續(xù)數(shù)據(jù)里的.sbn/.sbx如果損壞雖然不致命但會影響ArcGIS的空間索引讀取。在Windows下我更推薦用7-Zip的右鍵“解壓到當前文件夾”它對中文文件名和全角括號的處理比某些版本更穩(wěn)定。解壓完成后先看目錄樹確認是否存在1_studyarea、2a_DevelopmentIntensity等頂層文件夾以及每個文件夾內(nèi)是否都有成對的.shp和.dbf文件。這一步雖然瑣碎卻能在正式分析前暴露文件缺失。2.2 Shapefile的十個零件哪些文件必須一起拷貝Shapefile不是單文件而是多個文件集合。本數(shù)據(jù)集中每個矢量圖層都有至少8個同名文件例如EconomicDensity.shp、.shx、.dbf、.prj、.cpg、.sbn、.sbx、.shp.xml。其中.shp存幾何.shx存索引.dbf存屬性三者缺一不可.prj存投影信息.cpg存屬性表編碼.sbn/.sbx是ArcGIS維護的空間索引.shp.xml是元數(shù)據(jù)。擴展名作用缺失后的影響.shp幾何要素無法讀取.shx幾何索引部分庫會報錯.dbf屬性表沒有任何屬性字段.prj坐標系定義GIS軟件可能不識別坐標.cpg字符集聲明中文字段名或值亂碼.sbn/.sbx空間索引僅影響查詢性能.shp.xml元數(shù)據(jù)可忽略值得注意的是DevelopmentIntensity文件夾內(nèi).cpg被寫成了大寫.CPG。在Linux下文件系統(tǒng)區(qū)分大小寫如果從Windows復制到Linuxgeopandas自動搜索.cpg時可能找不到導致屬性編碼識別失敗。我遇到這種情況時一般會用一個for循環(huán)把所有文件改為小寫擴展名或者直接在讀取時顯式指定encoding。2.3 先看投影.prj與.tfw告訴你在哪種坐標下投影是一切空間計算的基準。用文本編輯器打開StudyArea.prj里面是一段WKT字符串。珠三角區(qū)域常見的坐標系有WGS84地理坐標、WGS84 UTM 50N投影以及國家2000坐標系具體要看文件內(nèi)容。我判斷投影時重點關(guān)注三個信息基準面datum、投影名稱projection、單位units。地理坐標系單位為度投影坐標系單位為米。NAGDP_1km_1998.tfw是GeoTIFF的世界文件用普通文本編輯器打開就能看到六行數(shù)字前兩行是像元尺寸和旋轉(zhuǎn)項后兩行是柵格原點。更直觀的方式是用gdalinfogdalinfo NAGDP_1km_1998.tif | head -n 40gdalinfo輸出里會顯示Coordinate System is、Origin和Pixel Size。如果三個年份片的坐標系不一致不要貿(mào)然做差值運算必須先統(tǒng)一投影。對于這份數(shù)據(jù)1_studyarea提供的邊界可以看作所有后續(xù)操作的空間基準。3. 三種密度指標的定義與屬性表讀取用geopandas把字段讀出來目錄里三個以數(shù)字開頭的子文件夾分別叫2c_EconomicDensity、2b_DevelopmentCompactness、2a_DevelopmentIntensity。這個“2c→2a”的字母順序有講究經(jīng)濟密度是結(jié)果緊湊度是形態(tài)開發(fā)強度是過程入口。摘要里解釋得很清楚經(jīng)濟密度度量單位面積經(jīng)濟產(chǎn)值緊湊度度量建成區(qū)集中程度開發(fā)強度則反映建設(shè)密度和人口密度。三者放在一起才能解釋城市擴張的機理。3.1 三個shapefile的字段差異與命名啟示我用ogrinfo或geopandas快速列出字段總能發(fā)現(xiàn)一些規(guī)律。通常EconomicDensity.dbf持有類似GDP_KM2、VAL1998這樣的字段DevelopmentCompactness.dbf里有類似SHAPE_INDEX的值DevelopmentIntensity.dbf里則可能是BUILD_PCT或POP_KM2。雖然名字不同但數(shù)值都是區(qū)域級統(tǒng)計量。命名前綴“2c、2b、2a”更像是一個生產(chǎn)流水線先算出經(jīng)濟產(chǎn)出密度再評估形態(tài)緊湊度最后落到開發(fā)強度。這樣設(shè)計的好處是屬性表里可以不冗余存儲幾何信息靠FID關(guān)聯(lián)。3.2 geopandas讀取示例import geopandas as gpd # 注意cpg編碼聲明如果亂碼就指定encodinggbk eco_1998 gpd.read_file(./2c_EconomicDensity/EconomicDensity.shp, encodingutf-8) print(eco_1998.columns) print(eco_1998.head()) print(eco_1998[[GDP_KM2, YEAR]].describe())這段代碼先用read_file加載 Shapefileencoding參數(shù)匹配.cpg文件聲明的字符集。然后打印所有列名和統(tǒng)計摘要。describe()會輸出count、mean、std等字段能快速判斷GDP_KM2是否存在異常負值或大范圍空值。拿到屬性表后最好順帶檢查YEAR字段因為一個shp文件里可能存儲多期數(shù)據(jù)靠YEAR區(qū)分。3.3 緊湊度與強度的計算口徑緊湊度計算公式在文獻里常用C 2√(πA)/P其中A是建成區(qū)面積P是周長結(jié)果落在0到1之間。本數(shù)據(jù)集里的DevelopmentCompactness如果直接是數(shù)值大概率就是這類形狀指數(shù)或緊湊比。至于DevelopmentIntensity我見過更合理的字段是“單位面積建設(shè)用地面積”或“人口密度”的網(wǎng)格化統(tǒng)計結(jié)果。比較年份時要注意如果行政邊界本身發(fā)生了調(diào)整比如2006年某鎮(zhèn)并入街道則必須先使用1_studyarea里統(tǒng)一后的邊界做裁剪或交集不然緊湊度會被區(qū)域面積變化干擾。4. 跨年份空間分析如何提取“開發(fā)密度”的變化信號三個年份的數(shù)據(jù)放一起最直接的分析是看同一行政單元在不同指標上的走勢。但前提是邊界一致。1_studyarea提供了研究區(qū)邊界而各年份shp的幾何可能來自不同來源必須先處理掉縫隙、重疊與坐標系差異。4.1 統(tǒng)一坐標系與屬性對齊先檢查三個文件的坐標系是否一致import geopandas as gpd study gpd.read_file(./1_studyarea/StudyArea.shp) eco_1998 gpd.read_file(./2c_EconomicDensity/EconomicDensity.shp, encodingutf-8) print(study.crs) print(eco_1998.crs) if eco_1998.crs ! study.crs: eco_1998 eco_1998.to_crs(study.crs)如果.prj文件缺失或定義不完整to_crs會拋錯。應(yīng)急方案是從一個已知正確的GeoDataFrame復制crs對象例如eco_1998.crs study.crs但這樣做的前提是你確認兩套數(shù)據(jù)的橢球和基準面一致否則會出現(xiàn)幾米的偏移。屬性對齊更麻煩。如果2006年的字段結(jié)構(gòu)變了比如新增了“空間GDP”字段而1998年沒有直接用pd.concat會生成大量NaN。我常用的做法是先只保留幾個核心字段統(tǒng)一命名為region、year、value再按區(qū)域名進行merge或groupby。# 讀取緊湊度數(shù)據(jù)并做空間連接 compact gpd.read_file(./2b_DevelopmentCompactness/DevelopmentCompactness.shp) joined gpd.sjoin(study, compact, howleft, predicateintersects) # 按城市名聚合 grouped joined.groupby(CITY_NAME)[COMPACT].mean().reset_index()代碼中sjoin的predicateintersects表示兩個要素只要任意部分相交就算匹配這能規(guī)避邊界微小縫隙。groupby取均值是因為一個城市可能由多個面要素組成。如果希望面積加權(quán)則需要先計算相交面積joined[inter_area] joined.geometry.intersection(study.unary_union).area weighted joined.groupby(CITY_NAME).apply( lambda df: (df[COMPACT] * df[inter_area]).sum() / df[inter_area].sum() )注意apply里的加權(quán)公式當區(qū)域面積差異大時簡單均值會偏向面積小但數(shù)值高的單元加權(quán)均值能反映整體狀態(tài)。4.2 生成三年對比表把三個年份的經(jīng)濟密度整合到一張寬表可以直接用pivot_tableimport pandas as pd # 假設(shè)三個gdf都有YEAR和VALUE字段 eco_1998[year] 1998 eco_2006[year] 2006 eco_2012[year] 2012 frames [eco_1998[[CITY_NAME, year, GDP_KM2]], eco_2006[[CITY_NAME, year, GDP_KM2]], eco_2012[[CITY_NAME, year, GDP_KM2]]] all_df pd.concat(frames) pivot all_df.pivot_table(indexCITY_NAME, columnsyear, valuesGDP_KM2, aggfuncmean) pivot[growth_06_98] pivot[2006] / pivot[1998] - 1 pivot[growth_12_06] pivot[2012] / pivot[2006] - 1pivot_table的aggfunc默認mean它會處理同一城市下多個面要素的重復記錄。用增長率列能直接識別哪幾年擴張最快。4.3 變化可視化可視化可以先用matplotlib畫折線趨勢再用geopandas給每個城市填充增長率顏色。需要注意用顏色分級圖時最好將增長率分位數(shù)分成5類避免極端值把配色拉平。5. 處理1km GDP柵格用rasterio讀取、重投影與統(tǒng)計柵格部分是這份數(shù)據(jù)里最硬的骨頭。NAGDP_1km_1998.tif這類文件是GeoTIFF除了.tif本體還有.tfw世界文件、.aux.xml輔助信息、.ovr金字塔。這套文件對新手不太友好但處理思路很固定先看元數(shù)據(jù)再處理空值最后做統(tǒng)計或重投影。5.1 柵格文件的組成.tif、.tfw、.ovr各干什么.tfw是文本世界文件記錄像元尺寸和左上角坐標.ovr是金字塔文件如果缺失大范圍縮放變慢但不影響數(shù)值精度.aux.xml包含統(tǒng)計信息和nodata標記。讀取時我習慣用rasterio而不是gdal命令行因為能直接返回numpy數(shù)組import rasterio with rasterio.open(./3_NAGDP_1km/NAGDP_1km_1998.tif) as src: data src.read(1) nodata src.nodata profile src.profile print(nodata , nodata) print(數(shù)據(jù)類型 , data.dtype) print(shape , data.shape)src.read(1)讀出第一波段nodata是無值標記。GDP數(shù)據(jù)通常用-9999或0表示無數(shù)據(jù)拿到數(shù)組后第一步要把這些值屏蔽掉否則計入均值會讓結(jié)果嚴重偏低import numpy as np masked np.ma.masked_where(data 0, data) print(masked.mean())5.2 做年度差值并重投影三個年份的tif如果像元尺寸和原點不完全一致不要直接做數(shù)組減法。正確做法是先統(tǒng)一到同一個格網(wǎng)。可以用reproject完成from rasterio.warp import reproject, Resampling, calculate_default_transform with rasterio.open(./3_NAGDP_1km/NAGDP_1km_1998.tif) as src: data_1998 src.read(1) src_crs src.crs src_transform src.transform nodata src.nodata # 目標坐標系如果源是WGS84地理坐標這里用EPSG:4326 transform, width, height calculate_default_transform( src_crs, EPSG:32650, src.width, src.height, *src.bounds ) dst np.zeros((height, width), dtypenp.float32) reproject( sourcedata_1998, destinationdst, src_transformsrc_transform, src_crssrc_crs, dst_transformtransform, dst_crsEPSG:32650, resamplingResampling.bilinear, src_nodatanodata, dst_nodata-9999 )上述代碼中calculate_default_transform根據(jù)源范圍和目標投影計算新尺寸bilinear重采樣適用于連續(xù)型GDP值如果處理分類柵格則改用nearest。src_nodata與dst_nodata保持一致保證重投影后空洞仍是同一標記。如果不想寫這么長的代碼優(yōu)先用rioxarray它將坐標和投影封裝得更好一行就能做差值import rioxarray ds_1998 rioxarray.open_rasterio(./3_NAGDP_1km/NAGDP_1km_1998.tif) ds_2006 rioxarray.open_rasterio(./3_NAGDP_1km/NAGDP_1km_2006.tif) diff ds_2006 - ds_1998但要先確認兩個對象的空間參考一致若不一致則調(diào)用ds_1998.rio.reproject_match(ds_2006)對齊。5.3 高分辨率GDP數(shù)據(jù)的局限與坑1km格網(wǎng)GDP并不是真實觀測值它一般基于夜間燈光、土地利用和統(tǒng)計年鑒做空間化插值因此“某個像元上的GDP”并不代表那個位置真有這么多產(chǎn)出。解讀時更適合用在區(qū)域總量或相對比上。我在做行政單元統(tǒng)計時會先用行政區(qū)邊界裁剪然后逐像元累加。這里最典型的坑有兩個一是邊緣像元被邊界切成窄條如果連通面積小于一個像元直接累加會高估邊界地帶數(shù)值二是.ovr金字塔文件版本不一致可能讓某些軟件讀到錯誤nodata。解決辦法是用rasterstats的all_touchedFalse只保留完整落入邊界的像元并在讀取時顯式指定maskedTrue。提示如果某一年份的GDP均值明顯異常優(yōu)先檢查nodata值是否被默認為0再檢查重投影時是否把邊緣的nodata像元插值成小數(shù)值。6. 收尾技巧用rasterstats一鍵把1km GDP匯總到行政區(qū)前面的操作分別處理了矢量和柵格實際交付成果時常需要“一張Excel表”每一行是一個區(qū)縣列是三個年份的GDP密度與總量。這時用rasterstats包可以省掉大量手寫循環(huán)。它專門用來把柵格數(shù)據(jù)按矢量多邊形進行分區(qū)統(tǒng)計pip install rasterstats然后寫一個極簡函數(shù)import pandas as pd from rasterstats import zonal_stats stats zonal_stats( ./1_studyarea/StudyArea.shp, # 矢量邊界 ./3_NAGDP_1km/NAGDP_1km_1998.tif, # 柵格文件 stats[sum, mean, min, max], nodata-9999, # 與tif元數(shù)據(jù)一致 all_touchedFalse # 只統(tǒng)計完整覆蓋的像元 ) df pd.DataFrame(stats)zonal_stats第一個參數(shù)接收矢量邊界第二個參數(shù)接收柵格路徑stats支持sum/mean/max/min/median等聚合統(tǒng)計。nodata參數(shù)能過濾柵格中的無值像元避免把-9999當作真實數(shù)值。all_touchedFalse表示只統(tǒng)計多邊形中心點位于內(nèi)部的像元這比True更嚴格也更能避免邊緣噪聲。三個年份只需包裝成一個循環(huán)def extract_year(year): raster f./3_NAGDP_1km/NAGDP_1km_{year}.tif st zonal_stats( ./1_studyarea/StudyArea.shp, raster, stats[sum, mean], nodata-9999, all_touchedFalse ) return pd.DataFrame(st) result pd.concat( [extract_year(1998), extract_year(2006), extract_year(2012)], axis1 ) result.columns [GDP_sum_1998, GDP_mean_1998, GDP_sum_2006, GDP_mean_2006, GDP_sum_2012, GDP_mean_2012]這個封裝直接生成一個DataFrame再配合study里的行政區(qū)名稱列即可保存為CSVstudy gpd.read_file(./1_studyarea/StudyArea.shp) result.insert(0, city_name, study[CITY_NAME]) result.to_csv(./prd_gdp_summary.csv, indexFalse)一個小技巧如果發(fā)現(xiàn)某年的sum異常低先回到第5章檢查nodata和重投影步驟或者換個思路用面積加權(quán)平均——因為柵格像元面積可能隨緯度變化而不同尤其在投影坐標系里不同緯度像元實際面積并不完全一致。遇到這種情況可以先用area屬性求出單位面積值再乘實際面積。這套流程跑通后后面再加2018年數(shù)據(jù)只需要復制這個函數(shù)模板替換年份即可。本文還有配套的精品資源點擊獲取 珠三角開發(fā)密度數(shù)據(jù)集全解析：從Shapefile到柵格處理實戰(zhàn)
簡介珠江三角洲城市群區(qū)域開發(fā)密度數(shù)據(jù)集面向城市規(guī)劃、經(jīng)濟地理與區(qū)域可持續(xù)發(fā)展研究者匯集1998年、2006年、2012年三個時點的地理空間信息可用于揭示近二十年間城市空間擴展、經(jīng)濟集聚與開發(fā)強度的時空演變規(guī)律。壓縮包共44個文件、約10.83MB以Shapefile矢量格式.shp、.shx、.dbf、.prj等和GeoTIFF柵格數(shù)據(jù).tif、.tfw、.ovr為主分層存儲研究區(qū)域邊界、經(jīng)濟密度、發(fā)展緊湊度、發(fā)展強度以及1公里網(wǎng)格市域GDP在ArcGIS、QGIS等常用GIS平臺中可直接打開和橫向比較。已有74人瀏覽學習。借助其中包含的城市行政邊界、河流海岸線、GDP空間分布、建成區(qū)緊湊度、建筑覆蓋率與人口密度等指標研究者能逐項對比不同年份的空間結(jié)構(gòu)變化識別經(jīng)濟熱點與城鄉(xiāng)差異亦可延伸至城市可持續(xù)性評估、交通規(guī)劃、環(huán)境保護及社會經(jīng)濟不平等分析為區(qū)域發(fā)展政策制定提供量化支撐。1. 這份珠三角密度數(shù)據(jù)集為什么值得拆開看拿到一個名為“珠江三角洲城市群區(qū)域開發(fā)密度數(shù)據(jù)集199820062012.rar”的壓縮包第一反應(yīng)不是解壓而是先問它能回答什么問題我給城市地理項目做前期數(shù)據(jù)梳理時常遇到這種混合了矢量邊界和柵格GDP的“研究型數(shù)據(jù)集”。它不像標準測繪產(chǎn)品有統(tǒng)一規(guī)范往往帶有一套自己的命名邏輯比如目錄里的2c_EconomicDensity、2b_DevelopmentCompactness、2a_DevelopmentIntensity以及后面的3_NAGDP_1km。這套數(shù)據(jù)覆蓋廣州、深圳、佛山、東莞、中山、珠海、江門、肇慶、惠州九市跨越三個時間斷面用兩種數(shù)據(jù)模型矢量面與1km柵格表達同一區(qū)域的發(fā)展演變。對規(guī)劃師和計量研究者來說最直接的用途是重現(xiàn)“2006是拐點”之類的空間經(jīng)濟命題對GIS工程師來說它則是一份練習時空數(shù)據(jù)清洗、投影轉(zhuǎn)換和柵格-矢量聯(lián)動的標準樣本。下面從解壓這一步開始按實際分析鏈路把每個文件的作用、讀取方式和常見坑逐層拆開盡量做到拿這份數(shù)據(jù)的人能直接照著操作。2. 先把.rar拆開文件結(jié)構(gòu)、坐標系統(tǒng)與柵格/矢量格式核對拿到壓縮包先不要急著在ArcGIS里拖進圖層。數(shù)據(jù)集使用的.rar壓縮格式在Windows上有WinRAR、7-Zip等常見工具在Linux服務(wù)器上則依賴unrar或7z命令行。解壓前先確認目標目錄不含中文路徑否則后續(xù)用Python處理時容易在路徑解析上踩編碼坑。2.1 解壓不是雙擊了事Linux命令與完整性校驗在Linux下常見的做法是先安裝unrar然后執(zhí)行解壓命令sudo apt install unrar # Debian/Ubuntu unrar x 珠江三角洲城市群區(qū)域開發(fā)密度數(shù)據(jù)集199820062012.rar ./prd_dataset/參數(shù)說明x表示保留壓縮包內(nèi)的目錄結(jié)構(gòu)./prd_dataset/是解壓目標目錄。解壓后建議用unrar t 壓縮包名.rar做一次完整性測試因為后續(xù)數(shù)據(jù)里的.sbn/.sbx如果損壞雖然不致命但會影響ArcGIS的空間索引讀取。在Windows下我更推薦用7-Zip的右鍵“解壓到當前文件夾”它對中文文件名和全角括號的處理比某些版本更穩(wěn)定。解壓完成后先看目錄樹確認是否存在1_studyarea、2a_DevelopmentIntensity等頂層文件夾以及每個文件夾內(nèi)是否都有成對的.shp和.dbf文件。這一步雖然瑣碎卻能在正式分析前暴露文件缺失。2.2 Shapefile的十個零件哪些文件必須一起拷貝Shapefile不是單文件而是多個文件集合。本數(shù)據(jù)集中每個矢量圖層都有至少8個同名文件例如EconomicDensity.shp、.shx、.dbf、.prj、.cpg、.sbn、.sbx、.shp.xml。其中.shp存幾何.shx存索引.dbf存屬性三者缺一不可.prj存投影信息.cpg存屬性表編碼.sbn/.sbx是ArcGIS維護的空間索引.shp.xml是元數(shù)據(jù)。擴展名作用缺失后的影響.shp幾何要素無法讀取.shx幾何索引部分庫會報錯.dbf屬性表沒有任何屬性字段.prj坐標系定義GIS軟件可能不識別坐標.cpg字符集聲明中文字段名或值亂碼.sbn/.sbx空間索引僅影響查詢性能.shp.xml元數(shù)據(jù)可忽略值得注意的是DevelopmentIntensity文件夾內(nèi).cpg被寫成了大寫.CPG。在Linux下文件系統(tǒng)區(qū)分大小寫如果從Windows復制到Linuxgeopandas自動搜索.cpg時可能找不到導致屬性編碼識別失敗。我遇到這種情況時一般會用一個for循環(huán)把所有文件改為小寫擴展名或者直接在讀取時顯式指定encoding。2.3 先看投影.prj與.tfw告訴你在哪種坐標下投影是一切空間計算的基準。用文本編輯器打開StudyArea.prj里面是一段WKT字符串。珠三角區(qū)域常見的坐標系有WGS84地理坐標、WGS84 UTM 50N投影以及國家2000坐標系具體要看文件內(nèi)容。我判斷投影時重點關(guān)注三個信息基準面datum、投影名稱projection、單位units。地理坐標系單位為度投影坐標系單位為米。NAGDP_1km_1998.tfw是GeoTIFF的世界文件用普通文本編輯器打開就能看到六行數(shù)字前兩行是像元尺寸和旋轉(zhuǎn)項后兩行是柵格原點。更直觀的方式是用gdalinfogdalinfo NAGDP_1km_1998.tif | head -n 40gdalinfo輸出里會顯示Coordinate System is、Origin和Pixel Size。如果三個年份片的坐標系不一致不要貿(mào)然做差值運算必須先統(tǒng)一投影。對于這份數(shù)據(jù)1_studyarea提供的邊界可以看作所有后續(xù)操作的空間基準。3. 三種密度指標的定義與屬性表讀取用geopandas把字段讀出來目錄里三個以數(shù)字開頭的子文件夾分別叫2c_EconomicDensity、2b_DevelopmentCompactness、2a_DevelopmentIntensity。這個“2c→2a”的字母順序有講究經(jīng)濟密度是結(jié)果緊湊度是形態(tài)開發(fā)強度是過程入口。摘要里解釋得很清楚經(jīng)濟密度度量單位面積經(jīng)濟產(chǎn)值緊湊度度量建成區(qū)集中程度開發(fā)強度則反映建設(shè)密度和人口密度。三者放在一起才能解釋城市擴張的機理。3.1 三個shapefile的字段差異與命名啟示我用ogrinfo或geopandas快速列出字段總能發(fā)現(xiàn)一些規(guī)律。通常EconomicDensity.dbf持有類似GDP_KM2、VAL1998這樣的字段DevelopmentCompactness.dbf里有類似SHAPE_INDEX的值DevelopmentIntensity.dbf里則可能是BUILD_PCT或POP_KM2。雖然名字不同但數(shù)值都是區(qū)域級統(tǒng)計量。命名前綴“2c、2b、2a”更像是一個生產(chǎn)流水線先算出經(jīng)濟產(chǎn)出密度再評估形態(tài)緊湊度最后落到開發(fā)強度。這樣設(shè)計的好處是屬性表里可以不冗余存儲幾何信息靠FID關(guān)聯(lián)。3.2 geopandas讀取示例import geopandas as gpd # 注意cpg編碼聲明如果亂碼就指定encodinggbk eco_1998 gpd.read_file(./2c_EconomicDensity/EconomicDensity.shp, encodingutf-8) print(eco_1998.columns) print(eco_1998.head()) print(eco_1998[[GDP_KM2, YEAR]].describe())這段代碼先用read_file加載 Shapefileencoding參數(shù)匹配.cpg文件聲明的字符集。然后打印所有列名和統(tǒng)計摘要。describe()會輸出count、mean、std等字段能快速判斷GDP_KM2是否存在異常負值或大范圍空值。拿到屬性表后最好順帶檢查YEAR字段因為一個shp文件里可能存儲多期數(shù)據(jù)靠YEAR區(qū)分。3.3 緊湊度與強度的計算口徑緊湊度計算公式在文獻里常用C 2√(πA)/P其中A是建成區(qū)面積P是周長結(jié)果落在0到1之間。本數(shù)據(jù)集里的DevelopmentCompactness如果直接是數(shù)值大概率就是這類形狀指數(shù)或緊湊比。至于DevelopmentIntensity我見過更合理的字段是“單位面積建設(shè)用地面積”或“人口密度”的網(wǎng)格化統(tǒng)計結(jié)果。比較年份時要注意如果行政邊界本身發(fā)生了調(diào)整比如2006年某鎮(zhèn)并入街道則必須先使用1_studyarea里統(tǒng)一后的邊界做裁剪或交集不然緊湊度會被區(qū)域面積變化干擾。4. 跨年份空間分析如何提取“開發(fā)密度”的變化信號三個年份的數(shù)據(jù)放一起最直接的分析是看同一行政單元在不同指標上的走勢。但前提是邊界一致。1_studyarea提供了研究區(qū)邊界而各年份shp的幾何可能來自不同來源必須先處理掉縫隙、重疊與坐標系差異。4.1 統(tǒng)一坐標系與屬性對齊先檢查三個文件的坐標系是否一致import geopandas as gpd study gpd.read_file(./1_studyarea/StudyArea.shp) eco_1998 gpd.read_file(./2c_EconomicDensity/EconomicDensity.shp, encodingutf-8) print(study.crs) print(eco_1998.crs) if eco_1998.crs ! study.crs: eco_1998 eco_1998.to_crs(study.crs)如果.prj文件缺失或定義不完整to_crs會拋錯。應(yīng)急方案是從一個已知正確的GeoDataFrame復制crs對象例如eco_1998.crs study.crs但這樣做的前提是你確認兩套數(shù)據(jù)的橢球和基準面一致否則會出現(xiàn)幾米的偏移。屬性對齊更麻煩。如果2006年的字段結(jié)構(gòu)變了比如新增了“空間GDP”字段而1998年沒有直接用pd.concat會生成大量NaN。我常用的做法是先只保留幾個核心字段統(tǒng)一命名為region、year、value再按區(qū)域名進行merge或groupby。# 讀取緊湊度數(shù)據(jù)并做空間連接 compact gpd.read_file(./2b_DevelopmentCompactness/DevelopmentCompactness.shp) joined gpd.sjoin(study, compact, howleft, predicateintersects) # 按城市名聚合 grouped joined.groupby(CITY_NAME)[COMPACT].mean().reset_index()代碼中sjoin的predicateintersects表示兩個要素只要任意部分相交就算匹配這能規(guī)避邊界微小縫隙。groupby取均值是因為一個城市可能由多個面要素組成。如果希望面積加權(quán)則需要先計算相交面積joined[inter_area] joined.geometry.intersection(study.unary_union).area weighted joined.groupby(CITY_NAME).apply( lambda df: (df[COMPACT] * df[inter_area]).sum() / df[inter_area].sum() )注意apply里的加權(quán)公式當區(qū)域面積差異大時簡單均值會偏向面積小但數(shù)值高的單元加權(quán)均值能反映整體狀態(tài)。4.2 生成三年對比表把三個年份的經(jīng)濟密度整合到一張寬表可以直接用pivot_tableimport pandas as pd # 假設(shè)三個gdf都有YEAR和VALUE字段 eco_1998[year] 1998 eco_2006[year] 2006 eco_2012[year] 2012 frames [eco_1998[[CITY_NAME, year, GDP_KM2]], eco_2006[[CITY_NAME, year, GDP_KM2]], eco_2012[[CITY_NAME, year, GDP_KM2]]] all_df pd.concat(frames) pivot all_df.pivot_table(indexCITY_NAME, columnsyear, valuesGDP_KM2, aggfuncmean) pivot[growth_06_98] pivot[2006] / pivot[1998] - 1 pivot[growth_12_06] pivot[2012] / pivot[2006] - 1pivot_table的aggfunc默認mean它會處理同一城市下多個面要素的重復記錄。用增長率列能直接識別哪幾年擴張最快。4.3 變化可視化可視化可以先用matplotlib畫折線趨勢再用geopandas給每個城市填充增長率顏色。需要注意用顏色分級圖時最好將增長率分位數(shù)分成5類避免極端值把配色拉平。5. 處理1km GDP柵格用rasterio讀取、重投影與統(tǒng)計柵格部分是這份數(shù)據(jù)里最硬的骨頭。NAGDP_1km_1998.tif這類文件是GeoTIFF除了.tif本體還有.tfw世界文件、.aux.xml輔助信息、.ovr金字塔。這套文件對新手不太友好但處理思路很固定先看元數(shù)據(jù)再處理空值最后做統(tǒng)計或重投影。5.1 柵格文件的組成.tif、.tfw、.ovr各干什么.tfw是文本世界文件記錄像元尺寸和左上角坐標.ovr是金字塔文件如果缺失大范圍縮放變慢但不影響數(shù)值精度.aux.xml包含統(tǒng)計信息和nodata標記。讀取時我習慣用rasterio而不是gdal命令行因為能直接返回numpy數(shù)組import rasterio with rasterio.open(./3_NAGDP_1km/NAGDP_1km_1998.tif) as src: data src.read(1) nodata src.nodata profile src.profile print(nodata , nodata) print(數(shù)據(jù)類型 , data.dtype) print(shape , data.shape)src.read(1)讀出第一波段nodata是無值標記。GDP數(shù)據(jù)通常用-9999或0表示無數(shù)據(jù)拿到數(shù)組后第一步要把這些值屏蔽掉否則計入均值會讓結(jié)果嚴重偏低import numpy as np masked np.ma.masked_where(data 0, data) print(masked.mean())5.2 做年度差值并重投影三個年份的tif如果像元尺寸和原點不完全一致不要直接做數(shù)組減法。正確做法是先統(tǒng)一到同一個格網(wǎng)。可以用reproject完成from rasterio.warp import reproject, Resampling, calculate_default_transform with rasterio.open(./3_NAGDP_1km/NAGDP_1km_1998.tif) as src: data_1998 src.read(1) src_crs src.crs src_transform src.transform nodata src.nodata # 目標坐標系如果源是WGS84地理坐標這里用EPSG:4326 transform, width, height calculate_default_transform( src_crs, EPSG:32650, src.width, src.height, *src.bounds ) dst np.zeros((height, width), dtypenp.float32) reproject( sourcedata_1998, destinationdst, src_transformsrc_transform, src_crssrc_crs, dst_transformtransform, dst_crsEPSG:32650, resamplingResampling.bilinear, src_nodatanodata, dst_nodata-9999 )上述代碼中calculate_default_transform根據(jù)源范圍和目標投影計算新尺寸bilinear重采樣適用于連續(xù)型GDP值如果處理分類柵格則改用nearest。src_nodata與dst_nodata保持一致保證重投影后空洞仍是同一標記。如果不想寫這么長的代碼優(yōu)先用rioxarray它將坐標和投影封裝得更好一行就能做差值import rioxarray ds_1998 rioxarray.open_rasterio(./3_NAGDP_1km/NAGDP_1km_1998.tif) ds_2006 rioxarray.open_rasterio(./3_NAGDP_1km/NAGDP_1km_2006.tif) diff ds_2006 - ds_1998但要先確認兩個對象的空間參考一致若不一致則調(diào)用ds_1998.rio.reproject_match(ds_2006)對齊。5.3 高分辨率GDP數(shù)據(jù)的局限與坑1km格網(wǎng)GDP并不是真實觀測值它一般基于夜間燈光、土地利用和統(tǒng)計年鑒做空間化插值因此“某個像元上的GDP”并不代表那個位置真有這么多產(chǎn)出。解讀時更適合用在區(qū)域總量或相對比上。我在做行政單元統(tǒng)計時會先用行政區(qū)邊界裁剪然后逐像元累加。這里最典型的坑有兩個一是邊緣像元被邊界切成窄條如果連通面積小于一個像元直接累加會高估邊界地帶數(shù)值二是.ovr金字塔文件版本不一致可能讓某些軟件讀到錯誤nodata。解決辦法是用rasterstats的all_touchedFalse只保留完整落入邊界的像元并在讀取時顯式指定maskedTrue。提示如果某一年份的GDP均值明顯異常優(yōu)先檢查nodata值是否被默認為0再檢查重投影時是否把邊緣的nodata像元插值成小數(shù)值。6. 收尾技巧用rasterstats一鍵把1km GDP匯總到行政區(qū)前面的操作分別處理了矢量和柵格實際交付成果時常需要“一張Excel表”每一行是一個區(qū)縣列是三個年份的GDP密度與總量。這時用rasterstats包可以省掉大量手寫循環(huán)。它專門用來把柵格數(shù)據(jù)按矢量多邊形進行分區(qū)統(tǒng)計pip install rasterstats然后寫一個極簡函數(shù)import pandas as pd from rasterstats import zonal_stats stats zonal_stats( ./1_studyarea/StudyArea.shp, # 矢量邊界 ./3_NAGDP_1km/NAGDP_1km_1998.tif, # 柵格文件 stats[sum, mean, min, max], nodata-9999, # 與tif元數(shù)據(jù)一致 all_touchedFalse # 只統(tǒng)計完整覆蓋的像元 ) df pd.DataFrame(stats)zonal_stats第一個參數(shù)接收矢量邊界第二個參數(shù)接收柵格路徑stats支持sum/mean/max/min/median等聚合統(tǒng)計。nodata參數(shù)能過濾柵格中的無值像元避免把-9999當作真實數(shù)值。all_touchedFalse表示只統(tǒng)計多邊形中心點位于內(nèi)部的像元這比True更嚴格也更能避免邊緣噪聲。三個年份只需包裝成一個循環(huán)def extract_year(year): raster f./3_NAGDP_1km/NAGDP_1km_{year}.tif st zonal_stats( ./1_studyarea/StudyArea.shp, raster, stats[sum, mean], nodata-9999, all_touchedFalse ) return pd.DataFrame(st) result pd.concat( [extract_year(1998), extract_year(2006), extract_year(2012)], axis1 ) result.columns [GDP_sum_1998, GDP_mean_1998, GDP_sum_2006, GDP_mean_2006, GDP_sum_2012, GDP_mean_2012]這個封裝直接生成一個DataFrame再配合study里的行政區(qū)名稱列即可保存為CSVstudy gpd.read_file(./1_studyarea/StudyArea.shp) result.insert(0, city_name, study[CITY_NAME]) result.to_csv(./prd_gdp_summary.csv, indexFalse)一個小技巧如果發(fā)現(xiàn)某年的sum異常低先回到第5章檢查nodata和重投影步驟或者換個思路用面積加權(quán)平均——因為柵格像元面積可能隨緯度變化而不同尤其在投影坐標系里不同緯度像元實際面積并不完全一致。遇到這種情況可以先用area屬性求出單位面積值再乘實際面積。這套流程跑通后后面再加2018年數(shù)據(jù)只需要復制這個函數(shù)模板替換年份即可。本文還有配套的精品資源點擊獲取