尧图网站建设 
简介2022年沈阳建筑轮廓GIS数据是一份面向城市规划者、建筑师和GIS分析人员的高质量矢量数据集。资源以标准Shapefile格式存储共6个文件包含建筑轮廓几何shp、属性信息dbf、空间索引shx、投影坐标prj及元数据xml等压缩包大小82.05MB。数据覆盖沈阳市主要建筑轮廓记录建筑位置、形状、高度、面积及用途等关键属性可直接在ArcGIS或QGIS中加载用于建筑密度分析、土地利用率评估、城市更新规划及日照模拟等应用场景。同时配套的元数据文件记录了采集方法与精度等关键信息为开展城市热岛效应研究、灾害风险评估等进阶分析奠定数据基础。目前已有288人学习下载对于希望获取真实城市建筑空间数据进行实操练习的GIS学习者与研究人员而言是一份难得的城市级样本数据。1. 一份建筑轮廓数据背后是一整条 GIS 处理链很多人在做城市分析时第一步就卡在数据上。想要 2022 年沈阳的建筑轮廓找到的要么是 OSM 导出的稀疏面要么是某个规划院出的局部成果再要么是遥感解译后精度不明的瓦片。真正的问题不是“哪里能下载”而是拿到一份轮廓数据后怎么确认它的坐标系、清洗掉异常几何、把散落的点坐标按规则成面再落到街道、缓冲区、路网上做分析最后批量出图和标注。这篇博客直接围绕“2022 年沈阳建筑轮廓 GIS 数据”展开先讲清楚一份轮廓数据在 GIS 里的结构再用 GDAL 和 GeoPandas 把数据读进来、清洗、转坐标、做叠加分析最后给出 QGIS 的批量出图和双标注配置。适合刚接触 GIS 数据处理的工程师也适合有基础但没系统性处理过建筑轮廓数据的人。2. 先读懂轮廓数据几何、字段和藏在里面的坐标系2.1 建筑轮廓的三种常见形态建筑轮廓数据在 GIS 里本质是一组闭合多边形。每个多边形代表一栋建筑在水平面上的投影而投影形状、节点密度、属性完整度决定了这份数据能做什么分析。常见的格式有三种。Shapefile 是旧但最普及的一套数据由 .shp、.dbf、.shx、.prj 等至少四个文件构成缺了 .prj 坐标系就丢了。GeoJSON 是 Web 和开源工具的首选结构直观一个文件搞定几何和属性但体积大、读取慢。FileGDB 是 ArcGIS 的原生格式支持复杂拓扑和更多数据类型适合企业级数据管理。另外要区分“面数据”和“线数据”。建筑轮廓的主数据是 Polygon 或 MultiPolygon但某些来源比如 CAD 转出给的是闭合 Polyline需要自己首尾相连生成面。不要一拿到线数据就直接做分析先确认几何类型是基本纪律。2.2 坐标系决定一切尤其是沈阳这种高纬度城市这是这份数据处理里最容易被忽视的一步。沈阳位于东经约 123 度、北纬约 42 度这意味着同样一套坐标放在不同坐标系里面积计算结果差异是巨大的。先看三个最常遇到的坐标系坐标系常见来源适合做什么不适合做什么WGS84 / CGCS2000 经纬度OSM 导出、GPS 采集、部分公开 API数据交换、在线地图显示面积和距离计算单位是度CGCS2000 3度带 GK中央经线 123E带号 41国内测绘成果、规划院数据面积、距离、缓冲区分析直接叠加互联网底图WGS84 / UTM 51NEPSG 32651国际数据格式、部分影像解译成果面积计算、空间分析与国内坐标成果混用其中重点是如果数据是经纬度直接buffer(30)得到的是 30 度而不是 30 米。而如果数据是 Web 墨卡托EPSG 3857在沈阳纬度做面积统计结果会被放大到真实面积的约 1.35 倍。很多分析报告里的“覆盖率”虚高根源就是没做投影。2.3 用 GDAL 把沈阳建筑轮廓数据读出来最常见的做法是先用 GDAL 命令行快速摸清数据底细再进入 Python 做处理。ogrinfo -so shenyang_building_2022.shp shenyang_building_2022这条命令不读取几何内部节点只输出图层级信息要素数量、几何类型、字段列表和空间参考。如果输出里空间参考为空或者在打开 .dbf 时中文属性字段乱码说明 .prj 缺失或字符集不对需要先修复这两项再进入后续步骤。用 Python 读取并检查几何类型代码如下from osgeo import ogr ds ogr.Open(shenyang_building_2022.shp, 0) layer ds.GetLayer(0) print(要素数量:, layer.GetFeatureCount()) print(空间参考:, layer.GetSpatialRef().ExportToWkt()[:200]) # 输出字段名与类型 defn layer.GetLayerDefn() for i in range(defn.GetFieldCount()): fdefn defn.GetFieldDefn(i) print(字段:, fdefn.GetName(), fdefn.GetTypeName()) # 检查第一个要素的几何 feat layer.GetNextFeature() geom feat.GetGeometryRef() if geom.GetGeometryName() MULTIPOLYGON: point_count sum(part.GetPointCount() for part in geom) else: point_count geom.GetPointCount() print(几何类型:, geom.GetGeometryName(), 节点数:, point_count)这段代码通过GetSpatialRef().ExportToWkt()拿到完整的坐标系定义文本如果该文本描述的是经纬度单位直接判断为“未投影”状态通过遍历字段名可以提前发现结构变化比如 2022 年的数据是否新增了SOURCE或CAPTURE_DATE字段避免后续按老字段名写死逻辑。3. 清洗与坐标转换从撞墙到能用的轮廓面3.1 几何修复自相交、重复点和 MultiPolygon 假象建筑轮廓数据里最普遍的三个脏问题是自相交、重复节点和退化面。自相交的多边形在计算面积时会得到错误结果甚至直接报错重复节点不影响视觉但影响性能退化面则通常是捕捉错误留下的零面积碎片。处理顺序建议是先查退化面再修自相交最后统一几何类型。用 Shapely 的make_valid可以一次性解决绝大多数几何合法性问题from shapely.validation import make_valid from shapely.geometry import shape import geopandas as gpd gdf gpd.read_file(shenyang_building_2022.shp, encodingutf-8) # 统计初始无效几何 invalid gdf[~gdf.geometry.is_valid] print(无效几何数量:, len(invalid)) # 修复并只保留面类型 gdf[geometry] gdf.geometry.apply(make_valid) gdf gdf[gdf.geometry.geom_type.isin([Polygon, MultiPolygon])] gdf gdf[gdf.geometry.area 1.0] # 去掉 1 平方米以下的碎片这里用了gpd.read_file(..., encodingutf-8)。如果 .dbf 是 GBK 编码需要把encoding改成gbk否则中文字段全变乱码。判断方法是读出来后打印前 5 行的NAME字段看到无法辨识的中文就换编码再试。make_valid会产生几种结果对自相交面它会拆成多个合法的面对某些复杂情况它会返回 GeometryCollection里面可能混着线和点。所以修复后必须过滤geom_type把非面要素剔除否则后续intersects、overlay等操作会持续报错。3.2 GIS 坐标转点质心提取与投影统一把建筑轮廓转成点是空间分析的常规动作。点要素可以用于密度分析、最近邻匹配、图表展示也能规避大量多边形叠加带来的性能开销。转点前必须确定目标坐标系。如果最终要算距离和面积我一般转成 UTM 51NEPSG 32651单位是米误差可控。用 pyproj 完成转换并提取质心import pyproj from shapely.ops import transform from shapely.geometry import Point # 从 WGS84 经纬度转到 UTM 51N wgs84 pyproj.CRS(EPSG:4326) utm51 pyproj.CRS(EPSG:32651) project pyproj.Transformer.from_crs(wgs84, utm51, always_xyTrue).transform gdf_utm gdf.to_crs(epsg32651) gdf_utm[centroid] gdf_utm.geometry.centroid # 提取成独立的点图层 points gpd.GeoDataFrame( gdf_utm.drop(columnsgeometry), geometrygdf_utm[centroid], crsEPSG:32651 ) points.to_file(shenyang_building_points.shp, encodingutf-8)参数说明always_xyTrue确保坐标顺序固定为经度、纬度而不是纬度、经度这一点在批量处理时特别容易踩坑to_crs(epsg32651)把整个 GeoDataFrame 的坐标系从原始状态转到 UTM 51N后面的centroid计算全部基于米制单位buffer 出来的距离也才是真实的米。注意如果原始数据本身是 CGCS2000 高斯克吕格投影它的坐标值看起来像八位数的平面坐标这时不要强行转 4326 再转 32651直接用to_crs(epsg32651)即可pyproj 会依据内置的基准面转换参数完成换算。沈阳地区 CGCS2000 与 WGS84 的差异通常在亚米级对建筑轮廓分析完全可接受。3.3 GIS 坐标成面两种最常见的成面场景标题里带“轮廓”但很多实际项目拿到的并不是现成的面而是离散的坐标点集或闭合线。这种情况下GIS 坐标成面就是核心步骤。第一种场景由边界线生成面。CAD 图纸转出的建筑轮廓通常是闭合多段线导入 GIS 后是 LineString。用 GeoPandas 的polygonize生成面from shapely.ops import polygonize from shapely.geometry import LineString # lines 是从 CAD 导入的 LineString 列表 # 先保证线是闭合的未闭合的要补终点 closed_lines [] for line in lines: coords list(line.coords) if coords[0] ! coords[-1]: coords.append(coords[0]) closed_lines.append(LineString(coords)) polygons list(polygonize(closed_lines)) print(生成面数量:, len(polygons))第二种场景由质心反算建筑范围。有些精简版数据只有坐标点和楼层数要生成轮廓面只能做缓冲区但这会与实际建筑形状不符。稳妥做法是先用质心做小半径缓冲区用于统计和可视化不要声称它代表真实轮廓。# 注意单位为米必须在投影坐标下操作 points_utm points.to_crs(epsg32651) buffered points_utm.copy() buffered[geometry] points_utm.geometry.buffer(8.0)参数说明8 米是按单栋建筑约 16 米见方估算的半径仅用于示意图。如果数据里有WIDTH或LENGTH字段可以用这两个字段动态生成矩形面比固定缓冲区更接近真实轮廓。具体做法是用方位角和半宽半高计算四个角点再构造成 Polygon但前提是数据本身含朝向字段否则只能按正南正北近似。3.4 常见错误对照为什么不能直接干这些事很多人会把点数据直接to_file成 Shapefile 拿去成图结果 ArcGIS 或 QGIS 打开是乱的。问题往往集中在三点坐标没统一。原始坐标是经纬度但目标图层是高斯投影叠加时错位几公里。投影有值但基准面不对。WGS84 与 CGCS2000 在沈阳的偏移通常在几十厘米到两三米对跨图层叠加影响不大但与影像配准时会出现边缘错位。编译环境问题。Windows 下 Python 读 Shapefile 时遇到中文路径会直接失败先把文件复制到纯英文路径再处理。提示处理前先对原始文件做一次备份。几何修复是不可逆操作make_valid会把原几何拆解重组保底文件可以避免修复结果不理想时无路可退。4. 空间分析覆盖率、体量与路网叠加4.1 按行政区统计建筑数量与覆盖率这一节把清洗好的轮廓数据用到真正的业务分析里。最常见的需求是按沈阳各区统计建筑数量和基底面积再除以区域面积得到建筑覆盖率。用 GeoPandas 的空间连接完成分组统计import geopandas as gpd buildings gpd.read_file(shenyang_building_2022_clean.shp, encodingutf-8).to_crs(epsg32651) districts gpd.read_file(shenyang_district.shp, encodingutf-8).to_crs(epsg32651) # 先投影再连接避免单位不一致导致面积统计失真 joined gpd.sjoin(buildings, districts, howinner, predicatewithin) districts[building_count] 0 districts[building_area] 0.0 count_series joined.groupby(index_right).size() area_series joined.groupby(index_right).apply( lambda g: g.geometry.area.sum(), include_groupsFalse ) districts.loc[count_series.index, building_count] count_series.values districts.loc[area_series.index, building_area] area_series.values districts[coverage_rate] districts[building_area] / districts.geometry.area print(districts[[NAME, building_count, coverage_rate]].head())参数说明predicatewithin表示建筑整体必须落在行政区边界内如果建筑横跨两个区within会将其剔除此时应改用intersects并按面积占比拆分但沈阳实际行政区边界与建筑轮廓误交叉的情况较少within的计算速度更快。include_groupsFalse是较新版本 GeoPandas 对apply行为的修正老版本会因分组列参与计算而报错。4.2 用层高字段推算总建筑面积建筑轮廓只有基底面积但很多业务指标需要总建筑面积。常见做法是用FLOORS层数字段乘以基底面积。# 清洗层数字段缺失或异常值处理 buildings[floors] pd.to_numeric(buildings[FLOORS], errorscoerce) buildings[floors] buildings[floors].fillna(1) buildings[floors] buildings[floors].clip(lower1, upper60) buildings[total_area] buildings.geometry.area * buildings[floors] print(buildings[total_area].sum())这里有个关键细节FLOORS字段如果是字符串类型to_numeric会把无法解析的值转成 NaN填充为 1 意味着默认按单层处理。对于明显偏离真实值的层数用clip限制上下限防止几栋异常高楼把全市统计结果拉偏。如果数据里只有HEIGHT而没有FLOORS可以按每层 3 米折算但沈阳老旧小区层高在 2.8 米到 3.1 米之间浮动折算误差无法避免。更可靠的做法是抽样比对随机挑 50 栋建筑用影像测量实际层数再和字段值对比确定该数据集的折算系数。4.3 与路网、POI 的叠加建筑到道路的距离与配套可达性建筑轮廓与路网的叠加能回答“这栋楼离最近道路多远”这类问题。计算最近道路距离的通用做法是对每栋建筑质心调用空间索引找最近的路段。from shapely.strtree import STRtree roads gpd.read_file(shenyang_roads.shp, encodingutf-8).to_crs(epsg32651) # 只保留主干道以上等级避免街坊路干扰计算结果 roads roads[roads[class].isin([trunk, primary, secondary])] tree STRtree(roads.geometry.values) distances [] for g in buildings.geometry.centroid: nearest tree.nearest(g) distances.append(g.distance(nearest)) buildings[dist_to_road] distances参数说明STRtree在构建时对所有道路几何建立空间索引nearest方法直接返回最近几何的索引复杂度接近 O(log n)几万栋建筑在几秒内完成计算。过滤道路等级这一步很重要如果把小区内部路和步行道混进去距离结果会系统性偏小尤其沈阳老城区胡同密集不筛选会影响“临街建筑”判定。到这里一份建筑轮廓数据已经经历了“读取确认 → 清洗修复 → 坐标统一 → 转点成面 → 空间统计 → 路网距离”的完整链路。最后一步是把这些结果变成能交给业务方的图件。5. 出图与双标注用 QGIS 把沈阳轮廓数据做成成果图5.1 批量出图按区分幅打印沈阳中心城区的建筑轮廓密度大一次导出全市范围的图几乎无法看清细节。常见做法是按行政区或图幅批量出图QGIS 的 Atlas 功能可以完成。在打印布局中新建布局添加地图项并设置为Atlas模式覆盖图层选择区分界数据。然后在布局里插入标签内容使用 Atlas 表达式引用区名和统计指标例如[% NAME %] 区建筑覆盖率 [% round(coverage_rate * 100, 1) %]%设置完成后导出 PDF 或 PNG 时会按每个区自动生成一页。批量导出 10 个区单次完成格式统一。5.2 GIS 怎么显示两个标注一个表达式解决单独的建筑数据往往既要显示编号又要显示楼层数。QGIS 的标注系统默认只允许一个标注样式但可以用表达式把两个属性拼接在一行concat(FID, \n, FLOORS, F)这样标注会显示成两行第一行是建筑编号第二行是楼层数加 F 后缀。如果希望两个标注分别用不同颜色或字号可以启用“标注工具箱”里的多标注支持或直接在表达式里用 HTML 标签控制颜色concat(font color#2b2b2b, FID, /font, \n, font color#d43f3a, FLOORS, F, /font)在随后的地图导出时勾选“允许 HTML 格式”才能让字号和颜色生效。相比做两个标注图层叠加这种方式在图例和导出时都更干净。5.3 验证把轮廓画到影像底图上检查出图前最后一道工序是位置验证。做法很简单加载天地图或 Esri World Imagery 底图把建筑轮廓图层设置为半透明红色缩放至 1:2000 以下逐片区目检。重点分辨两类错位一是整个图层整体偏移数米说明坐标系或基准面不匹配二是局部区域偏移说明原始数据本身存在分幅接边问题。对于后者不要指望统一处理能解决只能定位到具体区域后用 QGIS 的“几何编辑器”手工调整。最后检查单栋建筑轮廓是否与影像上的屋顶阴影边界重合。太阳方位造成的阴影会让影像里的建筑边缘比真实轮廓偏移不能拿阴影边缘当标准应该以墙根线为准。这一步做完这份 2022 年沈阳建筑轮廓数据才真正具备对外发布和支撑业务分析的资格。本文还有配套的精品资源点击获取 2022沈阳建筑轮廓GIS数据处理：从清洗修复到空间分析出图全流程
简介2022年沈阳建筑轮廓GIS数据是一份面向城市规划者、建筑师和GIS分析人员的高质量矢量数据集。资源以标准Shapefile格式存储共6个文件包含建筑轮廓几何shp、属性信息dbf、空间索引shx、投影坐标prj及元数据xml等压缩包大小82.05MB。数据覆盖沈阳市主要建筑轮廓记录建筑位置、形状、高度、面积及用途等关键属性可直接在ArcGIS或QGIS中加载用于建筑密度分析、土地利用率评估、城市更新规划及日照模拟等应用场景。同时配套的元数据文件记录了采集方法与精度等关键信息为开展城市热岛效应研究、灾害风险评估等进阶分析奠定数据基础。目前已有288人学习下载对于希望获取真实城市建筑空间数据进行实操练习的GIS学习者与研究人员而言是一份难得的城市级样本数据。1. 一份建筑轮廓数据背后是一整条 GIS 处理链很多人在做城市分析时第一步就卡在数据上。想要 2022 年沈阳的建筑轮廓找到的要么是 OSM 导出的稀疏面要么是某个规划院出的局部成果再要么是遥感解译后精度不明的瓦片。真正的问题不是“哪里能下载”而是拿到一份轮廓数据后怎么确认它的坐标系、清洗掉异常几何、把散落的点坐标按规则成面再落到街道、缓冲区、路网上做分析最后批量出图和标注。这篇博客直接围绕“2022 年沈阳建筑轮廓 GIS 数据”展开先讲清楚一份轮廓数据在 GIS 里的结构再用 GDAL 和 GeoPandas 把数据读进来、清洗、转坐标、做叠加分析最后给出 QGIS 的批量出图和双标注配置。适合刚接触 GIS 数据处理的工程师也适合有基础但没系统性处理过建筑轮廓数据的人。2. 先读懂轮廓数据几何、字段和藏在里面的坐标系2.1 建筑轮廓的三种常见形态建筑轮廓数据在 GIS 里本质是一组闭合多边形。每个多边形代表一栋建筑在水平面上的投影而投影形状、节点密度、属性完整度决定了这份数据能做什么分析。常见的格式有三种。Shapefile 是旧但最普及的一套数据由 .shp、.dbf、.shx、.prj 等至少四个文件构成缺了 .prj 坐标系就丢了。GeoJSON 是 Web 和开源工具的首选结构直观一个文件搞定几何和属性但体积大、读取慢。FileGDB 是 ArcGIS 的原生格式支持复杂拓扑和更多数据类型适合企业级数据管理。另外要区分“面数据”和“线数据”。建筑轮廓的主数据是 Polygon 或 MultiPolygon但某些来源比如 CAD 转出给的是闭合 Polyline需要自己首尾相连生成面。不要一拿到线数据就直接做分析先确认几何类型是基本纪律。2.2 坐标系决定一切尤其是沈阳这种高纬度城市这是这份数据处理里最容易被忽视的一步。沈阳位于东经约 123 度、北纬约 42 度这意味着同样一套坐标放在不同坐标系里面积计算结果差异是巨大的。先看三个最常遇到的坐标系坐标系常见来源适合做什么不适合做什么WGS84 / CGCS2000 经纬度OSM 导出、GPS 采集、部分公开 API数据交换、在线地图显示面积和距离计算单位是度CGCS2000 3度带 GK中央经线 123E带号 41国内测绘成果、规划院数据面积、距离、缓冲区分析直接叠加互联网底图WGS84 / UTM 51NEPSG 32651国际数据格式、部分影像解译成果面积计算、空间分析与国内坐标成果混用其中重点是如果数据是经纬度直接buffer(30)得到的是 30 度而不是 30 米。而如果数据是 Web 墨卡托EPSG 3857在沈阳纬度做面积统计结果会被放大到真实面积的约 1.35 倍。很多分析报告里的“覆盖率”虚高根源就是没做投影。2.3 用 GDAL 把沈阳建筑轮廓数据读出来最常见的做法是先用 GDAL 命令行快速摸清数据底细再进入 Python 做处理。ogrinfo -so shenyang_building_2022.shp shenyang_building_2022这条命令不读取几何内部节点只输出图层级信息要素数量、几何类型、字段列表和空间参考。如果输出里空间参考为空或者在打开 .dbf 时中文属性字段乱码说明 .prj 缺失或字符集不对需要先修复这两项再进入后续步骤。用 Python 读取并检查几何类型代码如下from osgeo import ogr ds ogr.Open(shenyang_building_2022.shp, 0) layer ds.GetLayer(0) print(要素数量:, layer.GetFeatureCount()) print(空间参考:, layer.GetSpatialRef().ExportToWkt()[:200]) # 输出字段名与类型 defn layer.GetLayerDefn() for i in range(defn.GetFieldCount()): fdefn defn.GetFieldDefn(i) print(字段:, fdefn.GetName(), fdefn.GetTypeName()) # 检查第一个要素的几何 feat layer.GetNextFeature() geom feat.GetGeometryRef() if geom.GetGeometryName() MULTIPOLYGON: point_count sum(part.GetPointCount() for part in geom) else: point_count geom.GetPointCount() print(几何类型:, geom.GetGeometryName(), 节点数:, point_count)这段代码通过GetSpatialRef().ExportToWkt()拿到完整的坐标系定义文本如果该文本描述的是经纬度单位直接判断为“未投影”状态通过遍历字段名可以提前发现结构变化比如 2022 年的数据是否新增了SOURCE或CAPTURE_DATE字段避免后续按老字段名写死逻辑。3. 清洗与坐标转换从撞墙到能用的轮廓面3.1 几何修复自相交、重复点和 MultiPolygon 假象建筑轮廓数据里最普遍的三个脏问题是自相交、重复节点和退化面。自相交的多边形在计算面积时会得到错误结果甚至直接报错重复节点不影响视觉但影响性能退化面则通常是捕捉错误留下的零面积碎片。处理顺序建议是先查退化面再修自相交最后统一几何类型。用 Shapely 的make_valid可以一次性解决绝大多数几何合法性问题from shapely.validation import make_valid from shapely.geometry import shape import geopandas as gpd gdf gpd.read_file(shenyang_building_2022.shp, encodingutf-8) # 统计初始无效几何 invalid gdf[~gdf.geometry.is_valid] print(无效几何数量:, len(invalid)) # 修复并只保留面类型 gdf[geometry] gdf.geometry.apply(make_valid) gdf gdf[gdf.geometry.geom_type.isin([Polygon, MultiPolygon])] gdf gdf[gdf.geometry.area 1.0] # 去掉 1 平方米以下的碎片这里用了gpd.read_file(..., encodingutf-8)。如果 .dbf 是 GBK 编码需要把encoding改成gbk否则中文字段全变乱码。判断方法是读出来后打印前 5 行的NAME字段看到无法辨识的中文就换编码再试。make_valid会产生几种结果对自相交面它会拆成多个合法的面对某些复杂情况它会返回 GeometryCollection里面可能混着线和点。所以修复后必须过滤geom_type把非面要素剔除否则后续intersects、overlay等操作会持续报错。3.2 GIS 坐标转点质心提取与投影统一把建筑轮廓转成点是空间分析的常规动作。点要素可以用于密度分析、最近邻匹配、图表展示也能规避大量多边形叠加带来的性能开销。转点前必须确定目标坐标系。如果最终要算距离和面积我一般转成 UTM 51NEPSG 32651单位是米误差可控。用 pyproj 完成转换并提取质心import pyproj from shapely.ops import transform from shapely.geometry import Point # 从 WGS84 经纬度转到 UTM 51N wgs84 pyproj.CRS(EPSG:4326) utm51 pyproj.CRS(EPSG:32651) project pyproj.Transformer.from_crs(wgs84, utm51, always_xyTrue).transform gdf_utm gdf.to_crs(epsg32651) gdf_utm[centroid] gdf_utm.geometry.centroid # 提取成独立的点图层 points gpd.GeoDataFrame( gdf_utm.drop(columnsgeometry), geometrygdf_utm[centroid], crsEPSG:32651 ) points.to_file(shenyang_building_points.shp, encodingutf-8)参数说明always_xyTrue确保坐标顺序固定为经度、纬度而不是纬度、经度这一点在批量处理时特别容易踩坑to_crs(epsg32651)把整个 GeoDataFrame 的坐标系从原始状态转到 UTM 51N后面的centroid计算全部基于米制单位buffer 出来的距离也才是真实的米。注意如果原始数据本身是 CGCS2000 高斯克吕格投影它的坐标值看起来像八位数的平面坐标这时不要强行转 4326 再转 32651直接用to_crs(epsg32651)即可pyproj 会依据内置的基准面转换参数完成换算。沈阳地区 CGCS2000 与 WGS84 的差异通常在亚米级对建筑轮廓分析完全可接受。3.3 GIS 坐标成面两种最常见的成面场景标题里带“轮廓”但很多实际项目拿到的并不是现成的面而是离散的坐标点集或闭合线。这种情况下GIS 坐标成面就是核心步骤。第一种场景由边界线生成面。CAD 图纸转出的建筑轮廓通常是闭合多段线导入 GIS 后是 LineString。用 GeoPandas 的polygonize生成面from shapely.ops import polygonize from shapely.geometry import LineString # lines 是从 CAD 导入的 LineString 列表 # 先保证线是闭合的未闭合的要补终点 closed_lines [] for line in lines: coords list(line.coords) if coords[0] ! coords[-1]: coords.append(coords[0]) closed_lines.append(LineString(coords)) polygons list(polygonize(closed_lines)) print(生成面数量:, len(polygons))第二种场景由质心反算建筑范围。有些精简版数据只有坐标点和楼层数要生成轮廓面只能做缓冲区但这会与实际建筑形状不符。稳妥做法是先用质心做小半径缓冲区用于统计和可视化不要声称它代表真实轮廓。# 注意单位为米必须在投影坐标下操作 points_utm points.to_crs(epsg32651) buffered points_utm.copy() buffered[geometry] points_utm.geometry.buffer(8.0)参数说明8 米是按单栋建筑约 16 米见方估算的半径仅用于示意图。如果数据里有WIDTH或LENGTH字段可以用这两个字段动态生成矩形面比固定缓冲区更接近真实轮廓。具体做法是用方位角和半宽半高计算四个角点再构造成 Polygon但前提是数据本身含朝向字段否则只能按正南正北近似。3.4 常见错误对照为什么不能直接干这些事很多人会把点数据直接to_file成 Shapefile 拿去成图结果 ArcGIS 或 QGIS 打开是乱的。问题往往集中在三点坐标没统一。原始坐标是经纬度但目标图层是高斯投影叠加时错位几公里。投影有值但基准面不对。WGS84 与 CGCS2000 在沈阳的偏移通常在几十厘米到两三米对跨图层叠加影响不大但与影像配准时会出现边缘错位。编译环境问题。Windows 下 Python 读 Shapefile 时遇到中文路径会直接失败先把文件复制到纯英文路径再处理。提示处理前先对原始文件做一次备份。几何修复是不可逆操作make_valid会把原几何拆解重组保底文件可以避免修复结果不理想时无路可退。4. 空间分析覆盖率、体量与路网叠加4.1 按行政区统计建筑数量与覆盖率这一节把清洗好的轮廓数据用到真正的业务分析里。最常见的需求是按沈阳各区统计建筑数量和基底面积再除以区域面积得到建筑覆盖率。用 GeoPandas 的空间连接完成分组统计import geopandas as gpd buildings gpd.read_file(shenyang_building_2022_clean.shp, encodingutf-8).to_crs(epsg32651) districts gpd.read_file(shenyang_district.shp, encodingutf-8).to_crs(epsg32651) # 先投影再连接避免单位不一致导致面积统计失真 joined gpd.sjoin(buildings, districts, howinner, predicatewithin) districts[building_count] 0 districts[building_area] 0.0 count_series joined.groupby(index_right).size() area_series joined.groupby(index_right).apply( lambda g: g.geometry.area.sum(), include_groupsFalse ) districts.loc[count_series.index, building_count] count_series.values districts.loc[area_series.index, building_area] area_series.values districts[coverage_rate] districts[building_area] / districts.geometry.area print(districts[[NAME, building_count, coverage_rate]].head())参数说明predicatewithin表示建筑整体必须落在行政区边界内如果建筑横跨两个区within会将其剔除此时应改用intersects并按面积占比拆分但沈阳实际行政区边界与建筑轮廓误交叉的情况较少within的计算速度更快。include_groupsFalse是较新版本 GeoPandas 对apply行为的修正老版本会因分组列参与计算而报错。4.2 用层高字段推算总建筑面积建筑轮廓只有基底面积但很多业务指标需要总建筑面积。常见做法是用FLOORS层数字段乘以基底面积。# 清洗层数字段缺失或异常值处理 buildings[floors] pd.to_numeric(buildings[FLOORS], errorscoerce) buildings[floors] buildings[floors].fillna(1) buildings[floors] buildings[floors].clip(lower1, upper60) buildings[total_area] buildings.geometry.area * buildings[floors] print(buildings[total_area].sum())这里有个关键细节FLOORS字段如果是字符串类型to_numeric会把无法解析的值转成 NaN填充为 1 意味着默认按单层处理。对于明显偏离真实值的层数用clip限制上下限防止几栋异常高楼把全市统计结果拉偏。如果数据里只有HEIGHT而没有FLOORS可以按每层 3 米折算但沈阳老旧小区层高在 2.8 米到 3.1 米之间浮动折算误差无法避免。更可靠的做法是抽样比对随机挑 50 栋建筑用影像测量实际层数再和字段值对比确定该数据集的折算系数。4.3 与路网、POI 的叠加建筑到道路的距离与配套可达性建筑轮廓与路网的叠加能回答“这栋楼离最近道路多远”这类问题。计算最近道路距离的通用做法是对每栋建筑质心调用空间索引找最近的路段。from shapely.strtree import STRtree roads gpd.read_file(shenyang_roads.shp, encodingutf-8).to_crs(epsg32651) # 只保留主干道以上等级避免街坊路干扰计算结果 roads roads[roads[class].isin([trunk, primary, secondary])] tree STRtree(roads.geometry.values) distances [] for g in buildings.geometry.centroid: nearest tree.nearest(g) distances.append(g.distance(nearest)) buildings[dist_to_road] distances参数说明STRtree在构建时对所有道路几何建立空间索引nearest方法直接返回最近几何的索引复杂度接近 O(log n)几万栋建筑在几秒内完成计算。过滤道路等级这一步很重要如果把小区内部路和步行道混进去距离结果会系统性偏小尤其沈阳老城区胡同密集不筛选会影响“临街建筑”判定。到这里一份建筑轮廓数据已经经历了“读取确认 → 清洗修复 → 坐标统一 → 转点成面 → 空间统计 → 路网距离”的完整链路。最后一步是把这些结果变成能交给业务方的图件。5. 出图与双标注用 QGIS 把沈阳轮廓数据做成成果图5.1 批量出图按区分幅打印沈阳中心城区的建筑轮廓密度大一次导出全市范围的图几乎无法看清细节。常见做法是按行政区或图幅批量出图QGIS 的 Atlas 功能可以完成。在打印布局中新建布局添加地图项并设置为Atlas模式覆盖图层选择区分界数据。然后在布局里插入标签内容使用 Atlas 表达式引用区名和统计指标例如[% NAME %] 区建筑覆盖率 [% round(coverage_rate * 100, 1) %]%设置完成后导出 PDF 或 PNG 时会按每个区自动生成一页。批量导出 10 个区单次完成格式统一。5.2 GIS 怎么显示两个标注一个表达式解决单独的建筑数据往往既要显示编号又要显示楼层数。QGIS 的标注系统默认只允许一个标注样式但可以用表达式把两个属性拼接在一行concat(FID, \n, FLOORS, F)这样标注会显示成两行第一行是建筑编号第二行是楼层数加 F 后缀。如果希望两个标注分别用不同颜色或字号可以启用“标注工具箱”里的多标注支持或直接在表达式里用 HTML 标签控制颜色concat(font color#2b2b2b, FID, /font, \n, font color#d43f3a, FLOORS, F, /font)在随后的地图导出时勾选“允许 HTML 格式”才能让字号和颜色生效。相比做两个标注图层叠加这种方式在图例和导出时都更干净。5.3 验证把轮廓画到影像底图上检查出图前最后一道工序是位置验证。做法很简单加载天地图或 Esri World Imagery 底图把建筑轮廓图层设置为半透明红色缩放至 1:2000 以下逐片区目检。重点分辨两类错位一是整个图层整体偏移数米说明坐标系或基准面不匹配二是局部区域偏移说明原始数据本身存在分幅接边问题。对于后者不要指望统一处理能解决只能定位到具体区域后用 QGIS 的“几何编辑器”手工调整。最后检查单栋建筑轮廓是否与影像上的屋顶阴影边界重合。太阳方位造成的阴影会让影像里的建筑边缘比真实轮廓偏移不能拿阴影边缘当标准应该以墙根线为准。这一步做完这份 2022 年沈阳建筑轮廓数据才真正具备对外发布和支撑业务分析的资格。本文还有配套的精品资源点击获取