深圳38万栋建筑体块Shapefile数据解析：从加载到空间分析

简介《2020年深圳市建筑体块数据》是一套基于GIS的SHP格式建筑空间数据集覆盖深圳全市约38万栋建筑采用WGS84坐标适用于城市规划、建筑密度分析、房地产研究等场景读者需具备基础GIS软件操作能力。压缩包共6个文件包含.shp几何数据、.dbf属性表、.prj坐标参考、.shx索引及配套元数据合计约34.01MB结构精简便于直接加载。目前已有405人学习下载。数据可支撑空间分析、叠加分析、缓冲区与聚类分析等操作借助属性表中的建筑面积字段还能快速统计总量、排序或对比2020年城市建筑格局为研究深圳城市化进程和建筑趋势提供可靠基础数据。1. 从一份38万栋的数据包说起2020年深圳市建筑体块数据到底能做什么当我第一次从下载链接拿到2020年深圳市建筑体块数据.zip解压之后发现文件比印象里多Shenzhen.shp、Shenzhen.dbf、Shenzhen.shx、Shenzhen.prj、Shenzhen.CPG还有Shenzhen.shp.xml。这份数据记录的是2020年深圳38万栋建筑的体块轮廓用WGS84地理坐标系存储每个多边形对应一栋建筑的占地面积。对城市规划工程师来说价值是直接用多边形做建筑密度、容积率和街区尺度分析对做GIS开发的人而言它又是测试空间叠加、缓冲区分析和三维体块拉伸的基础样例。下面先从Shapefile的文件结构拆起再走一遍加载、查询和面积统计的关键步骤最后聊一聊环方向与zip包处理的收尾细节。2. Shapefile文件族六个文件各自的功能边界和读取逻辑Shapefile不是单一文件。Shenzhen.shp只存几何坐标属性全部放在Shenzhen.dbf里。dbf是dBASE时代留传下来的表格式文件字段名长度被限制在10个字节内读取时稍不注意就会出现字段名截断和类型不匹配。Shenzhen.shx是对应的几何索引记录每个要素在shp文件中的字节偏移缺少它时多数桌面软件仍能打开图层但会明显变慢。Shenzhen.prj保存的是坐标系描述的WKT文本QGIS和ArcGIS都靠它判断图层应该用哪种坐标参考系。Shenzhen.CPG通常只有一行字指定dbf的编码。Shenzhen.shp.xml属于元数据写清了数据说明和属性项含义不影响空间数据本身的加载。需要澄清一个常见的解释shapefile不是“一个文件一所建筑”Shenzhen.shp这一个文件就包含了38万条建筑要素。后面所有空间操作都在这个单一图层上完成不再涉及一对多的文件拆分这并不影响分析方式只是避免把图层与要素混为一谈。2.1 六个文件的职责与缺失后果先列一个总表后续所有操作都围绕这个文件族展开。文件作用缺失或被改名的后果Shenzhen.shp存储建筑多边形的坐标无法加载几何Shenzhen.shx记录几何记录的偏移量随机读取变慢Shenzhen.dbf存储建筑属性如高度、面积属性表为空Shenzhen.prj声明WGS84坐标参考被当作未知坐标系Shenzhen.CPG声明dbf的字符编码中文属性乱码Shenzhen.shp.xml元数据含坐标系和属性说明仅影响元数据浏览在ArcGIS中如果同目录下存在Shenzhen.sbn和Shenzhen.sbx那是ArcGIS自己生成的空间索引不是Shapefile规范的必要组成部分。用QGIS读取时这两个文件会被忽略因为QGIS默认建立自己的格网索引。这份数据里没有.sbn/.sbx不是问题但拿到ArcGIS里做叠加分析时软件可能会提示重新构建空间索引选择确认即可。2.2 WGS84与.prj为什么经纬度坐标系下不能直接算面积打开Shenzhen.prj可以看到WGS_1984的字样。WGS84是地心坐标系经纬度单位为度适合做全球定位和跨区域分析。深圳的公开数据里还有相当一部分采用CGCS2000或深圳地方坐标和WGS84在平面位置上的差异通常只有几米到几十厘米。做城市级建筑密度这类大范围统计时这种差异在结果里体现得不明显但要算单栋建筑面积、建筑间距就必须先把数据转到投影坐标。判断图层当前坐标范围最快的方式是用ogrinfo查看Extentogrinfo -ro -so Shenzhen.shp Shenzhen | grep -E Extent|Geometry|Feature Count输出里Feature Count对应约38万Geometry是PolygonExtent的X范围应在113.7到114.6附近Y范围在22.4到22.8附近这说明数据还是经纬度。如果X范围是3260000、Y是2500000那就是已经做好的投影坐标面积计算无需再次投影。Windows终端没有grep时直接去掉管道部分运行即可。这里用到-ro只读模式-so只输出摘要信息避免把38万条要素全部打出来。后续做文件转换时我最常用的投影参数是EPSG:32650也就是WGS 84 / UTM zone 50N深圳全域都在这个分带内计算面积的最大变形控制在千分之一以内。2.3 多边形体块环、内环与顶点顺序Shenzhen.shp中的每个建筑要素都是一个Polygon由外环和内环组成。普通住宅只有外环带天井或庭院的大型公建则多一个内环内环在几何上表示该范围内的面积不应计入建筑基底。读取几何信息时直接使用pyshp可以绕过桌面GIS快速验证数据结构import shapefile sf shapefile.Reader(Shenzhen.shp) print(sf.shapeTypeName) # 输出 POLYGON first sf.shape(0) print(len(first.points)) # 第一个建筑包含的顶点数 print(first.parts) # 环的起始索引shapeTypeName返回POLYGON说明这份数据是面要素。parts数组存储每个环在大坐标列表中的起点位置长度大于1就意味着该建筑有内环。环的顶点顺序在不同生产方那里并不统一Shapefile规范建议外环顺时针、内环逆时针但很多采集工具输出的方向完全相反单独看不影响渲染参与布尔运算和面积计算时则可能带来问题。这个坑在最后一章会专门处理。3. 加载与测绘把38万栋建筑读进GIS之前先处理Shenzhen.shp的编码与字段这一章直接进入实际操作。无论之后要用QGIS还是ArcGIS加载shapefile之前都要先确认三件事坐标参考是否是WGS84、属性表能不能正常显示中文、字段类型是否适合后面做统计。这三件事在前处理阶段解决比在分析做到一半再回头处理省事得多。3.1 在QGIS中添加Shapefile并验证坐标系QGIS打开新项目后通过“图层”菜单找到“添加图层”并选择“添加矢量图层”数据源类型选“文件”定位到解压目录里的Shenzhen.shp。添加完成后图层面板里的图层名称右侧会显示当前坐标参考正常情况下是WGS 84 / EPSG:4326。这个参考信息来自Shenzhen.prj如果软件显示“未知CRS”就不要继续做面积和距离相关的操作。加载后第一件事是快速查看多边形是否存在明显错误。双击图层打开“图层属性”在“符号化”里把填充改成半透明再用“信息”工具点击任意建筑弹出的属性框会同时显示几何类型和字段值。如果点击时提示“无效几何”说明原始要素存在自相交或环断裂需要先修复再分析。也可以用第二章提到的ogrinfo做加载前检查但命令行方式看不到图形的实际位置。更稳妥的组合是先用ogrinfo确认四至坐标再用QGIS看图形两个结果相互印证。3.2 读取dbf属性字段不要凭文件名猜字段双击打开属性表之前先读取一下Shenzhen.dbf的表结构避免后面写查询条件时把字段名写错。命令行方式ogrinfo -al -so Shenzhen.shp输出里会列出每一个字段的名字、类型和宽度。国内发布的shp中这段属性可能包含建筑高度、层数、竣工年份、建筑面积等但具体字段名和类型应当以这里的实际输出为准。很多分析脚本出错都是因为把字段名写成了英文缩写而实际字段是中文。用pyshp读取的等价操作是这样的import shapefile sf shapefile.Reader(Shenzhen.shp) for field in sf.fields[1:]: print(field[0], field[1], field[2], field[3])字段信息打印出来是四元组依次为字段名、字段类型、长度、小数位。常见类型中C代表字符串N代表整数F代表浮点数D代表日期。查询条件里如果对字符串字段写数字比较软件不会报错但结果一定是空。3.3 中文乱码CPG文件与编码修复2020年深圳市建筑体块数据中的属性如果包含中文Shenzhen.CPG就决定了dbf按什么编码读取。解压后的.CPG内容可能是UTF-8也可能是GBK。如果QGIS属性表里出现“淇℃伅”这类乱码基本可以断定CPG文件缺失或内容与实际编码不一致。修复方法非常直接在shp文件同一目录下重写CPG文件echo UTF-8 Shenzhen.CPG写完后重新加载Shenzhen.shp。如果乱码仍然存在把UTF-8改成GBK再试。注意CPG文件必须和shp主文件名完全一致扩展名是.CPG放在同一目录。这个文件虽然小但它比dbf里的编码标记更早被GIS软件读取改错方向会影响后续所有中文属性查询。4. 空间分析从建筑体块算建筑面积、密度与覆盖范围拿到Shenzhen.shp后最常见的需求是统计建筑面积、筛选大型建筑、计算片区建筑密度。这三类操作分别对应三个技术点面积量算的投影前提、属性过滤的语法、以及多边形叠加的空间关联。下面逐个说明。4.1 先投影再算面积UTM 50N是默认选择WGS84坐标系下的多边形面积用GIS直接算出来的单位是“平方度”平方度不是真正的面积单位不同维度下同一平方度代表的实际面积差别很大。要得到平方米第一步是把Shenzhen.shp转到投影坐标系。命令行转换ogr2ogr -t_srs EPSG:32650 -lco ENCODINGUTF-8 Shenzhen_UTM50N.shp Shenzhen.shp这条命令生成一个新的Shenzhen_UTM50N.shp坐标系是WGS 84 / UTM zone 50N。参数-t_srs EPSG:32650指定目标投影-lco ENCODINGUTF-8确保输出dbf继续保持UTF-8编码。转换之后在Python里统算所有建筑几何面积from osgeo import ogr ds ogr.Open(Shenzhen_UTM50N.shp) layer ds.GetLayer(0) total 0.0 feature layer.GetNextFeature() while feature: geom feature.GetGeometryRef() if geom is not None: total geom.GetArea() feature layer.GetNextFeature() print(total geometry area:, round(total, 2))这段代码遍历图层用GetArea()获取每个多边形的平面面积。输出单位是平方米因为输入Shenzhen_UTM50N.shp的坐标单位是米。需要注意的是这里算的是所有建筑多边形的投影面积总和和官方统计里的“总建筑面积”未必一致后者通常指各楼层的建筑面积之和体块数据只能表达楼层平面轮廓。4.2 属性过滤用SQL或表达式筛出大建筑如果只想保留面积超过10000平方米的大型建筑用ogr2ogr加where从句生成子集ogr2ogr -t_srs EPSG:32650 -where area 10000 Shenzhen_large.shp Shenzhen.shpwhere后面是标准SQL表达式字段名用双引号数字不加引号。字段名如果带中文在部分终端里需要额外处理字符集建议先把dbf字段复制为英文字段再过滤。在QGIS字段计算器里可以用$area 10000直接筛选不依赖具体字段名需求表达式说明几何面积大于1万平米$area 10000当前坐标系下计算缩放多边形做分析buffer($geometry, 50)以米为单位的缓冲区中心点坐标x($geometry)/y($geometry)取多边形质心在QGIS中使用这些表达式时如果没有先投影$area会沿用图层的坐标单位结果会变成难以理解的“平方度”。城市级分析一般先另存为UTM投影文件再做计算这样表达式结果才能输出为平方米。4.3 叠加分析把建筑密度落到社区或地块如果想计算某个行政单元内建筑占地面积占单元面积的比例需要把建筑体块同行政边界做空间连接。最直接的方式是使用QGIS Processing算法在QGIS Python控制台中执行processing.run(qgis:joinattributesbylocation, { INPUT: Shenzhen_UTM50N.shp, JOIN: community_boundary.shp, PREDICATE: [0], OUTPUT: joined.shp })PREDICATE参数中的0表示intersects即建筑多边形与目标边界有重叠就分组。这个操作会把行政单元的唯一标识写到每个建筑要素上之后用“按属性统计”工具对同一标识做面积求和再除以单元总面积就是建筑密度。对于更大规模的数据我会把shp导入PostGIS用空间SQL一次完成分组汇总SELECT c.name, sum(ST_Area(b.geom)) AS building_area, sum(ST_Area(c.geom)) AS unit_area FROM buildings b JOIN communities c ON ST_Intersects(b.geom, c.geom) GROUP BY c.name;这段SQL用ST_Intersects把建筑数据映射到社区ST_Area在投影坐标下得到平方米最后分组输出每个单元的建筑占地和单元面积。使用前要先建空间索引否则几十万要素的关联在常规机器上会跑很久CREATE INDEX idx_buildings_geom ON buildings USING gist(geom);CREATE INDEX ... USING gist(geom)为建筑几何字段建立GiST空间索引PostGIS执行空间连接时会优先走索引大幅缩短查询时间。若导入时未做投影两表坐标不一致会导致ST_Area结果错误这一点比索引问题更隐蔽。5. 环方向修复与zip包版本管理交付Shenzhen建筑体块数据的两个收尾技巧5.1 环方向与无效几何的检查修复基础面积计算不会暴露环方向问题但做布尔运算、三维体块拉伸或拓扑检查时环方向不一致会让结果异常。用shapely统一环方向是轻量做法import geopandas as gpd from shapely.ops import orient gdf gpd.read_file(Shenzhen_UTM50N.shp) print(gdf.is_valid.sum(), gdf.shape[0]) # 有效要素数 / 总要素数 gdf[geometry] gdf.geometry.map(lambda g: orient(g, sign1.0)) gdf.to_file(Shenzhen_oriented.shp, encodingutf-8)orient(g, sign1.0)强制外环按逆时针排列内环自动反向sign-1.0则相反。多数场景推荐sign1.0它在很多几何算法里被当作标准方向。is_valid.sum()统计有效要素数量用于评估原始数据是否需要进一步修复。如果is_valid结果出现大量False先做零缓冲修复gdf[geometry] gdf.geometry.buffer(0)buffer(0)不改变图形主体但会把自相交和重复顶点一并消除。经过这步修复后重新检查is_valid直到有效要素数量接近总数再进入分析流程。5.2 zip包内文件同名与ArcGIS互操作最后收尾是文件交付。Shenzhen.shp只是一整套文件的领衔角色改名时如果只改主文件名其他五个同伴不跟进加载必然报错。比如把Shenzhen.shp改成sz_buildings_2020.shpShenzhen.dbf、Shenzhen.shx、Shenzhen.prj、Shenzhen.CPG也必须同步改名前缀保持完全一致。压缩成zip时建议把六个文件放在压缩包根目录不要在包内再建一层文件夹。ArcGIS在没有安装Data Interoperability扩展时不会把zip包当作可浏览数据源QGIS虽然支持直接读取zip内的shapefile但对中文路径和中文文件名的兼容性仍然不佳所以我一般先解压再把文件名改成拼音或英文最后重新打包。如果下载的zip报错“Error reading zip archive”优先考虑重新下载并比对文件大小而不是急着用第三方修复工具。zip是顺序读取结构下载中断最容易损伤中央目录区修复工具只能挽救部分尾部数据对shp这种二进制定长记录格式来说修复成功率不高。建议在原始副本上先做环方向修复再执行投影和zip打包这样每个交付版本都可以独立重复校验。本文还有配套的精品资源点击获取