澇建模：下墊面與人口暴露量化)
簡(jiǎn)介北京市建筑物面數(shù)據(jù)提供城市與農(nóng)村建筑輪廓的矢量要素并附帶面積、人口等屬性信息可支撐內(nèi)澇治理、下墊面建模、建筑能耗分析及城鄉(xiāng)規(guī)劃等工作。整套資源打包為rar格式共含六個(gè)文件具備shp主文件、dbf屬性表、prj坐標(biāo)系文件、shx空間索引等完整結(jié)構(gòu)屬于可直接使用的地理信息系統(tǒng)數(shù)據(jù)要素集。壓縮包大小約136.15MB從體量可初步判斷建筑輪廓細(xì)節(jié)較為完整。目前已有七百六十四人學(xué)習(xí)下載常見應(yīng)用包括城市雨洪模擬、建筑密度統(tǒng)計(jì)和可持續(xù)發(fā)展評(píng)估。獲取后可快速提取建筑面要素計(jì)算面積并與人口字段關(guān)聯(lián)開展空間統(tǒng)計(jì)也可作為城市內(nèi)澇模型的基礎(chǔ)底圖節(jié)省資料收集與格式轉(zhuǎn)換時(shí)間適合地理信息分析、環(huán)境工程和城市規(guī)劃從業(yè)者及研究人員使用。1. 內(nèi)澇治理里的下墊面缺口北京建筑面shp數(shù)據(jù)能補(bǔ)上什么做內(nèi)澇治理和SWMM建模的人最頭疼的不是模型參數(shù)率定而是下墊面數(shù)據(jù)。要算清屋頂不透水面積、單元人口暴露量先得有可靠的建筑物輪廓。但遙感影像提取輪廓要折騰一晚公開路網(wǎng)推算不透水面又屬于玄學(xué)。這份北京市建筑物面數(shù)據(jù)是標(biāo)準(zhǔn)的shp格式矢量面圖層建筑面要素、面積、人口信息三樣都齊了。拿到手按研究范圍裁剪再跟SWMM子匯水區(qū)做空間連接屋頂不透水比例和暴露人口就都有了。正在做城市排水模型、內(nèi)澇風(fēng)險(xiǎn)圖或洪澇評(píng)估又不想在矢量化上重復(fù)造輪子的從業(yè)者可以直接往下看怎么落地。2. 打開這份shp之前把坐標(biāo)系、字段口徑、面積基準(zhǔn)三件事釘死shp格式看起來只是拖進(jìn)GIS就能用但真正決定這份數(shù)據(jù)能不能對(duì)接SWMM項(xiàng)目的是加載之后前三分鐘你核查了什么。很多人翻車在坐標(biāo)系對(duì)不上、面積字段單位不符合預(yù)期、人口信息口徑搞混最后結(jié)果直接失真。前面花十分鐘確認(rèn)基準(zhǔn)比后面返工一天強(qiáng)得多。2.1 用GIS加載shp后先看“三樣基礎(chǔ)信息”先把這份北京建筑面數(shù)據(jù)拖進(jìn)QGIS或者ArcGIS。以QGIS為例我從Python控制臺(tái)加載順便把圖層基本狀態(tài)打印出來這樣后面寫腳本處理時(shí)能復(fù)用同一個(gè)圖層對(duì)象。from qgis.core import QgsVectorLayer, QgsProject shp_path /data/beijing_buildings.shp # 按你本機(jī)實(shí)際路徑改 layer QgsVectorLayer(shp_path, beijing_buildings, ogr) if not layer.isValid(): print(加載失敗檢查shp路徑和文件權(quán)限) else: QgsProject.instance().addMapLayer(layer) print(要素?cái)?shù)量:, layer.featureCount()) print(字段列表:, [f.name() for f in layer.fields()]) print(坐標(biāo)系:, layer.crs().authid())這段代碼用QgsVectorLayer把shp包進(jìn)來參數(shù)里第三個(gè)ogr表示讓GDAL/OGR驅(qū)動(dòng)去讀文件。加載成功后重點(diǎn)看三樣?xùn)|西要素?cái)?shù)量決定統(tǒng)計(jì)樣本規(guī)模字段列表決定你之后能用什么屬性坐標(biāo)系決定能不能直接和其他圖層疊加。常見做法是先打開圖層屬性面板再核對(duì)坐標(biāo)系。北京地區(qū)的數(shù)據(jù)經(jīng)常落在CGCS2000 3度帶高斯投影或者WGS 84 UTM 50N這類投影坐標(biāo)系上也有少數(shù)情況給了WGS 84地理坐標(biāo)。我一般會(huì)記住一個(gè)原則洗數(shù)據(jù)時(shí)先把所有圖層統(tǒng)一到同一個(gè)投影坐標(biāo)系。如果你的管網(wǎng)、DEM是UTM建筑面卻停在WGS84經(jīng)緯度后面做緩沖區(qū)、面積計(jì)算都會(huì)出偏差。2.2 屬性表里三類關(guān)鍵字段幾何、面積、人口這份shp的屬性表里至少會(huì)包含三類信息建筑面要素本身、面積字段、人口信息字段。用QGIS打開屬性表你會(huì)看到一行代表一棟建筑或一個(gè)建筑面的要素幾何信息隱含在圖形里屬性字段則把面積和人口掛在邊上。字段類型常見作用容易踩的坑建筑面要素FID/id/OBJECTID關(guān)聯(lián)幾何和屬性作為主鍵數(shù)據(jù)經(jīng)過多次裁剪后ID不連續(xù)外連接時(shí)容易丟行面積字段area/AREA/shape_area描述建筑面面積單位可能是平方米也可能是公頃或平方千米人口字段pop/POP/人口描述建筑對(duì)應(yīng)人口有的按樓棟級(jí)給有的按小區(qū)級(jí)給口徑差異大拿到手如果你發(fā)現(xiàn)字段名是拼音縮寫比如JZMJ、RK也不要慌先用屬性表排序看一眼數(shù)值。面積字段如果最大值上億那單位大概率不是平方米人口字段如果出現(xiàn)很多0和空值說明數(shù)據(jù)生產(chǎn)者只對(duì)居住類建筑做了人口估算。我在實(shí)際項(xiàng)目里會(huì)把字段名先改成自己熟悉的名字避免后面寫SQL或者Python時(shí)搞混。用QGIS圖層屬性里的字段計(jì)算器新建字段即可-- 字段計(jì)算器表達(dá)式新建 area_m2 字段單位轉(zhuǎn)換為平方米 -- 如果原面積字段 JZMJ 單位是平方米直接用 $area $area -- 如果原字段單位是公頃乘 10000 JZMJ * 10000QGIS表達(dá)式里$area返回的是投影坐標(biāo)系下的橢球面積單位取決于圖層當(dāng)前坐標(biāo)系。投影坐標(biāo)系選得好$area就是平方米。如果圖層是地理坐標(biāo)$area的值會(huì)是經(jīng)緯度算出來的度平方那就要先重投影再計(jì)算。2.3 面積口徑建筑基底面積、建筑面積、投影面積不是一個(gè)東西建筑面shp里的面積一般是建筑基底面積也就是建筑落地的那一圈多邊形面積。它不等于容積率口徑下的建筑面積也不同于遙感影像上帶屋頂凸出部分的投影面積。對(duì)SWMM下墊面而言我們要的是“屋頂攔住降雨的水平投影面積”建筑基底面積是最接近這個(gè)物理含義的字段。所以別把shp里的area字段當(dāng)成能直接算容積率的東西。它只能用于地表匯流計(jì)算。如果屬性表里既有面積字段又有層數(shù)字段用來估算屋頂可收集面積時(shí)不要乘層數(shù)屋頂降雨只有一層。人口字段則可以結(jié)合層數(shù)估算一戶幾人那是做內(nèi)澇暴露評(píng)估時(shí)用的和產(chǎn)流計(jì)算是兩回事。我在拿到這種數(shù)據(jù)時(shí)會(huì)先隨機(jī)抽幾棟肉眼可見的大體量建筑和普通住宅對(duì)比屬性面積和圖上量測(cè)面積。如果差在5%以內(nèi)說明面積字段可信如果差得離譜就檢查是不是投影坐標(biāo)系沒配對(duì)。3. 把建筑物輪廓轉(zhuǎn)成SWMM下墊面裁剪、幾何修復(fù)、不透水比例一條線SWMM里面需要的是子匯水區(qū)、不透水比例、特征寬度這些物理參數(shù)。建筑物屋頂本質(zhì)上是不透水地表所以從建筑面數(shù)據(jù)到SWMM參數(shù)的鏈路是把建筑面落到每個(gè)SWMM子匯水區(qū)上算清楚每塊匯水區(qū)里屋頂面積占比然后折算成不透水率。這條鏈路中間有三個(gè)環(huán)節(jié)最容易出問題逐個(gè)說。3.1 先按研究區(qū)裁剪建筑面避免全市數(shù)據(jù)參與統(tǒng)計(jì)SWMM建模范圍通常是某個(gè)排水分區(qū)或某個(gè)內(nèi)澇易發(fā)片區(qū)。如果不裁剪直接把全市建筑面拿來連接會(huì)導(dǎo)致分區(qū)統(tǒng)計(jì)結(jié)果失真。用QGIS的處理框架做一次裁剪順便把下一步要用到的輸出路徑定好。import processing # 輸入已加載的建筑面圖層 buildings_layer研究區(qū)圖層 study_layer # 輸出裁剪后的建筑面 shp processing.run(native:clip, { INPUT: buildings_layer, OVERLAY: study_layer, OUTPUT: /output/buildings_clip.shp })native:clip是QGIS內(nèi)置的裁剪算法INPUT是要被裁剪的圖層OVERLAY是邊界圖層。裁剪完再打印一下要素?cái)?shù)量看是不是比原來少了一大截順便確認(rèn)邊界沒有把建筑切成離破碎。這里有個(gè)細(xì)節(jié)跨在研究區(qū)邊界上的建筑裁剪后會(huì)變成一塊殘片面積比原來小很多。如果你做的是內(nèi)澇風(fēng)險(xiǎn)圖殘片面積會(huì)直接影響匯水區(qū)不透水率常見處理方案是把中心點(diǎn)落在研究區(qū)內(nèi)的建筑整體保留中心點(diǎn)在外的建筑剔除。我是先用質(zhì)心判斷再裁剪這個(gè)邏輯放到避坑章細(xì)說。3.2 修復(fù)幾何和碎多邊形別讓空間連接翻車建筑面數(shù)據(jù)來源多樣常見問題是多邊形自相交、相鄰建筑之間有微小縫隙、共邊不完全重合。直接拿去做空間連接可能出現(xiàn)“這個(gè)面少一塊”“那個(gè)面面積變負(fù)數(shù)”的情況這類問題特別隱蔽不仔細(xì)挑數(shù)據(jù)根本看不見。QGIS里有現(xiàn)成的修復(fù)工具也可以在Python里直接調(diào)用from qgis import processing processing.run(native:fixgeometries, { INPUT: /output/buildings_clip.shp, OUTPUT: /output/buildings_fixed.shp })native:fixgeometries調(diào)用GEOS庫修復(fù)自相交和無效幾何。跑完之后我還習(xí)慣跑一遍“檢查有效幾何”看修復(fù)后還有沒有報(bào)錯(cuò)要素。對(duì)于碎多邊形我會(huì)用“消除”工具按面積閾值合并小于比如10平方米的小碎面直接歸并到相鄰最大面。我一般會(huì)把閾值設(shè)為5到10平方米具體看研究區(qū)建筑密度。老城區(qū)建筑密碎面多閾值可以掛高一點(diǎn)新區(qū)建筑規(guī)整閾值低一些不然會(huì)誤吞合法的小建筑。3.3 用空間連接把屋頂面積攤到SWMM子匯水區(qū)建筑面修復(fù)完成后下一步就是把建筑面和SWMM子匯水區(qū)疊加。這里的關(guān)鍵不是畫一張好看的疊加圖而是按子匯水區(qū)ID做空間連接匯總出每塊匯水區(qū)里的屋頂總面積。import geopandas as gpd build gpd.read_file(/output/buildings_fixed.shp) zones gpd.read_file(/data/swmm_zones.shp) # SWMM劃分好的子匯水區(qū) # 空間連接把子匯水區(qū)信息掛到每個(gè)建筑面上 joined gpd.sjoin(build, zones, howleft, predicateintersects) print(未匹配到子匯水區(qū)的建筑數(shù):, joined[zone_id].isna().sum()) # 按子匯水區(qū)匯總屋頂面積 roof_by_zone joined.groupby(zone_id).agg( roof_area_m2(area_m2, sum) ).reset_index() # 合并回子匯水區(qū)算不透水比例 zones zones.merge(roof_by_zone, onzone_id, howleft) zones[roof_area_m2] zones[roof_area_m2].fillna(0) zones[roof_ratio] zones[roof_area_m2] / zones[zone_area_m2] zones[roof_ratio] zones[roof_ratio].clip(upper0.9) # 保護(hù)性截?cái)噙@段代碼用GeoPandas做空間連接intersects表示只要建筑面和子匯水區(qū)邊界相交就算命中。howleft保證子匯水區(qū)一個(gè)不落。groupby按zone_id匯總屋頂面積最后除以子匯水區(qū)自己的面積得到roof_ratio。參數(shù)說明clip(upper0.9)是防止某些子匯水區(qū)建筑冗余統(tǒng)計(jì)導(dǎo)致比例超過100%。理論上建筑占比不可能超過90%超出多半是邊界裁剪或重疊面沒清理干凈。SWMM里的不透水比例最終要按屋頂、道路、廣場(chǎng)分別累加屋頂只是其中一項(xiàng)所以不要直接把roof_ratio當(dāng)成總不透水比例喂給模型。如果屬性表里人口字段也要參與后續(xù)暴露分析可以在這個(gè)groupby里加上人口匯總pop_by_zone joined.groupby(zone_id).agg( pop_sum(pop, sum) )注意多余的人口字段需要先確認(rèn)不是樓棟級(jí)重復(fù)統(tǒng)計(jì)否則一個(gè)小區(qū)十幾棟樓的人口加起來會(huì)虛高好幾倍。4. 避坑記錄北京建筑面轉(zhuǎn)SWMM下墊面的五個(gè)典型翻車這部分?jǐn)?shù)據(jù)我在不同城市換著花樣踩過坑北京建筑面數(shù)據(jù)的特點(diǎn)是內(nèi)容全、屬性雜、邊界情況多。下面五條只要碰上一次就夠你折騰半天直接抄驗(yàn)證方法。4.1 坐標(biāo)系不一致導(dǎo)致建筑面和管網(wǎng)相差幾公里跑通流程后第一次疊加發(fā)現(xiàn)建筑面和排水管網(wǎng)完全錯(cuò)位距離差到肉眼可見。原因是建筑面shp停在CGCS2000高斯投影而管網(wǎng)數(shù)據(jù)是WGS84經(jīng)緯度兩個(gè)圖層坐標(biāo)系不一致。解決是把所有圖層統(tǒng)一到同一個(gè)投影坐標(biāo)系再疊加用QGIS菜單“重投影圖層”指定 EPSG:4528CGCS2000 / 3-degree Gauss-Kruger zone 29這類投影適合北京經(jīng)度范圍。從那以后我每次建項(xiàng)目第一件事就是建一個(gè)統(tǒng)一的CRS變量所有讀取的圖層強(qiáng)制重投影。4.2 面積字段單位沒核對(duì)算出個(gè)天文數(shù)字在子匯水區(qū)匯總時(shí)不透水比例算出來200%多檢查發(fā)現(xiàn)面積字段原單位是公頃我卻按平方米算。$area返回平方米但你手上字段可能是從原始ArcGIS文件遷移來的AREA單位依然是公頃。解決是看字段值量級(jí)如果一層樓占地面積動(dòng)輒幾萬就要懷疑是平方米寫成千瓦時(shí)的那種低級(jí)錯(cuò)誤。正確做法是先統(tǒng)計(jì)面積字段的min/max再和實(shí)際建筑體量比對(duì)確認(rèn)無誤后再進(jìn)入計(jì)算。4.3 人口字段是樓棟級(jí)不是每棟建筑都有值屬性表里pop字段有大量空值有的樓是0有的樓是幾百。因?yàn)槿丝诮y(tǒng)計(jì)到樓棟級(jí)只有住宅樓掛數(shù)配套公建、車棚沒有。直接用groupby().sum()匯總會(huì)嚴(yán)重低估實(shí)際人口因?yàn)榭罩当粊G棄。解決是先用fillna(0)看總量再疊加社區(qū)戶數(shù)校正系數(shù)。碰到明顯低于街道統(tǒng)計(jì)的時(shí)候就把研究區(qū)內(nèi)的人口字段按建筑面積加權(quán)攤回去而不是用shp里的原始值直接進(jìn)結(jié)論。4.4 跨邊界建筑被裁剪后形成碎面污染匯水區(qū)統(tǒng)計(jì)前面講跨邊界建筑要用中心點(diǎn)判斷實(shí)際操作中直接裁剪會(huì)留下一條條極薄的殘片這些殘片面積可能只有2、3平方米。匯總時(shí)殘片本身影響不大但它們會(huì)把相鄰子匯水區(qū)的邊界“縫”起來空間連接的時(shí)候把這些殘片錯(cuò)分到隔壁匯水區(qū)。解決方式是在裁剪之前先給建筑面計(jì)算質(zhì)心坐標(biāo)按質(zhì)心落入哪個(gè)研究區(qū)來分配整棟建筑不做什么裁剪減少一個(gè)大坑。4.5 建筑面里混著非居住類建筑人口字段全空學(xué)校操場(chǎng)輪廓、廠房、車棚、配電房都在建筑面圖層里它們面積很大但人口字段是空的。SWMM計(jì)算不透水面積時(shí)操場(chǎng)和廠房其實(shí)也算不透水或半透水保留沒問題。但做人口暴露評(píng)估時(shí)這些“空人口”建筑會(huì)導(dǎo)致平均人口密度被拉低。解決是加一個(gè)building_type判斷字段把居住類建筑單獨(dú)篩出來做人口匯總其他建筑只參與面積統(tǒng)計(jì)。5. 把人口字段接進(jìn)內(nèi)澇風(fēng)險(xiǎn)按建筑面估算單元人口暴露量建筑面數(shù)據(jù)最后能產(chǎn)出一個(gè)很有價(jià)值的東西不是屋頂面積而是“最容易淹到人的地方到底有多少人”。SWMM只算水算出積水深度后還缺一張“受淹人口”圖這時(shí)人口字段就派上用場(chǎng)了。做法是先用建筑面的pop字段按子匯水區(qū)匯總暴露人口再把SWMM模擬出的淹沒范圍疊加到人口分布上就能得到不同重現(xiàn)期下受影響人口量級(jí)。import geopandas as gpd # SWMM模擬得到的淹沒區(qū)保存為 inundation.shp inundation gpd.read_file(/output/inundation.shp) pop_by_zone gpd.read_file(/output/pop_by_zone.shp) # 上一章算出的匯水區(qū)人口 # 淹沒區(qū)與人口統(tǒng)計(jì)區(qū)做空間連接按淹沒面積比例折算人口 hit gpd.overlay(inundation, pop_by_zone, howintersection) hit[exposed_pop] hit[pop_sum] * (hit.geometry.area / hit[zone_area_m2]) print(受影響人口:, hit[exposed_pop].sum()) print(受影響匯水區(qū)數(shù)量:, hit[zone_id].nunique())這段代碼核心就一個(gè)分配邏輯淹沒區(qū)面積占這個(gè)子匯水區(qū)的比例乘以該區(qū)總?cè)丝诘玫奖┞度丝?。如果淹沒區(qū)只淹了一小片人口按面積均勻攤?cè)绻麄€(gè)匯水區(qū)都在水里暴露人口就等于總?cè)丝?。精度夠用不需要做到單棟樓逐戶?jí)別。最后說一個(gè)習(xí)慣腳本跑完我不會(huì)直接拿去匯報(bào)而是把結(jié)果和街道人口普查數(shù)量級(jí)做一次對(duì)比。建筑面數(shù)據(jù)加工出來的人口只是相對(duì)分布絕對(duì)值和普查數(shù)差30%以上就說明字段口徑有問題需要回去看是按戶籍算還是常住算。從那以后我每次交付前都強(qiáng)制走一遍“面積字段單位核查——人口口徑核查——淹沒比例驗(yàn)證”這條天路少一步都不敢出圖。這方法救了我好幾次數(shù)據(jù)事故希望幫到你。本文還有配套的精品資源點(diǎn)擊獲取 從北京建筑面shp到SWMM內(nèi)澇建模：下墊面與人口暴露量化
簡(jiǎn)介北京市建筑物面數(shù)據(jù)提供城市與農(nóng)村建筑輪廓的矢量要素并附帶面積、人口等屬性信息可支撐內(nèi)澇治理、下墊面建模、建筑能耗分析及城鄉(xiāng)規(guī)劃等工作。整套資源打包為rar格式共含六個(gè)文件具備shp主文件、dbf屬性表、prj坐標(biāo)系文件、shx空間索引等完整結(jié)構(gòu)屬于可直接使用的地理信息系統(tǒng)數(shù)據(jù)要素集。壓縮包大小約136.15MB從體量可初步判斷建筑輪廓細(xì)節(jié)較為完整。目前已有七百六十四人學(xué)習(xí)下載常見應(yīng)用包括城市雨洪模擬、建筑密度統(tǒng)計(jì)和可持續(xù)發(fā)展評(píng)估。獲取后可快速提取建筑面要素計(jì)算面積并與人口字段關(guān)聯(lián)開展空間統(tǒng)計(jì)也可作為城市內(nèi)澇模型的基礎(chǔ)底圖節(jié)省資料收集與格式轉(zhuǎn)換時(shí)間適合地理信息分析、環(huán)境工程和城市規(guī)劃從業(yè)者及研究人員使用。1. 內(nèi)澇治理里的下墊面缺口北京建筑面shp數(shù)據(jù)能補(bǔ)上什么做內(nèi)澇治理和SWMM建模的人最頭疼的不是模型參數(shù)率定而是下墊面數(shù)據(jù)。要算清屋頂不透水面積、單元人口暴露量先得有可靠的建筑物輪廓。但遙感影像提取輪廓要折騰一晚公開路網(wǎng)推算不透水面又屬于玄學(xué)。這份北京市建筑物面數(shù)據(jù)是標(biāo)準(zhǔn)的shp格式矢量面圖層建筑面要素、面積、人口信息三樣都齊了。拿到手按研究范圍裁剪再跟SWMM子匯水區(qū)做空間連接屋頂不透水比例和暴露人口就都有了。正在做城市排水模型、內(nèi)澇風(fēng)險(xiǎn)圖或洪澇評(píng)估又不想在矢量化上重復(fù)造輪子的從業(yè)者可以直接往下看怎么落地。2. 打開這份shp之前把坐標(biāo)系、字段口徑、面積基準(zhǔn)三件事釘死shp格式看起來只是拖進(jìn)GIS就能用但真正決定這份數(shù)據(jù)能不能對(duì)接SWMM項(xiàng)目的是加載之后前三分鐘你核查了什么。很多人翻車在坐標(biāo)系對(duì)不上、面積字段單位不符合預(yù)期、人口信息口徑搞混最后結(jié)果直接失真。前面花十分鐘確認(rèn)基準(zhǔn)比后面返工一天強(qiáng)得多。2.1 用GIS加載shp后先看“三樣基礎(chǔ)信息”先把這份北京建筑面數(shù)據(jù)拖進(jìn)QGIS或者ArcGIS。以QGIS為例我從Python控制臺(tái)加載順便把圖層基本狀態(tài)打印出來這樣后面寫腳本處理時(shí)能復(fù)用同一個(gè)圖層對(duì)象。from qgis.core import QgsVectorLayer, QgsProject shp_path /data/beijing_buildings.shp # 按你本機(jī)實(shí)際路徑改 layer QgsVectorLayer(shp_path, beijing_buildings, ogr) if not layer.isValid(): print(加載失敗檢查shp路徑和文件權(quán)限) else: QgsProject.instance().addMapLayer(layer) print(要素?cái)?shù)量:, layer.featureCount()) print(字段列表:, [f.name() for f in layer.fields()]) print(坐標(biāo)系:, layer.crs().authid())這段代碼用QgsVectorLayer把shp包進(jìn)來參數(shù)里第三個(gè)ogr表示讓GDAL/OGR驅(qū)動(dòng)去讀文件。加載成功后重點(diǎn)看三樣?xùn)|西要素?cái)?shù)量決定統(tǒng)計(jì)樣本規(guī)模字段列表決定你之后能用什么屬性坐標(biāo)系決定能不能直接和其他圖層疊加。常見做法是先打開圖層屬性面板再核對(duì)坐標(biāo)系。北京地區(qū)的數(shù)據(jù)經(jīng)常落在CGCS2000 3度帶高斯投影或者WGS 84 UTM 50N這類投影坐標(biāo)系上也有少數(shù)情況給了WGS 84地理坐標(biāo)。我一般會(huì)記住一個(gè)原則洗數(shù)據(jù)時(shí)先把所有圖層統(tǒng)一到同一個(gè)投影坐標(biāo)系。如果你的管網(wǎng)、DEM是UTM建筑面卻停在WGS84經(jīng)緯度后面做緩沖區(qū)、面積計(jì)算都會(huì)出偏差。2.2 屬性表里三類關(guān)鍵字段幾何、面積、人口這份shp的屬性表里至少會(huì)包含三類信息建筑面要素本身、面積字段、人口信息字段。用QGIS打開屬性表你會(huì)看到一行代表一棟建筑或一個(gè)建筑面的要素幾何信息隱含在圖形里屬性字段則把面積和人口掛在邊上。字段類型常見作用容易踩的坑建筑面要素FID/id/OBJECTID關(guān)聯(lián)幾何和屬性作為主鍵數(shù)據(jù)經(jīng)過多次裁剪后ID不連續(xù)外連接時(shí)容易丟行面積字段area/AREA/shape_area描述建筑面面積單位可能是平方米也可能是公頃或平方千米人口字段pop/POP/人口描述建筑對(duì)應(yīng)人口有的按樓棟級(jí)給有的按小區(qū)級(jí)給口徑差異大拿到手如果你發(fā)現(xiàn)字段名是拼音縮寫比如JZMJ、RK也不要慌先用屬性表排序看一眼數(shù)值。面積字段如果最大值上億那單位大概率不是平方米人口字段如果出現(xiàn)很多0和空值說明數(shù)據(jù)生產(chǎn)者只對(duì)居住類建筑做了人口估算。我在實(shí)際項(xiàng)目里會(huì)把字段名先改成自己熟悉的名字避免后面寫SQL或者Python時(shí)搞混。用QGIS圖層屬性里的字段計(jì)算器新建字段即可-- 字段計(jì)算器表達(dá)式新建 area_m2 字段單位轉(zhuǎn)換為平方米 -- 如果原面積字段 JZMJ 單位是平方米直接用 $area $area -- 如果原字段單位是公頃乘 10000 JZMJ * 10000QGIS表達(dá)式里$area返回的是投影坐標(biāo)系下的橢球面積單位取決于圖層當(dāng)前坐標(biāo)系。投影坐標(biāo)系選得好$area就是平方米。如果圖層是地理坐標(biāo)$area的值會(huì)是經(jīng)緯度算出來的度平方那就要先重投影再計(jì)算。2.3 面積口徑建筑基底面積、建筑面積、投影面積不是一個(gè)東西建筑面shp里的面積一般是建筑基底面積也就是建筑落地的那一圈多邊形面積。它不等于容積率口徑下的建筑面積也不同于遙感影像上帶屋頂凸出部分的投影面積。對(duì)SWMM下墊面而言我們要的是“屋頂攔住降雨的水平投影面積”建筑基底面積是最接近這個(gè)物理含義的字段。所以別把shp里的area字段當(dāng)成能直接算容積率的東西。它只能用于地表匯流計(jì)算。如果屬性表里既有面積字段又有層數(shù)字段用來估算屋頂可收集面積時(shí)不要乘層數(shù)屋頂降雨只有一層。人口字段則可以結(jié)合層數(shù)估算一戶幾人那是做內(nèi)澇暴露評(píng)估時(shí)用的和產(chǎn)流計(jì)算是兩回事。我在拿到這種數(shù)據(jù)時(shí)會(huì)先隨機(jī)抽幾棟肉眼可見的大體量建筑和普通住宅對(duì)比屬性面積和圖上量測(cè)面積。如果差在5%以內(nèi)說明面積字段可信如果差得離譜就檢查是不是投影坐標(biāo)系沒配對(duì)。3. 把建筑物輪廓轉(zhuǎn)成SWMM下墊面裁剪、幾何修復(fù)、不透水比例一條線SWMM里面需要的是子匯水區(qū)、不透水比例、特征寬度這些物理參數(shù)。建筑物屋頂本質(zhì)上是不透水地表所以從建筑面數(shù)據(jù)到SWMM參數(shù)的鏈路是把建筑面落到每個(gè)SWMM子匯水區(qū)上算清楚每塊匯水區(qū)里屋頂面積占比然后折算成不透水率。這條鏈路中間有三個(gè)環(huán)節(jié)最容易出問題逐個(gè)說。3.1 先按研究區(qū)裁剪建筑面避免全市數(shù)據(jù)參與統(tǒng)計(jì)SWMM建模范圍通常是某個(gè)排水分區(qū)或某個(gè)內(nèi)澇易發(fā)片區(qū)。如果不裁剪直接把全市建筑面拿來連接會(huì)導(dǎo)致分區(qū)統(tǒng)計(jì)結(jié)果失真。用QGIS的處理框架做一次裁剪順便把下一步要用到的輸出路徑定好。import processing # 輸入已加載的建筑面圖層 buildings_layer研究區(qū)圖層 study_layer # 輸出裁剪后的建筑面 shp processing.run(native:clip, { INPUT: buildings_layer, OVERLAY: study_layer, OUTPUT: /output/buildings_clip.shp })native:clip是QGIS內(nèi)置的裁剪算法INPUT是要被裁剪的圖層OVERLAY是邊界圖層。裁剪完再打印一下要素?cái)?shù)量看是不是比原來少了一大截順便確認(rèn)邊界沒有把建筑切成離破碎。這里有個(gè)細(xì)節(jié)跨在研究區(qū)邊界上的建筑裁剪后會(huì)變成一塊殘片面積比原來小很多。如果你做的是內(nèi)澇風(fēng)險(xiǎn)圖殘片面積會(huì)直接影響匯水區(qū)不透水率常見處理方案是把中心點(diǎn)落在研究區(qū)內(nèi)的建筑整體保留中心點(diǎn)在外的建筑剔除。我是先用質(zhì)心判斷再裁剪這個(gè)邏輯放到避坑章細(xì)說。3.2 修復(fù)幾何和碎多邊形別讓空間連接翻車建筑面數(shù)據(jù)來源多樣常見問題是多邊形自相交、相鄰建筑之間有微小縫隙、共邊不完全重合。直接拿去做空間連接可能出現(xiàn)“這個(gè)面少一塊”“那個(gè)面面積變負(fù)數(shù)”的情況這類問題特別隱蔽不仔細(xì)挑數(shù)據(jù)根本看不見。QGIS里有現(xiàn)成的修復(fù)工具也可以在Python里直接調(diào)用from qgis import processing processing.run(native:fixgeometries, { INPUT: /output/buildings_clip.shp, OUTPUT: /output/buildings_fixed.shp })native:fixgeometries調(diào)用GEOS庫修復(fù)自相交和無效幾何。跑完之后我還習(xí)慣跑一遍“檢查有效幾何”看修復(fù)后還有沒有報(bào)錯(cuò)要素。對(duì)于碎多邊形我會(huì)用“消除”工具按面積閾值合并小于比如10平方米的小碎面直接歸并到相鄰最大面。我一般會(huì)把閾值設(shè)為5到10平方米具體看研究區(qū)建筑密度。老城區(qū)建筑密碎面多閾值可以掛高一點(diǎn)新區(qū)建筑規(guī)整閾值低一些不然會(huì)誤吞合法的小建筑。3.3 用空間連接把屋頂面積攤到SWMM子匯水區(qū)建筑面修復(fù)完成后下一步就是把建筑面和SWMM子匯水區(qū)疊加。這里的關(guān)鍵不是畫一張好看的疊加圖而是按子匯水區(qū)ID做空間連接匯總出每塊匯水區(qū)里的屋頂總面積。import geopandas as gpd build gpd.read_file(/output/buildings_fixed.shp) zones gpd.read_file(/data/swmm_zones.shp) # SWMM劃分好的子匯水區(qū) # 空間連接把子匯水區(qū)信息掛到每個(gè)建筑面上 joined gpd.sjoin(build, zones, howleft, predicateintersects) print(未匹配到子匯水區(qū)的建筑數(shù):, joined[zone_id].isna().sum()) # 按子匯水區(qū)匯總屋頂面積 roof_by_zone joined.groupby(zone_id).agg( roof_area_m2(area_m2, sum) ).reset_index() # 合并回子匯水區(qū)算不透水比例 zones zones.merge(roof_by_zone, onzone_id, howleft) zones[roof_area_m2] zones[roof_area_m2].fillna(0) zones[roof_ratio] zones[roof_area_m2] / zones[zone_area_m2] zones[roof_ratio] zones[roof_ratio].clip(upper0.9) # 保護(hù)性截?cái)噙@段代碼用GeoPandas做空間連接intersects表示只要建筑面和子匯水區(qū)邊界相交就算命中。howleft保證子匯水區(qū)一個(gè)不落。groupby按zone_id匯總屋頂面積最后除以子匯水區(qū)自己的面積得到roof_ratio。參數(shù)說明clip(upper0.9)是防止某些子匯水區(qū)建筑冗余統(tǒng)計(jì)導(dǎo)致比例超過100%。理論上建筑占比不可能超過90%超出多半是邊界裁剪或重疊面沒清理干凈。SWMM里的不透水比例最終要按屋頂、道路、廣場(chǎng)分別累加屋頂只是其中一項(xiàng)所以不要直接把roof_ratio當(dāng)成總不透水比例喂給模型。如果屬性表里人口字段也要參與后續(xù)暴露分析可以在這個(gè)groupby里加上人口匯總pop_by_zone joined.groupby(zone_id).agg( pop_sum(pop, sum) )注意多余的人口字段需要先確認(rèn)不是樓棟級(jí)重復(fù)統(tǒng)計(jì)否則一個(gè)小區(qū)十幾棟樓的人口加起來會(huì)虛高好幾倍。4. 避坑記錄北京建筑面轉(zhuǎn)SWMM下墊面的五個(gè)典型翻車這部分?jǐn)?shù)據(jù)我在不同城市換著花樣踩過坑北京建筑面數(shù)據(jù)的特點(diǎn)是內(nèi)容全、屬性雜、邊界情況多。下面五條只要碰上一次就夠你折騰半天直接抄驗(yàn)證方法。4.1 坐標(biāo)系不一致導(dǎo)致建筑面和管網(wǎng)相差幾公里跑通流程后第一次疊加發(fā)現(xiàn)建筑面和排水管網(wǎng)完全錯(cuò)位距離差到肉眼可見。原因是建筑面shp停在CGCS2000高斯投影而管網(wǎng)數(shù)據(jù)是WGS84經(jīng)緯度兩個(gè)圖層坐標(biāo)系不一致。解決是把所有圖層統(tǒng)一到同一個(gè)投影坐標(biāo)系再疊加用QGIS菜單“重投影圖層”指定 EPSG:4528CGCS2000 / 3-degree Gauss-Kruger zone 29這類投影適合北京經(jīng)度范圍。從那以后我每次建項(xiàng)目第一件事就是建一個(gè)統(tǒng)一的CRS變量所有讀取的圖層強(qiáng)制重投影。4.2 面積字段單位沒核對(duì)算出個(gè)天文數(shù)字在子匯水區(qū)匯總時(shí)不透水比例算出來200%多檢查發(fā)現(xiàn)面積字段原單位是公頃我卻按平方米算。$area返回平方米但你手上字段可能是從原始ArcGIS文件遷移來的AREA單位依然是公頃。解決是看字段值量級(jí)如果一層樓占地面積動(dòng)輒幾萬就要懷疑是平方米寫成千瓦時(shí)的那種低級(jí)錯(cuò)誤。正確做法是先統(tǒng)計(jì)面積字段的min/max再和實(shí)際建筑體量比對(duì)確認(rèn)無誤后再進(jìn)入計(jì)算。4.3 人口字段是樓棟級(jí)不是每棟建筑都有值屬性表里pop字段有大量空值有的樓是0有的樓是幾百。因?yàn)槿丝诮y(tǒng)計(jì)到樓棟級(jí)只有住宅樓掛數(shù)配套公建、車棚沒有。直接用groupby().sum()匯總會(huì)嚴(yán)重低估實(shí)際人口因?yàn)榭罩当粊G棄。解決是先用fillna(0)看總量再疊加社區(qū)戶數(shù)校正系數(shù)。碰到明顯低于街道統(tǒng)計(jì)的時(shí)候就把研究區(qū)內(nèi)的人口字段按建筑面積加權(quán)攤回去而不是用shp里的原始值直接進(jìn)結(jié)論。4.4 跨邊界建筑被裁剪后形成碎面污染匯水區(qū)統(tǒng)計(jì)前面講跨邊界建筑要用中心點(diǎn)判斷實(shí)際操作中直接裁剪會(huì)留下一條條極薄的殘片這些殘片面積可能只有2、3平方米。匯總時(shí)殘片本身影響不大但它們會(huì)把相鄰子匯水區(qū)的邊界“縫”起來空間連接的時(shí)候把這些殘片錯(cuò)分到隔壁匯水區(qū)。解決方式是在裁剪之前先給建筑面計(jì)算質(zhì)心坐標(biāo)按質(zhì)心落入哪個(gè)研究區(qū)來分配整棟建筑不做什么裁剪減少一個(gè)大坑。4.5 建筑面里混著非居住類建筑人口字段全空學(xué)校操場(chǎng)輪廓、廠房、車棚、配電房都在建筑面圖層里它們面積很大但人口字段是空的。SWMM計(jì)算不透水面積時(shí)操場(chǎng)和廠房其實(shí)也算不透水或半透水保留沒問題。但做人口暴露評(píng)估時(shí)這些“空人口”建筑會(huì)導(dǎo)致平均人口密度被拉低。解決是加一個(gè)building_type判斷字段把居住類建筑單獨(dú)篩出來做人口匯總其他建筑只參與面積統(tǒng)計(jì)。5. 把人口字段接進(jìn)內(nèi)澇風(fēng)險(xiǎn)按建筑面估算單元人口暴露量建筑面數(shù)據(jù)最后能產(chǎn)出一個(gè)很有價(jià)值的東西不是屋頂面積而是“最容易淹到人的地方到底有多少人”。SWMM只算水算出積水深度后還缺一張“受淹人口”圖這時(shí)人口字段就派上用場(chǎng)了。做法是先用建筑面的pop字段按子匯水區(qū)匯總暴露人口再把SWMM模擬出的淹沒范圍疊加到人口分布上就能得到不同重現(xiàn)期下受影響人口量級(jí)。import geopandas as gpd # SWMM模擬得到的淹沒區(qū)保存為 inundation.shp inundation gpd.read_file(/output/inundation.shp) pop_by_zone gpd.read_file(/output/pop_by_zone.shp) # 上一章算出的匯水區(qū)人口 # 淹沒區(qū)與人口統(tǒng)計(jì)區(qū)做空間連接按淹沒面積比例折算人口 hit gpd.overlay(inundation, pop_by_zone, howintersection) hit[exposed_pop] hit[pop_sum] * (hit.geometry.area / hit[zone_area_m2]) print(受影響人口:, hit[exposed_pop].sum()) print(受影響匯水區(qū)數(shù)量:, hit[zone_id].nunique())這段代碼核心就一個(gè)分配邏輯淹沒區(qū)面積占這個(gè)子匯水區(qū)的比例乘以該區(qū)總?cè)丝诘玫奖┞度丝?。如果淹沒區(qū)只淹了一小片人口按面積均勻攤?cè)绻麄€(gè)匯水區(qū)都在水里暴露人口就等于總?cè)丝?。精度夠用不需要做到單棟樓逐戶?jí)別。最后說一個(gè)習(xí)慣腳本跑完我不會(huì)直接拿去匯報(bào)而是把結(jié)果和街道人口普查數(shù)量級(jí)做一次對(duì)比。建筑面數(shù)據(jù)加工出來的人口只是相對(duì)分布絕對(duì)值和普查數(shù)差30%以上就說明字段口徑有問題需要回去看是按戶籍算還是常住算。從那以后我每次交付前都強(qiáng)制走一遍“面積字段單位核查——人口口徑核查——淹沒比例驗(yàn)證”這條天路少一步都不敢出圖。這方法救了我好幾次數(shù)據(jù)事故希望幫到你。本文還有配套的精品資源點(diǎn)擊獲取