★ FEATURED ARTICLE 简介北京市行政区划矢量数据包面向GIS初学者、城市规划与人口统计等空间分析人员提供北京市及下辖区县行政边界的shp格式数据可直接用于制图、空间查询与专题分析。压缩包共11个文件约82KB其中shp存储几何边界dbf保存属性字段shx提供空间索引prj定义坐标投影xml记录元数据lock为软件运行时的临时锁文件配套结构完整便于在ArcGIS、QGIS等工具中直接加载。已有3677人学习下载说明其在行政区划底图需求中具有较高参考价值。读者可获得北京市行政区域的多边形边界与属性表用于叠加人口密度、交通流量等图层开展区县对比、规划决策支持和基础地图制作也可作为教学练习数据帮助理解shp多文件协同机制与坐标投影信息。数据体量轻巧适合快速预览与二次整理。1. 北京市 shp 文件到手先看三样dbf 编码、prj 坐标系、shx 索引上周有人拿着标题写着「北京市 shp 文件」的压缩包来问QGIS 里属性表全是问号拖到在线底图上整体偏出去几十公里按字段算出来的区县面积还是负数。三个症状对应三件完全不同的事——dbf 里的字符编码、prj 里声明的坐标系、多边形环的绕向。市级行政区划的 shp 格式数据几乎是所有空间分析的起点选址要拿它切商圈统计要拿它把点位归到区县可视化要拿它当轮廓底图。这份资源适合三类人需要区县、街道两级边界能直接用的数据分析同学要把它灌进 PostGIS 或者切成瓦片的后端以及被编码和坐标系反复折腾过的 GIS 新手。下面按「先体检、再动刀、最后验收」的顺序走命令和参数都给到能直接抄的程度。2. shp 格式不是单个文件七个附属件各自管什么很多人对 shp 格式的第一印象是「一个 .shp 文件」实际上它是一套同名文件的组合。压缩包里如果只挑了 .shp 拷出来几何确实还在但属性、索引、坐标系、编码全丢了后面每一步都会出问题。这一章先把文件清单和各自的职责讲清楚再给一条能把家底一次性翻出来的命令。2.1 .shp 只装几何属性和索引都在隔壁扩展名职责缺失后的表现.shp存几何坐标点是点、面是顶点序列文件直接打不开.shx几何的定长索引按 fid 快速定位能打开但随机读取极慢部分库直接报错.dbfdBase 格式的属性表只剩几何没有任何名称和代码字段.prjWKT 形式的坐标系声明crs 为空面积、距离、叠加全部失去依据.cpg一个纯文本文件写UTF-8或CP936中文变问号或者被按错误的代码页解码.sbn / .sbxEsri 自己的空间索引只影响查询性能不影响正确性.qixQGIS 生成的空间索引同上重新生成即可.shp.xml元数据只影响溯源不影响使用这套格式的设计年代很早限制也一并继承下来了单个文件不能超过 2GB字段名最长 10 个字符中文名基本等于 5 个字一个图层只能有一种几何类型点线面不能混存不支持真正的曲线圆弧都要先离散成折线面要素靠环的绕向来区分外环和洞外环顺时针、内环逆时针绕反了面积就是负的日期字段的格式在 dbf 里也不统一。这些限制不是陈年旧账只要你打算把数据交给别人它们就会以「字段名被截断」「面积算反」的形式冒出来。2.2 一条 ogrinfo 把家底翻出来拿到压缩包先解压到一个纯英文路径的目录然后别急着开 QGIS用 GDAL 自带的ogrinfo看一眼。-so表示 summary only不遍历要素几秒钟出结果-al表示所有图层。两个参数一起用就是一份完整体检报告。# 概览要素数、几何类型、范围、字段、坐标系 ogrinfo -al -so ./beijing_admin.shp # 单独把 .prj 翻译成 EPSG 代号比肉眼读 WKT 快 gdalsrsinfo -o epsg ./beijing_admin.prj # 抽一条要素看属性值确认中文没坏 ogrinfo -al -fid 1 ./beijing_admin.shp正常输出大致长这样重点看四行Layer name: beijing_admin Geometry: Polygon Feature Count: 16 Extent: (115.423411, 39.442571) - (117.514343, 41.060799) Layer SRS WKT: GEOGCRS[China Geodetic Coordinate System 2000, ...] adcode: String (6.0) name: String (30.0)Geometry说明这一层是面Feature Count是要素数区县一级通常是十几个街道一级会到几百个两个数量级对不上就说明拿错了层Extent的数值落在 115118 和 3941 之间说明坐标存的是经纬度单位是度Layer SRS WKT里出现China Geodetic Coordinate System 2000就说明 prj 在而且声明的是 2000 国家大地坐标系。如果这一行是空的直接跳到 2.3 补 prj不要往下做任何计算。字段那一行也值得多看两眼。String (6.0)里的 6 是宽度行政区划代码正好 6 位但 dbf 的宽度是写入时定死的超过就会被截断。中文名如果被压成district_n这种说明源数据写入时就已经截过了后面靠字段名做映射的代码最好改成按位置索引或者建一张映射表。2.3 缺 .prj 和缺 .cpg 的两种补救补 prj 有一个前提必须先说清楚赋值和转换是两件事。-a_srs是给数据贴标签坐标数值一个字节都不动适用于「我知道它是什么坐标系只是声明丢了」-t_srs是真的做投影换算适用于「它现在是 A 坐标系我要 B 坐标系」。贴错标签比没标签更麻烦因为后续所有计算都会在错误的假设下静默跑完。# 1) 只知道是 CGCS2000 经纬度、prj 丢了赋值不动坐标 ogr2ogr -a_srs EPSG:4490 fixed.shp raw_no_prj.shp # 2) 直接把 prj 文本写出来省一次全量拷贝 gdalsrsinfo -o wkt EPSG:4490 raw_no_prj.prj # 3) 缺 .cpg新建一个纯文本文件内容就一行 printf UTF-8 raw_no_prj.cpg # 4) 源数据是 CP936想统一重写成 UTF-8 的 dbf ogr2ogr -f ESRI Shapefile -lco ENCODINGUTF-8 utf8_out.shp gbk_in.shp第 1 条里的-a_srs只接受坐标系定义不校验数据本身写错了不会报错所以赋值前最好拿一个已知点交叉验证——比如取区政府驻地的经纬度和属性表里的坐标对一下量级。第 3 条的.cpg文件内容必须是标准写法UTF-8、CP936、GBK都能被 GDAL 和 QGIS 识别但写成utf8小写有时会被忽略。第 4 条的-lco ENCODING是图层创建选项只对输出生效源数据的编码它不会自动猜源端还是得靠.cpg或者显式参数说明。3. 4326、4490、4548 别混着用面积和长度算错多半栽在这坐标系是这份数据里最容易出事、又最不容易被发现的一环。因为投影错了不会报错它只会让结果偏一点点而你往往要等到出图或者汇报时才看出来。这一章把三个最常出现的代号摆在一起对比再给重投影和面积计算的可执行写法。3.1 地理坐标存位置投影坐标算距离EPSG:4326 是 WGS84 的经纬度单位是度EPSG:4490 是 2000 国家大地坐标系的经纬度单位也是度。两者在椭球定义上差异很小实际业务里经常被混用但混用的代价是米级的系统偏差做小比例尺展示无所谓做地块核算就要命。EPSG:4548 是 2000 坐标系下的三度带高斯克吕格投影中央经线 117°E不带带号单位是米。这三个的关系是前两个告诉你「这个点在地球上的位置」后一个告诉你「这个点在平面上离中央经线多少米」只有后者能直接算面积和长度。选带号还是不带号看东坐标的第一位。带号版本三度带第 39 带EPSG:4527的东坐标会是 39 开头的八位数形如 39500000 左右不带号版本EPSG:4548是六位数500000 上下。这两个混了图会瞬间跑到几千公里外但如果你只做相对计算可能一直发现不了。拿不准就用projinfo EPSG:4548查一下 central meridian 和 false easting 两个参数比翻文档快。提示城市级数据用三度带足够。经度跨度超过 3 度、或者要做全市横向对比时改用等积投影Albers更稳下面 3.3 给了参数。3.2 重投影怎么写ogr2ogr、GeoPandas、PostGIS 三选一命令行适合批量处理Python 适合和后续分析串起来数据库适合数据量上去了以后反复查询。三种写法给的是同一个结果挑顺手的用。# 命令行一次转换输出到新文件 ogr2ogr -t_srs EPSG:4548 -f ESRI Shapefile bj_4548.shp beijing_admin.shp-t_srs触发真正的坐标重算输出的 .prj 会自动写成目标坐标系。注意输出仍然是 shapefile字段名会被再截一次如果原字段名已经接近 10 个字符最好在这一步换成 GeoPackage 输出。import geopandas as gpd from shapely.validation import make_valid gdf gpd.read_file(beijing_admin.shp, encodingutf-8) print(gdf.crs) # 为空说明 prj 丢了先按 2.3 补 # 缺 crs 时显式赋值allow_override 允许覆盖已有的错误声明 if gdf.crs is None: gdf gdf.set_crs(EPSG:4490, allow_overrideTrue) # 先修几何再投影顺序反了自相交的环会被算得更歪 gdf[geometry] gdf.geometry.apply(make_valid) gdf_m gdf.to_crs(EPSG:4548) # 单位从度变成米read_file的encoding参数只影响 dbf 解码不影响几何。to_crs会调用 pyproj 做逐点换算几万个顶点是秒级百万级顶点就会明显卡这时候换ogr2ogr或者 PostGIS 更好。-- PostGIS入库时就带上 SRID后面查询不用反复转 SELECT name, ST_Area(ST_Transform(geom, 4548)) / 1e6 AS area_km2 FROM bj_admin ORDER BY area_km2 DESC;ST_Transform是逐要素调用配合ST_Area得到的是平方米除以 1e6 换成平方公里。如果这张表要反复算面积建议加一个投影列并建空间索引比每次查询都转要快一个数量级。3.3 面积、周长、等积投影的参数在经纬度坐标下调用.area返回的单位是「度的平方」数值大概在 0.00x 这个量级既不是平方米也不是平方公里拿来排序勉强能用拿来汇报一定翻车。正确做法是转到投影坐标再算。gdf_m gdf.to_crs(EPSG:4548) gdf_m[area_km2] gdf_m.geometry.area / 1e6 # 平方米转平方公里 gdf_m[peri_km] gdf_m.geometry.length / 1e3 # 米转公里如果要做跨城市的横向对比或者数据本身跨了几个三度带用等积投影更公平。Albers 的两条标准纬线决定变形分布中央经线决定形状的对称轴CHINA_ALBERS ( projaea lat_125 lat_247 lat_00 lon_0105 x_00 y_00 ellpsGRS80 unitsm no_defs ) gdf_ea gdf.to_crs(CHINA_ALBERS) gdf_ea[area_km2] gdf_ea.geometry.area / 1e6lat_125和lat_247是常用的一组标准纬线落在这两条纬线上的长度比精确为 1中间区域变形最小lon_0105让中央经线穿过国土中部东西两端的角度变形比较均衡ellpsGRS80和 2000 坐标系的椭球一致不要图省事写成WGS84。这套参数在城市级尺度上与高斯投影算出的面积差异通常在万分之一以内做统计口径对比完全够用。4. 属性表清洗与空间操作从街道边界到区县汇总几何对了接下来是属性和空间关系。这一章的三段内容是一条完整的流水线先把字段读对再做合并裁剪最后挑一个合适的落盘格式。中间任何一个环节偷懒后面都要返工。4.1 DBF 的 10 字符字段名和前导零中文乱码最省事的排查顺序是先看.cpg里写了什么没有就猜编码猜不中再上gb18030兜底。写个能自动试的函数比每次手动改参数快得多。import geopandas as gpd def read_shp_any(path): 按 utf-8 → gbk → gb18030 依次尝试返回 (gdf, 命中的编码) for enc in (utf-8, gbk, gb18030): try: return gpd.read_file(path, encodingenc), enc except UnicodeDecodeError: continue # gb18030 是超集兜底基本不会失败 raise RuntimeError(三种编码都解不开检查 .cpg 内容是否损坏) towns, enc read_shp_any(beijing_town.shp) print(命中的编码, enc)行政区划代码是另一个高频坑。dbf 里它可能被存成数值型pandas 读进来是int64或者float64一旦有前导零就没了街道一级是 9 位村级是 12 位超出 int32 范围后还可能变成科学计数法。统一按字符串处理顺手把浮点尾巴切掉towns[adcode] ( towns[adcode].astype(str) .str.replace(r\.0$, , regexTrue) # 去掉 float 读进来的 .0 .str.zfill(9) # 街道级补足 9 位 )zfill的位数按你的数据层级定区县级 6 位街道级 9 位村级 12 位。这一步做完后面按代码前缀做层级汇总才不会串。字段名被截断的问题没法在读取端解决只能建一张映射表把district_n之类的名字映射回业务含义或者干脆在这一步就用 GeoPackage 重新输出一次。4.2 dissolve、clip、sjoin 三段式街道合并成区县、按范围裁剪、把点位归到区县是这份数据最常见的三种用法。三个函数的分工很清楚dissolve按字段合并几何clip用另一个几何做裁切sjoin只判断空间关系、不改几何。# 街道合并成区县按区县代码分组几何做并集 districts towns.dissolve( bydistrict_code, as_indexFalse, aggfunc{district_name: first} # 非几何列怎么取必须显式指定 ) # 用一个矩形裁出感兴趣的范围 from shapely.geometry import box clip_box gpd.GeoDataFrame( geometry[box(116.20, 39.70, 116.60, 40.10)], crsEPSG:4490 # 必须和待裁剪数据同坐标系 ) sub gpd.clip(towns, clip_box) # 点位归到区县 points gpd.GeoDataFrame( points_df, geometrygpd.points_from_xy(points_df[lon], points_df[lat]), crsEPSG:4490 ) hit gpd.sjoin(points, districts[[district_name, geometry]], howleft, predicatewithin) print(没有落进任何区县的点, hit[district_name].isna().sum())dissolve的aggfunc不写会丢掉所有非几何列写first表示每组取第一个值前提是这些列在组内本来就一致。clip前不做to_crs对齐是最常见的静默失败坐标系不一致时 GeoPandas 会直接抛异常反倒是种保护。sjoin的predicate默认是intersects边界上的点会被重复计入相邻两个区县用within可以避免重复但正好压在边界线上的点会被丢掉。两种口径都合理关键是要在结果里说清楚用的哪一种别让下游以为总数对得上。4.3 落地格式选 GeoPackage 还是 GeoJSON交付给人看、给系统读、给前端用三种场景对应三种格式。shapefile 适合当原始档案留存不适合当交付件。# GeoPackage单文件、UTF-8 原生、字段名不限 10 字符适合当交付件 ogr2ogr -f GPKG beijing.gpkg beijing_admin.shp -nln admin -lco SPATIAL_INDEXYES # GeoJSON按 RFC 7946 规定必须是 WGS84 经纬度所以要显式转一次 ogr2ogr -f GeoJSON -t_srs EPSG:4326 beijing.geojson beijing_admin.shp # 多个图层打进同一个 gpkg靠 -nln 区分 ogr2ogr -f GPKG beijing.gpkg beijing_town.shp -nln town -update-lco SPATIAL_INDEXYES会顺手建空间索引后面按范围查询不用全表扫。-nln指定图层名不加就会沿用源文件名。-update表示往已有文件里追加图层第一次创建时不用加。GeoJSON 那条的-t_srs EPSG:4326不是可选项很多前端地图库读到非 4326 的 GeoJSON 会直接把坐标当经纬度用图会跑到南极洲附近。5. 避坑与排查五条亲眼见过的翻车记录下面这五条不是我编的是这几年在别人机器上真实复现过的。每一条都按「现象 → 原因 → 解决」写遇到类似的先对号入座比从头排查省时间。5.1 属性表里中文全是问号或者方块现象QGIS 打开图层属性表里所有中文都是???或者一排方块但用 Excel 打开同目录的 .dbf 却是正常的。原因dbf 本身不带编码信息全靠旁边的.cpg声明。压缩包里.cpg丢了或者内容写的是936这种非标准写法QGIS 会按 Latin-1 兜底解码中文自然全废。用 Excel 打开正常是因为 Excel 在中文环境下默认按 CP936 猜。解决先看.cpg在不在不在就新建一个内容写CP936或者UTF-8重启 QGIS还是不行就在图层属性的「源」里手动指定编码然后「另存为」一份 UTF-8 的新文件一劳永逸。5.2 叠到在线底图上整体偏移几百米现象shp 里的区县边界和在线底图的路网对不上整体平移了几百米而且偏移方向在各个区之间还不完全一致。原因两边的坐标系不是同一套。一种情况是 prj 声明错了比如数据实际是 2000 坐标系却被写成了 WGS84另一种情况是底图服务用的坐标系和 shp 里的不是一套叠加时没有做转换。精度要求高的场景下这种偏移足以让一条街落进隔壁区。解决先用 3.1 的方法确认 prj 声明的到底是哪个代号再用gdalsrsinfo输出 WKT 和已知点交叉验证确认无误后在 QGIS 里统一工程坐标系然后对底图和矢量分别检查「图层 CRS」是否被正确识别。别用肉眼拖拽去对齐那只是把误差藏起来了。5.3 面积是负数或者比预期大出两个数量级现象算出来的区县面积有正有负或者全都是零点零零几。原因负数来自面的环绕向反了外环逆时针、内环顺时针鞋带公式算出来就是负值零点零零几来自直接在经纬度坐标下调用.area单位是「度的平方」。解决先to_crs到投影坐标系再算负值自然消失如果转完还是负的用make_valid修一遍几何再算。批量处理时加一句断言assert (areas 0).all()让问题在流水线里就暴露出来别等到出图才发现。5.4 buffer 和 clip 报拓扑异常现象对全市面数据做缓冲区或者叠加跑到某一个区就抛TopologyException提示某某坐标处 self-intersection。原因源数据是拼接来的相邻区县的边界不严丝合缝或者单个面里存在自相交的环。shapefile 不存拓扑关系这种错误在写入时不会被拦截只在做几何运算时爆出来。解决运算前统一走一遍gdf.geometry gdf.geometry.apply(make_valid)仍失败就先把数据转成 GeoPackage 或者 GeoJSON 再读回来格式转换过程会顺手修掉一批无效环。数据量大时用 PostGIS 的ST_MakeValid配合ST_Buffer(geom, 0)后者能把自相交的环「熨平」。5.5 读到一半报错或者字段被悄悄截断现象read_file报解析错误或者读进来发现某个字段的所有值都只剩前几个字。原因前一种是 dbf 和 shp 的记录数不一致通常是拷贝时只拷了 .shp 没拷 .dbf或者两者被不同版本的软件先后改过后一种是 dbf 的 10 字符字段名限制源数据写入时就截断了读取端无能为力。解决先用 2.2 的ogrinfo -al -so核对Feature Count再和属性表行数对一下不一致就找回原始压缩包重新解压别在残缺数据上做修补。字段截断只能通过建立映射表解决或者让上游改用 GeoPackage 交付。6. 交付前跑一遍拓扑校验简化与切片参数怎么定数据能读、能算不等于能交付。区县边界这种数据最容易在「相邻两块之间多出一条缝」或者「两块互相压盖」上出问题肉眼看不出来一放到大比例尺下就露馅。6.1 有效性、压盖、缝隙三道检查有效性用is_valid压盖用overlay自相交缝隙靠合并前后面积差反推。import geopandas as gpd from shapely.validation import make_valid gdf gpd.read_file(beijing.gpkg, layeradmin).to_crs(EPSG:4548) print(无效几何, (~gdf.geometry.is_valid).sum()) gdf[geometry] gdf.geometry.apply(make_valid) # 压盖检查自己和自己求交去掉同一要素的自身相交 ov gpd.overlay(gdf, gdf, howintersection, keep_geom_typeTrue) ov ov[ov[adcode_1] ! ov[adcode_2]] print(互相压盖面积(km2), round(ov.geometry.area.sum() / 1e6, 4)) # 缝隙检查合并后的总面积应小于等于各块面积之和 total gdf.geometry.area.sum() / 1e6 merged gdf.dissolve().geometry.area.sum() / 1e6 print(合并差值(km2), round(total - merged, 4))压盖面积在 0.01 平方公里以下通常是浮点误差可以放行超过 1 就要人工看图。合并差值同理差值明显为正说明块与块之间有空隙差值明显为负说明有压盖。这三个数我一般会写进交付说明里比口头保证可靠。6.2 简化容差与瓦片参数做前端展示时几万个顶点的边界文件直接扔给浏览器很浪费。简化要在投影坐标下做容差用米preserve_topologyTrue可以避免简化后出现自相交。# 在米制坐标下简化容差 5 米肉眼在大比例尺下看不出差别 gdf[geometry] gdf.geometry.simplify(5, preserve_topologyTrue) gdf.to_file(beijing_lite.gpkg, layeradmin, driverGPKG)切成瓦片时简化交给 GDAL 做单位是像素SIMPLIFICATION1表示容差约 1 个像素ogr2ogr -f MVT beijing.mbtiles beijing.gpkg admin \ -dsco MINZOOM6 -dsco MAXZOOM13 \ -dsco EXTENT4096 \ -dsco SIMPLIFICATION1MINZOOM决定最小的那一级市级数据从 6 起步够用MAXZOOM13对应街区级再往下切文件会膨胀得很快EXTENT4096是瓦片内的坐标分辨率比默认的 4096 更细会给前端增加负担一般不用改。切完记得抽查两三个缩放级别重点看边界在低层级是否断裂——那是简化容差设太大的典型症状。这份数据我前后用过好几轮最深刻的教训是任何一次交付前必须把「有效性、压盖、缝隙」三个数跑出来并且和上一版对比。有一次省了这一步一个街道的边界少了两个顶点导致相邻两个区的面积一个多了一个少了下游按面积做的配额算完才发现对不上回滚重做花了整整两天。从那以后我每次导出前都强制走一遍这套校验宁可多花五分钟。希望帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
