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