简介以重庆市城市与农村建筑物面数据为核心的shp格式资源面要素内附建筑面积与人口等属性字段适用于内涝治理、SWMM建模、智慧城市评估与碳中和背景下的下垫面分析。数据整体约309MB共6个文件包含.shp空间几何、.dbf属性表、.prj坐标参考、.shx索引及配套元数据结构标准可直接加载至ArcGIS/QGIS。已有919人浏览学习适合城乡规划、市政水务和建筑能耗研究方向的初学者与工程师使用。借助该数据可批量提取屋顶轮廓与建筑密度为海绵城市模拟、片区人口估算及建筑能源模型提供基础底图省去自行采集与矢量化处理的环节。1. 重庆建筑物面数据城市农村shp不止是轮廓还带着面积和人口属性在GIS项目里拿到“重庆市建筑物面数据城市农村shp格式包含建筑面要素、面积、人口信息”这类数据时最常见的错误是把它当成普通建筑轮廓直接叠加底图画两笔专题图就完事。实际上这份数据里最有价值的是两列属性面积和人口。面积能帮你算容积率、户均面积、建造成本估算人口字段如果可靠可以整合进应急疏散、网格化治理、甚至“15分钟生活圈”补短板分析。它覆盖城市和农村意味着坐标系、精度、字段完整性在两类区域往往不一致不能一套参数走到底。这篇文章就从结构读取、面积修正、人口匹配到导出验证走完整条你应该对这份数据做的事。2. 读懂shp里的建筑面、面积字段和人口字段先解决字段混乱的问题2.1 用geopandas读出重庆建筑数据的基本结构与常见属性拿到shp第一件事不是打开QGIS看图层而是先看字段字典和空间范围。常见做法是用Python的geopandas读取打印字段名、前几行和坐标系用30秒判断这份数据有没有经过拆分或投影转换。import geopandas as gpd # 读取重庆建筑面shp假定路径为cq_buildings.shp gdf gpd.read_file(cq_buildings.shp, encodingutf-8) # 打印全部字段名、行数和坐标系 print(字段列表:, gdf.columns.tolist()) print(要素数量:, len(gdf)) print(坐标系:, gdf.crs) # 查看面积和人口相关字段的前5行判断字段实际内容 cols [c for c in gdf.columns if any(k in c.lower() for k in [area, mj, pop, rk, jzmj, hjs, hhs])] print(疑似面积/人口字段:, cols) print(gdf[cols].head())这段代码逻辑先用encoding参数应对图层属性表的乱码因为shp的dbf编码常见GBK或UTF-8再通过关键词筛出面积和人口字段避免自己一个个猜。注意read_file默认会读取全部几何大型文件可以先只读属性gdf gpd.read_file(path, ignore_geometryTrue)但那样就拿不到几何信息只适合快速查字段。参数说明encoding设错会导致中文乱码建议先试utf-8出现乱码再换gbkignore_geometryTrue时返回的GeoDataFrame没有geometry列后续无法做空间操作。实际项目中我一般先把属性读出来记录字段名后再用完整读取去做空间分析。2.2 面积字段的三种来源字段现成、几何计算、单位换算这份数据里面积可能来自三处一是数据生产方在dbf里直接写好的area/mj/jzmj字段二是AREA字段是ArcGIS自动生成但不带单位三是有字段但单位是平方千米而你要的是平方米。读数据后先做一致性检查用几何算一遍面积跟字段比差值。如果差值超过5%说明字段是原始等级路网时代的估算值不可直接用。常见的面积字段名称对照字段名抄录真实含义单位常见坑AREA图形面积平方米或平方千米ArcGIS导出时继承的旧值常单位错JZMJ建筑面积平方米可能与占地面积混淆CQMJ产权面积平方米多用于房产农村字段多为空POP/RK人口数人可能是户籍人口不是常住人口HJS户籍户数户农村常用城市缺失率高先别急着删除或覆盖字段建议新增一列area_geom用几何计算得出标准面积与原始area字段并存。这样后续分析可以随时回溯。2.3 人口信息字段的常见形态与缺失处理人口信息在这类shp里有两种存在方式一种是每个建筑面要素都挂一个人口数另一种是只有少数街区/村庄面要素带人口建筑面里是空的。拿重庆农村房屋来说很多数据源以“户数”代替“人口数”因为村民宅基地确权时登记的是户主与家庭成员数但统计口径五花八门有的含老人有的不含儿童。城市建筑则常见“居住人数估算值”由容积率乘以户均人数推出。如果人口字段全空你需要放弃“按建筑直接统计人口”的幻想改用一个靠谱的替代方案把重庆乡镇或街道边界的人口普查数据通过空间连接分摊到每个建筑面上。这个方案在第4章展开。目前这个阶段你要做的是把人口字段类型统一为float并记录缺失率# 统一人口字段名避免大小写和别名干扰 pop_cols [c for c in gdf.columns if c.upper() in (POP, POPULATION, RK, ZRK)] if pop_cols: gdf[pop_clean] gdf[pop_cols[0]].astype(float) else: gdf[pop_clean] None missing_rate gdf[pop_clean].isna().mean() print(f人口字段缺失率: {missing_rate:.1%}) # 缺失率超过20%就必须考虑外部人口数据做空间匹配这段代码先把疑似人口字段统一成float算出缺失率。缺失率超过20%就不建议用户直接按字段汇总了。参数说明astype(float)会把1,024这类带千分位字符串报错需要先gdf[pop_cols[0]] gdf[pop_cols[0]].astype(str).str.replace(,, )再转换。人口字段的清洗必须先于任何统计否则后面所有按人口汇总的结果都会带着隐性错误。3. 用正确的投影重算建筑面面积别被shp自带面积骗了3.1 为什么重庆地理坐标系下面积会偏小重庆地处东经105°17110°11、北纬28°1032°13如果shp的坐标系是WGS84或CGCS2000地理坐标系度为单位那么计算面积时任意GIS软件都会采用经纬度球面算法结果会是平方米级别的小数或者软件默认给一个虚假的“度”单位没人能看懂。更麻烦的是很多shp虽然显示为GCS_WGS_1984实际数据却来自高斯-克吕格投影的成果只是导出时丢了投影信息。正确的做法是把它投影到重庆当地常用的高斯-克吕格3度带投影。国家标准中重庆主城多用CGCS2000 / 3-degree Gauss-Kruger zone 35EPSG代码为4547也可用4548针对区县。有些团队用WGS 84 / UTM zone 48NEPSG:32648也是可行的但高斯投影同一条带内面积变形更小。3.2 用Python按高斯-克吕格投影计算椭球面积GeoPandas里最简单的做法是to_crs到EPSG:4547然后用area属性。area的值是该投影平面下的面积单位为平方米但它不是真正的椭球面积而是投影面积。对于50平方公里以下的县级尺度投影面积与椭球面积误差小于千分之一对重庆这种山地地形如果只做相对排名和比例分析完全够用。import geopandas as gpd gdf gpd.read_file(cq_buildings.shp, encodingutf-8) # 先统一几何有效性避免自相交和多边形环错误导致面积异常 gdf[geometry] gdf.geometry.buffer(0) # 投影到重庆高斯3度带CGCS2000, zone 35 gdf_proj gdf.to_crs(epsg4547) # 计算投影平面面积单位平方米 gdf_proj[area_m2_proj] gdf_proj.geometry.area # 再算椭球面积用pyproj的Geod from pyproj import Geod geod Geod(ellpsWGS84) def geodesic_area(geom): from shapely.geometry import Polygon, MultiPolygon if geom.is_empty: return 0.0 if geom.geom_type Polygon: x, y geom.exterior.xy area_m2, _ geod.polygon_area_perimeter(list(y), list(x)) return abs(area_m2) elif geom.geom_type MultiPolygon: total 0.0 for part in geom.geoms: total geodesic_area(part) return total else: return 0.0 gdf_proj[area_m2_geo] gdf_proj.apply(lambda r: geodesic_area(r[geometry]), axis1)这段代码有两组面积area_m2_proj是投影平面面积area_m2_geo是WGS84椭球下的测地线面积。因为坐标系是CGCS2000和WGS84在米级精度下差异极小直接用WGS84椭球参数没问题。注意geod.polygon_area_perimeter的输入是(纬, 经)的顺序也就是第一个参数是ylat第二个是xlon写反会得到荒谬值。buffer(0)这一步不是给几何加缓冲而是利用Shapely的容错机制修复细小的自相交和无效环使面积计算不会返回NaN。提示如果shp里有多部件建筑比如一排联排房屋被合并为一个MultiPolygon上述geodesic_area会把每个子多边形面积相加正合适。不要先dissolve否则总面积相同但会丢掉单栋建筑的独立性。# 如果存在无效几何对每个要素单独检查 invalid ~gdf_proj.geometry.is_valid print(f无效几何数量: {invalid.sum()}) if invalid.sum() 0: gdf_proj.loc[invalid, geometry] gdf_proj.loc[invalid, geometry].buffer(0)这段检查代码放在面积计算的任意阶段都可以主要确保area_m2_geo计算时不会因为单条自相交几何报错。注意buffer(0)会替换原几何如果后续还要做空间连接建议保留一个原始几何副本比如gdf_proj[geom_orig] gdf_proj.geometry因为buffer可能把极窄的建筑缝隙合并掉。3.3 按城市/农村分组汇总面积验证数据是否自洽拿到两种面积后顺手做一次分组汇总既验证数据完整性也给后续人口分摊提供分母。假如shp里带一个“城乡分类”字段如CLSD或XZQMC按它分组如果没有分类字段就按村级行政区划点做空间叠加来给每栋建筑打标。# 假设存在city_rural字段取值urban/rural if city_rural in gdf_proj.columns: df_group gdf_proj.groupby(city_rural).agg( 建筑数(geometry, count), 投影面积_m2(area_m2_proj, sum), 椭球面积_m2(area_m2_geo, sum), 人口_缺失率(pop_clean, lambda s: s.isna().mean()) ).reset_index() print(df_group) else: # 没有城乡字段时用土地利用或村庄点数据做空间连接 gdf_proj[city_rural] gdf_proj[XZQMC].str.contains(村|组, naFalse).map({True: rural, False: urban})这个聚合结果能告诉你三件事农村建筑数量是否远多于城市、人口字段缺失率是否在城乡间差异巨大、投影面积和椭球面积差值是否一致。如果投影面积总和比椭球面积小超过2%说明数据跨了两个分带需要改用按建筑所在地中央经线动态投影而不是整套数据固定4547。重庆面积不小主城和渝东北分属不同高斯带这时应该用geopandas的estimate_utm_crs或者按几何中心点动态选择投影。动态投影的一个快速实现对每个街道分组取组内所有几何的union中心点用geodesic反算经纬度再按3度带规则计算出中央经线构造一个Transverse Mercator投影。这个做法比固定EPSG更稳健但要保证同一街道不会被拆到两个带建议按乡镇级行政面做分区而不是按单栋建筑。4. 把人口信息匹配到建筑面上空间连接与面积权重分摊4.1 先分清人口字段的粒度别把社区人口直接塞给每一栋楼重庆建筑shp里的人口属性最理想的情况是每一栋建筑都有“同栋人口”字段但现实是很多数据源只有村/居委会一级的户籍人口甚至只有乡镇街道一级。这时你要有一个清晰的判断人口是“落在这些建筑里”的不是“属于这些建筑”的。如果你直接把街道人口均分到该街道每栋建筑上等于默认每栋楼住一样多的人这不符合城市建筑高度和容积率差异。正确思路是把人口总量从高层级向低层级“按面积权重”分摊。分摊公式某栋建筑分摊人口 街道总人口 × (该建筑建筑面积 / 该街道所有建筑建筑面积之和)这个公式有个隐含前提人均建筑面积在街道内基本一致。这个前提在城市不一定成立但在农村和城市分开做比混合在一起更靠谱。做之前先按城乡分组城、乡分别使用各自的分摊系数。4.2 用空间连接把社区人口挂到建筑面上第一步先把重庆的乡镇/街道边界shp读进来用空间连接给每栋建筑打上街道编码。常见做法是GeoPandas的sjoin把建筑面作为左表街道面作为右表操作类型是within或intersects。由于建筑边界可能贴着街道边界用intersects会有极少数建筑被匹配到两个街道用within又可能漏掉因边界重叠微差而未完全落入的建筑。稳妥方案是先用within然后对未匹配的建筑用representative_point()再算一次。import geopandas as gpd build gpd.read_file(cq_buildings_proj.shp) # 已投影修复 street gpd.read_file(cq_street.shp) # 街道边界 # 确保两边坐标系一致 street street.to_crs(build.crs) # 空间连接建筑在街道内 joined gpd.sjoin(build, street, howleft, predicatewithin) # 修正未匹配的把geometry换成内部代表点再算一次 cent build.copy() cent[geometry] cent.geometry.representative_point() joined_cent gpd.sjoin(cent, street, howleft, predicatewithin) # 用cent结果填充未匹配部分 joined[STREET_ID] joined[STREET_ID].fillna(joined_cent[STREET_ID])代码说明sjoin的predicatewithin要求左要素完全落在右要素内部处理建筑面靠谱但建筑几何可能因为采集误差稍微越界导致匹配失败。第二次用representative_point()取多边形内部的一个点再做within能修复大部分边界问题。STREET_ID是从街道属性来的主键如果不存在就自己构建street[temp_id] street.index然后join。参数说明howleft表示保留所有建筑面即使没有匹配到街道也会保留为NaN便于你检查漏匹predicate参数取决于GeoPandas版本旧版本用opwithin。代表点representative_point跟centroid不同centroid可能落在多边形外尤其对重庆这种带天井的传统院落式建筑所以别用centroid做匹配点。4.3 按建筑面积分摊人口的具体算法连接街道编码之后把街道人口字段合并进来然后按建筑面的area_m2_proj计算占比。注意这里要按街道内所有建筑的面积求和而不是按街道总面积求和。因为建筑面数据可能只覆盖房屋本体院落、棚房、空地在shp里不是建筑面用街道总面积会低估人口密度。# 假设joined已含街道总人口字段pop_street和街道IDSTREET_ID joined[area] joined[area_m2_proj] # 计算每栋建筑在所属街道内的面积占比 joined[area_share] joined.groupby(STREET_ID)[area].transform(lambda x: x / x.sum()) # 分摊人口结果保留两位小数 joined[pop_est] (joined[pop_street] * joined[area_share]).round(2) # 验证每个街道分摊人口总和应等于街道原总人口允许舍入差 check joined.groupby(STREET_ID).agg( 原人口(pop_street, first), 分摊合计(pop_est, sum) ) check[差值] check[原人口] - check[分摊合计] print(check[check[差值].abs() 0.5].head(10))这段代码先算面积占比再人口反推。transform(lambda x: x / x.sum())在groupby后对每个街道组做归一化结果与组原顺序一致不会错位。验证差值一般不会超过0.5人如果差值大可能是因为部分建筑area_m2_proj为0或NaN这些建筑会把分母污染导致占比失真。你要在计算前把面积NaN的行填为0或剔除并标记出来。注意如果shp里已有原始人口字段且缺失率不高那就用原始字段不要强行分摊。分摊只服务于没有建筑级人口的数据。如果街道人口只有户籍口径而你想拿到常住人口重庆的农村房屋普遍人户分离外出务工比例高直接分摊会高估实际居住密度。常见做法是根据第七次人口普查的镇街常住人口与户籍人口比值修正这一步就不是shp文件本身能解决的了需要外部统计表格配合。5. 落地技巧批量修复字段并导出在QGIS里用$area验证面积可信度面积和人口字段修复后别急着直接拿去做分析先用一个小脚本把成果固化成一个干净版本再用QGIS做一次独立验证。下面这个函数可以批量处理同名结构的shp文件把所有几何统一投影到EPSG:4547重算面积并输出精简属性表。def fix_and_export(in_file, out_file): gdf gpd.read_file(in_file, encodingutf-8) gdf gdf.to_crs(epsg4547) gdf[geometry] gdf.geometry.buffer(0) gdf[area_calc] gdf.geometry.area # 保留必要字段减少shp体积 keep [c for c in [OBJECTID, XZQMC, city_rural, pop_est, area_calc] if c in gdf.columns] gdf gdf[keep [geometry]] gdf.to_file(out_file, driverESRI Shapefile, encodingutf-8)压缩技巧shp中大量冗余字段会让dbf体积膨胀到几十MB。去掉不必要的浮点字段把双精度面积改成INT或保留两位小数的REAL能明显减小文件体积。如果只是做分析不下发直接转GeoPackage更合适ogr2ogr -f GPKG cq_buildings.gpkg cq_buildings.shp -nlt PROMOTE_TO_MULTI -lco ENCODINGUTF-8这条命令把shp转成GeoPackagePROMOTE_TO_MULTI保证所有MultiPolygon结构不变ENCODINGUTF-8防止中文属性乱码。GeoPackage只一个文件比shp一堆配套文件更利于归档。验证面积时在QGIS中打开导出的shp进入字段计算器新建字段d_area使用表达式$area。然后新建一个虚拟层用SQL比较字段与几何计算面积的差SELECT area_calc, $area AS qgis_area, abs(area_calc - $area) / $area AS diff_ratio FROM cq_buildings WHERE $area 0 AND abs(area_calc - $area) / $area 0.03如果查出来超过总数的1%说明投影选择仍有问题或者原始几何存在拓扑错误。再进一步可以把导出的建筑面跟道路shp叠加用“删除重复几何”检查是否有重叠面对重庆山地城市还要关注相邻建筑间的极窄缝隙buffer(0)之后用sym_difference检查成对边界。完全干净的数据在QGIS里应该能看到每个建筑面独立、无重叠并且area_calc与$area的差异控制在可接受范围内。本文还有配套的精品资源点击获取 重庆建筑面shp数据面积计算与人口匹配实战
简介以重庆市城市与农村建筑物面数据为核心的shp格式资源面要素内附建筑面积与人口等属性字段适用于内涝治理、SWMM建模、智慧城市评估与碳中和背景下的下垫面分析。数据整体约309MB共6个文件包含.shp空间几何、.dbf属性表、.prj坐标参考、.shx索引及配套元数据结构标准可直接加载至ArcGIS/QGIS。已有919人浏览学习适合城乡规划、市政水务和建筑能耗研究方向的初学者与工程师使用。借助该数据可批量提取屋顶轮廓与建筑密度为海绵城市模拟、片区人口估算及建筑能源模型提供基础底图省去自行采集与矢量化处理的环节。1. 重庆建筑物面数据城市农村shp不止是轮廓还带着面积和人口属性在GIS项目里拿到“重庆市建筑物面数据城市农村shp格式包含建筑面要素、面积、人口信息”这类数据时最常见的错误是把它当成普通建筑轮廓直接叠加底图画两笔专题图就完事。实际上这份数据里最有价值的是两列属性面积和人口。面积能帮你算容积率、户均面积、建造成本估算人口字段如果可靠可以整合进应急疏散、网格化治理、甚至“15分钟生活圈”补短板分析。它覆盖城市和农村意味着坐标系、精度、字段完整性在两类区域往往不一致不能一套参数走到底。这篇文章就从结构读取、面积修正、人口匹配到导出验证走完整条你应该对这份数据做的事。2. 读懂shp里的建筑面、面积字段和人口字段先解决字段混乱的问题2.1 用geopandas读出重庆建筑数据的基本结构与常见属性拿到shp第一件事不是打开QGIS看图层而是先看字段字典和空间范围。常见做法是用Python的geopandas读取打印字段名、前几行和坐标系用30秒判断这份数据有没有经过拆分或投影转换。import geopandas as gpd # 读取重庆建筑面shp假定路径为cq_buildings.shp gdf gpd.read_file(cq_buildings.shp, encodingutf-8) # 打印全部字段名、行数和坐标系 print(字段列表:, gdf.columns.tolist()) print(要素数量:, len(gdf)) print(坐标系:, gdf.crs) # 查看面积和人口相关字段的前5行判断字段实际内容 cols [c for c in gdf.columns if any(k in c.lower() for k in [area, mj, pop, rk, jzmj, hjs, hhs])] print(疑似面积/人口字段:, cols) print(gdf[cols].head())这段代码逻辑先用encoding参数应对图层属性表的乱码因为shp的dbf编码常见GBK或UTF-8再通过关键词筛出面积和人口字段避免自己一个个猜。注意read_file默认会读取全部几何大型文件可以先只读属性gdf gpd.read_file(path, ignore_geometryTrue)但那样就拿不到几何信息只适合快速查字段。参数说明encoding设错会导致中文乱码建议先试utf-8出现乱码再换gbkignore_geometryTrue时返回的GeoDataFrame没有geometry列后续无法做空间操作。实际项目中我一般先把属性读出来记录字段名后再用完整读取去做空间分析。2.2 面积字段的三种来源字段现成、几何计算、单位换算这份数据里面积可能来自三处一是数据生产方在dbf里直接写好的area/mj/jzmj字段二是AREA字段是ArcGIS自动生成但不带单位三是有字段但单位是平方千米而你要的是平方米。读数据后先做一致性检查用几何算一遍面积跟字段比差值。如果差值超过5%说明字段是原始等级路网时代的估算值不可直接用。常见的面积字段名称对照字段名抄录真实含义单位常见坑AREA图形面积平方米或平方千米ArcGIS导出时继承的旧值常单位错JZMJ建筑面积平方米可能与占地面积混淆CQMJ产权面积平方米多用于房产农村字段多为空POP/RK人口数人可能是户籍人口不是常住人口HJS户籍户数户农村常用城市缺失率高先别急着删除或覆盖字段建议新增一列area_geom用几何计算得出标准面积与原始area字段并存。这样后续分析可以随时回溯。2.3 人口信息字段的常见形态与缺失处理人口信息在这类shp里有两种存在方式一种是每个建筑面要素都挂一个人口数另一种是只有少数街区/村庄面要素带人口建筑面里是空的。拿重庆农村房屋来说很多数据源以“户数”代替“人口数”因为村民宅基地确权时登记的是户主与家庭成员数但统计口径五花八门有的含老人有的不含儿童。城市建筑则常见“居住人数估算值”由容积率乘以户均人数推出。如果人口字段全空你需要放弃“按建筑直接统计人口”的幻想改用一个靠谱的替代方案把重庆乡镇或街道边界的人口普查数据通过空间连接分摊到每个建筑面上。这个方案在第4章展开。目前这个阶段你要做的是把人口字段类型统一为float并记录缺失率# 统一人口字段名避免大小写和别名干扰 pop_cols [c for c in gdf.columns if c.upper() in (POP, POPULATION, RK, ZRK)] if pop_cols: gdf[pop_clean] gdf[pop_cols[0]].astype(float) else: gdf[pop_clean] None missing_rate gdf[pop_clean].isna().mean() print(f人口字段缺失率: {missing_rate:.1%}) # 缺失率超过20%就必须考虑外部人口数据做空间匹配这段代码先把疑似人口字段统一成float算出缺失率。缺失率超过20%就不建议用户直接按字段汇总了。参数说明astype(float)会把1,024这类带千分位字符串报错需要先gdf[pop_cols[0]] gdf[pop_cols[0]].astype(str).str.replace(,, )再转换。人口字段的清洗必须先于任何统计否则后面所有按人口汇总的结果都会带着隐性错误。3. 用正确的投影重算建筑面面积别被shp自带面积骗了3.1 为什么重庆地理坐标系下面积会偏小重庆地处东经105°17110°11、北纬28°1032°13如果shp的坐标系是WGS84或CGCS2000地理坐标系度为单位那么计算面积时任意GIS软件都会采用经纬度球面算法结果会是平方米级别的小数或者软件默认给一个虚假的“度”单位没人能看懂。更麻烦的是很多shp虽然显示为GCS_WGS_1984实际数据却来自高斯-克吕格投影的成果只是导出时丢了投影信息。正确的做法是把它投影到重庆当地常用的高斯-克吕格3度带投影。国家标准中重庆主城多用CGCS2000 / 3-degree Gauss-Kruger zone 35EPSG代码为4547也可用4548针对区县。有些团队用WGS 84 / UTM zone 48NEPSG:32648也是可行的但高斯投影同一条带内面积变形更小。3.2 用Python按高斯-克吕格投影计算椭球面积GeoPandas里最简单的做法是to_crs到EPSG:4547然后用area属性。area的值是该投影平面下的面积单位为平方米但它不是真正的椭球面积而是投影面积。对于50平方公里以下的县级尺度投影面积与椭球面积误差小于千分之一对重庆这种山地地形如果只做相对排名和比例分析完全够用。import geopandas as gpd gdf gpd.read_file(cq_buildings.shp, encodingutf-8) # 先统一几何有效性避免自相交和多边形环错误导致面积异常 gdf[geometry] gdf.geometry.buffer(0) # 投影到重庆高斯3度带CGCS2000, zone 35 gdf_proj gdf.to_crs(epsg4547) # 计算投影平面面积单位平方米 gdf_proj[area_m2_proj] gdf_proj.geometry.area # 再算椭球面积用pyproj的Geod from pyproj import Geod geod Geod(ellpsWGS84) def geodesic_area(geom): from shapely.geometry import Polygon, MultiPolygon if geom.is_empty: return 0.0 if geom.geom_type Polygon: x, y geom.exterior.xy area_m2, _ geod.polygon_area_perimeter(list(y), list(x)) return abs(area_m2) elif geom.geom_type MultiPolygon: total 0.0 for part in geom.geoms: total geodesic_area(part) return total else: return 0.0 gdf_proj[area_m2_geo] gdf_proj.apply(lambda r: geodesic_area(r[geometry]), axis1)这段代码有两组面积area_m2_proj是投影平面面积area_m2_geo是WGS84椭球下的测地线面积。因为坐标系是CGCS2000和WGS84在米级精度下差异极小直接用WGS84椭球参数没问题。注意geod.polygon_area_perimeter的输入是(纬, 经)的顺序也就是第一个参数是ylat第二个是xlon写反会得到荒谬值。buffer(0)这一步不是给几何加缓冲而是利用Shapely的容错机制修复细小的自相交和无效环使面积计算不会返回NaN。提示如果shp里有多部件建筑比如一排联排房屋被合并为一个MultiPolygon上述geodesic_area会把每个子多边形面积相加正合适。不要先dissolve否则总面积相同但会丢掉单栋建筑的独立性。# 如果存在无效几何对每个要素单独检查 invalid ~gdf_proj.geometry.is_valid print(f无效几何数量: {invalid.sum()}) if invalid.sum() 0: gdf_proj.loc[invalid, geometry] gdf_proj.loc[invalid, geometry].buffer(0)这段检查代码放在面积计算的任意阶段都可以主要确保area_m2_geo计算时不会因为单条自相交几何报错。注意buffer(0)会替换原几何如果后续还要做空间连接建议保留一个原始几何副本比如gdf_proj[geom_orig] gdf_proj.geometry因为buffer可能把极窄的建筑缝隙合并掉。3.3 按城市/农村分组汇总面积验证数据是否自洽拿到两种面积后顺手做一次分组汇总既验证数据完整性也给后续人口分摊提供分母。假如shp里带一个“城乡分类”字段如CLSD或XZQMC按它分组如果没有分类字段就按村级行政区划点做空间叠加来给每栋建筑打标。# 假设存在city_rural字段取值urban/rural if city_rural in gdf_proj.columns: df_group gdf_proj.groupby(city_rural).agg( 建筑数(geometry, count), 投影面积_m2(area_m2_proj, sum), 椭球面积_m2(area_m2_geo, sum), 人口_缺失率(pop_clean, lambda s: s.isna().mean()) ).reset_index() print(df_group) else: # 没有城乡字段时用土地利用或村庄点数据做空间连接 gdf_proj[city_rural] gdf_proj[XZQMC].str.contains(村|组, naFalse).map({True: rural, False: urban})这个聚合结果能告诉你三件事农村建筑数量是否远多于城市、人口字段缺失率是否在城乡间差异巨大、投影面积和椭球面积差值是否一致。如果投影面积总和比椭球面积小超过2%说明数据跨了两个分带需要改用按建筑所在地中央经线动态投影而不是整套数据固定4547。重庆面积不小主城和渝东北分属不同高斯带这时应该用geopandas的estimate_utm_crs或者按几何中心点动态选择投影。动态投影的一个快速实现对每个街道分组取组内所有几何的union中心点用geodesic反算经纬度再按3度带规则计算出中央经线构造一个Transverse Mercator投影。这个做法比固定EPSG更稳健但要保证同一街道不会被拆到两个带建议按乡镇级行政面做分区而不是按单栋建筑。4. 把人口信息匹配到建筑面上空间连接与面积权重分摊4.1 先分清人口字段的粒度别把社区人口直接塞给每一栋楼重庆建筑shp里的人口属性最理想的情况是每一栋建筑都有“同栋人口”字段但现实是很多数据源只有村/居委会一级的户籍人口甚至只有乡镇街道一级。这时你要有一个清晰的判断人口是“落在这些建筑里”的不是“属于这些建筑”的。如果你直接把街道人口均分到该街道每栋建筑上等于默认每栋楼住一样多的人这不符合城市建筑高度和容积率差异。正确思路是把人口总量从高层级向低层级“按面积权重”分摊。分摊公式某栋建筑分摊人口 街道总人口 × (该建筑建筑面积 / 该街道所有建筑建筑面积之和)这个公式有个隐含前提人均建筑面积在街道内基本一致。这个前提在城市不一定成立但在农村和城市分开做比混合在一起更靠谱。做之前先按城乡分组城、乡分别使用各自的分摊系数。4.2 用空间连接把社区人口挂到建筑面上第一步先把重庆的乡镇/街道边界shp读进来用空间连接给每栋建筑打上街道编码。常见做法是GeoPandas的sjoin把建筑面作为左表街道面作为右表操作类型是within或intersects。由于建筑边界可能贴着街道边界用intersects会有极少数建筑被匹配到两个街道用within又可能漏掉因边界重叠微差而未完全落入的建筑。稳妥方案是先用within然后对未匹配的建筑用representative_point()再算一次。import geopandas as gpd build gpd.read_file(cq_buildings_proj.shp) # 已投影修复 street gpd.read_file(cq_street.shp) # 街道边界 # 确保两边坐标系一致 street street.to_crs(build.crs) # 空间连接建筑在街道内 joined gpd.sjoin(build, street, howleft, predicatewithin) # 修正未匹配的把geometry换成内部代表点再算一次 cent build.copy() cent[geometry] cent.geometry.representative_point() joined_cent gpd.sjoin(cent, street, howleft, predicatewithin) # 用cent结果填充未匹配部分 joined[STREET_ID] joined[STREET_ID].fillna(joined_cent[STREET_ID])代码说明sjoin的predicatewithin要求左要素完全落在右要素内部处理建筑面靠谱但建筑几何可能因为采集误差稍微越界导致匹配失败。第二次用representative_point()取多边形内部的一个点再做within能修复大部分边界问题。STREET_ID是从街道属性来的主键如果不存在就自己构建street[temp_id] street.index然后join。参数说明howleft表示保留所有建筑面即使没有匹配到街道也会保留为NaN便于你检查漏匹predicate参数取决于GeoPandas版本旧版本用opwithin。代表点representative_point跟centroid不同centroid可能落在多边形外尤其对重庆这种带天井的传统院落式建筑所以别用centroid做匹配点。4.3 按建筑面积分摊人口的具体算法连接街道编码之后把街道人口字段合并进来然后按建筑面的area_m2_proj计算占比。注意这里要按街道内所有建筑的面积求和而不是按街道总面积求和。因为建筑面数据可能只覆盖房屋本体院落、棚房、空地在shp里不是建筑面用街道总面积会低估人口密度。# 假设joined已含街道总人口字段pop_street和街道IDSTREET_ID joined[area] joined[area_m2_proj] # 计算每栋建筑在所属街道内的面积占比 joined[area_share] joined.groupby(STREET_ID)[area].transform(lambda x: x / x.sum()) # 分摊人口结果保留两位小数 joined[pop_est] (joined[pop_street] * joined[area_share]).round(2) # 验证每个街道分摊人口总和应等于街道原总人口允许舍入差 check joined.groupby(STREET_ID).agg( 原人口(pop_street, first), 分摊合计(pop_est, sum) ) check[差值] check[原人口] - check[分摊合计] print(check[check[差值].abs() 0.5].head(10))这段代码先算面积占比再人口反推。transform(lambda x: x / x.sum())在groupby后对每个街道组做归一化结果与组原顺序一致不会错位。验证差值一般不会超过0.5人如果差值大可能是因为部分建筑area_m2_proj为0或NaN这些建筑会把分母污染导致占比失真。你要在计算前把面积NaN的行填为0或剔除并标记出来。注意如果shp里已有原始人口字段且缺失率不高那就用原始字段不要强行分摊。分摊只服务于没有建筑级人口的数据。如果街道人口只有户籍口径而你想拿到常住人口重庆的农村房屋普遍人户分离外出务工比例高直接分摊会高估实际居住密度。常见做法是根据第七次人口普查的镇街常住人口与户籍人口比值修正这一步就不是shp文件本身能解决的了需要外部统计表格配合。5. 落地技巧批量修复字段并导出在QGIS里用$area验证面积可信度面积和人口字段修复后别急着直接拿去做分析先用一个小脚本把成果固化成一个干净版本再用QGIS做一次独立验证。下面这个函数可以批量处理同名结构的shp文件把所有几何统一投影到EPSG:4547重算面积并输出精简属性表。def fix_and_export(in_file, out_file): gdf gpd.read_file(in_file, encodingutf-8) gdf gdf.to_crs(epsg4547) gdf[geometry] gdf.geometry.buffer(0) gdf[area_calc] gdf.geometry.area # 保留必要字段减少shp体积 keep [c for c in [OBJECTID, XZQMC, city_rural, pop_est, area_calc] if c in gdf.columns] gdf gdf[keep [geometry]] gdf.to_file(out_file, driverESRI Shapefile, encodingutf-8)压缩技巧shp中大量冗余字段会让dbf体积膨胀到几十MB。去掉不必要的浮点字段把双精度面积改成INT或保留两位小数的REAL能明显减小文件体积。如果只是做分析不下发直接转GeoPackage更合适ogr2ogr -f GPKG cq_buildings.gpkg cq_buildings.shp -nlt PROMOTE_TO_MULTI -lco ENCODINGUTF-8这条命令把shp转成GeoPackagePROMOTE_TO_MULTI保证所有MultiPolygon结构不变ENCODINGUTF-8防止中文属性乱码。GeoPackage只一个文件比shp一堆配套文件更利于归档。验证面积时在QGIS中打开导出的shp进入字段计算器新建字段d_area使用表达式$area。然后新建一个虚拟层用SQL比较字段与几何计算面积的差SELECT area_calc, $area AS qgis_area, abs(area_calc - $area) / $area AS diff_ratio FROM cq_buildings WHERE $area 0 AND abs(area_calc - $area) / $area 0.03如果查出来超过总数的1%说明投影选择仍有问题或者原始几何存在拓扑错误。再进一步可以把导出的建筑面跟道路shp叠加用“删除重复几何”检查是否有重叠面对重庆山地城市还要关注相邻建筑间的极窄缝隙buffer(0)之后用sym_difference检查成对边界。完全干净的数据在QGIS里应该能看到每个建筑面独立、无重叠并且area_calc与$area的差异控制在可接受范围内。本文还有配套的精品资源点击获取