简介:这份合肥功能区区划 shp 数据面向城市规划、地理研究与GIS制图学习者,提供超出常规行政区县的功能区空间边界,涵盖湖滨新区、经济开发区、经济开发区北区、高新区和新站区等单元,可用于空间结构分析、专题制图与规划政策研究。压缩包共8个文件、约276KB,以 shp 矢量主体配合 dbf 属性表、prj 坐标定义及索引与元数据文件组成,便于在 ArcGIS、QGIS 等工具中直接加载并保留区域名称、面积等属性信息。已有298人学习,说明该数据在区域研究与制图练习中具有一定参考价值。读者可获得较完整的功能区划分底图,用于叠加交通、生态或产业数据,观察各功能片区的分布关系,也能借此熟悉 shapefile 多文件协同、坐标投影与属性关联的处理流程,适合作为课程作业、论文前期分析或城市专题地图的基础数据。
1. 功能区划 shp 数据:它替掉的是两周矢量化工时
做空间统计,最贵的从来不是写代码,是找底图。上个月帮某实验室做商业网点覆盖分析,需求单上只有一句话:按功能区口径出一版统计。开工才发现,网格数据一抓一大把,能跟业务表对上号的区划面却要自己描——描完还要校准、修拓扑、跟业务方确认口径,两周就搭进去了。这份合肥功能区区划 shp 数据解决的正是这一步:把功能区这一层的面要素做成标准 shapefile,GIS 软件能直接拖进去,GeoPandas、PostGIS 能直接读,省掉的是从零矢量化那段最枯燥的工期。
它适合做城市研究、选址评估、网点覆盖、公共服务设施均衡性分析的人。手上有一张带经纬度的点位表,或者一份按片区归集的业务明细,需要落到“面”上做汇总,这层区划就是底座。反过来,要找实时更新的行政边界、或者精确到宗地级别的,得另找来源,区划面是宏观口径,别拿它做微观确权。下面按认清结构、跑通加载、改成业务可用、排查现场问题、进阶验证的顺序走一遍。
2. 拆开一份区划 shp:图层组成、字段口径与坐标系
2.1 一个 shp 其实是一组文件
新手最容易翻车的地方在这儿:从压缩包里只把.shp拖出来用,打开一看图形在、属性表空了,或者图形在、位置飘到海上去。Shapefile 从来不是单文件格式,它是一个文件族,几何、属性、投影声明各自躺在不同文件里,少一个就残。拿到包先解压到整目录再动,拷贝的时候整个文件夹一起走。
| 扩展名 | 装的是什么 | 缺失后的典型症状 |
|---|---|---|
.shp | 几何体(点/线/面) | 什么都打不开 |
.shx | 几何索引 | 部分软件能开但极慢,或报格式错误 |
.dbf | 属性表 | 图形正常,属性表全空 |
.prj | 坐标系 WKT 声明 | 图形偏到 0,0 附近或被当成经纬度错位 |
.cpg | 字符编码声明 | 中文字段变问号或乱码方块 |
.sbn/.sbx | ESRI 空间索引 | 不影响打开,只是首次渲染变慢 |
.qix | MapInfo 空间索引 | 同上,缺了自动重建 |
.shp.xml | 元数据 | 丢的是来源说明,不影响几何 |
.cpg和.prj这两个是重灾区。.cpg里通常就一行UTF-8或GBK,很多导出工具根本不生成它,导致 QGIS 按系统默认编码猜,字段一打开就是乱码。.prj更隐蔽,缺了之后软件往往默认按 WGS84 处理,图层跟底图叠上去差几百米,你肉眼还未必看得出来。
我一般习惯先把整个包解压到纯英文的短路径,比如D:/gis/hf_zone/,而不是丢在桌面上带中文和空格的目录里。老版本 GDAL 和部分桌面软件对中文路径的处理一直不干净,报错信息还特别含糊,“文件不存在”或者“无法识别的格式”,查半天才发现是路径问题。这条纯属玄学级别的血泪经验,但确实省事。
2.2 属性表先看清哪几类字段
几何是骨架,属性才是能干活的部分。区划类数据打开后,属性表里的字段大致分四类,看懂了才知道后面怎么用:
- 标识类:区划名称、区划代码。这是跟业务表 join 的主键,最要紧。
- 分类类:功能区类型、等级。做分级配色、分组统计都靠它。
- 量算类:面积、周长。注意这是导出时算好的静态值,口径可能是投影面积也可能是椭球面积,跟后面自己算的常常对不上。
- 辅助类:备注、来源、更新日期。看着没用,追溯数据版本的时候救命。
.dbf有三个天然限制,绕不过去。字段名最长 10 个字符,超了就被截断,中文名更惨,直接变字段1;文本字段最大长度 254 字符,长描述会被砍;数值精度有限,存经纬度这种要小数点后七八位的,容易丢精度。
最坑的是编码字段的类型。区划代码这类东西,一旦导出方把字段定义成数值型,010203读进来就变成10203,前导零没了。你拿它去 join 业务表,一条都对不上,还不报错,就是匹配率为零。读的时候必须显式指定类型:
import geopandas as gpd import fiona # 先只看字段结构,不加载几何,速度快 with fiona.open("D:/gis/hf_zone/zones.shp", encoding="GBK") as src: print(src.schema) # 字段名 -> 类型,重点看代码字段是不是 int print(src.crs) # 打印出来是 None 就说明 .prj 丢了 print(len(src)) # 要素条数 # 带类型读,把标识类字段强制成字符串,避免前导零丢失 gdf = gpd.read_file( "D:/gis/hf_zone/zones.shp", encoding="GBK", # 跟 .cpg 声明保持一致 dtype={"qhdm": str, "code": str} # 换成你实际的字段名 ) print(gdf.dtypes)fiona.open那句只读元数据不开几何,几十万要素也是秒回,适合先摸结构。dtype参数是保命的,字典的键要换成属性表里真实的字段名——不同来源的区划数据字段名五花八门,QHDM、XZQDM、CODE都见过,先print(src.schema)看清楚再写。
2.3 坐标系与范围:投影没定,后面全白干
坐标系这事,先分两类。地理坐标系用经纬度,单位是度,EPSG:4326(WGS84)和EPSG:4490(CGCS2000)在国内数据里最常见,两者差异在厘米到米级,做宏观统计基本无所谓,做高精度量算就得问清来源。投影坐标系直接用米,算面积、算距离、建缓冲区都必须先把坐标系切到投影。
不确认坐标系就直接算面积的,结果能差出三成以上。举个具体的:Web Mercator(EPSG:3857)在纬度 31.8° 附近,面积变形因子约1/cos²(φ) ≈ 1.39,也就是说同一块地在这个投影下量出来会虚胖将近四成。这种偏差混在一堆报表里,你看数字是合理的,直到跟别的口径一对账才发现全歪了。
先用命令行把底摸清,ogrinfo是 GDAL 自带的,装了 QGIS 就有:
# -so 只输出摘要,-al 应用到所有图层,不打印几何坐标 ogrinfo -so -al D:/gis/hf_zone/zones.shp # 输出里重点关注这几行: # Geometry: Polygon # Feature Count: 23 # Extent: (117.1xxxxx, 31.6xxxxx) - (117.5xxxxx, 32.0xxxxx) # Layer SRS WKT: ... 没有这一行 = .prj 缺失Extent那四个数是判断坐标系对不对最快的办法。数值在 117、31 这个量级,说明是经纬度;数值在几十万、几百万,说明已经是投影坐标;如果看到(0, 0) - (1, 1)或者-180这种,基本可以判定坐标系声明出了问题,先别往下走。
.prj确实丢了的情况,别用to_crs(),那是给已有坐标系的图层做转换用的。正确做法是先set_crs(epsg=4326, allow_override=True)把坐标系“贴”上去,再to_crs()转到目标投影。顺序反了,数据会被转两次,偏到怀疑人生。
3. 跑通最小验证:QGIS 目视加 GeoPandas 体检的六步
3.1 QGIS 里先看图形,再看属性
拿到数据别急着写代码,先用 QGIS 目视一遍,五分钟能发现一堆脚本跑一天也发现不了的问题。步骤固定下来:
- 解压到英文路径,确认目录里
.shp/.shx/.dbf/.prj四件套齐全。 - 把
.shp直接拖进画布,或者图层 → 添加矢量图层 → 选文件。 - 看右下角状态栏的坐标系标识,跟
ogrinfo的输出是否一致。 - 右键图层 → 属性 → 信息,核对要素数、几何类型、范围。
- 打开属性表,切到字段视图看类型,重点盯编码字段和面积字段。
- 叠一层在线底图核对位置关系——这一步放到最后,别一上来就加载底图,底图的位置误差会干扰你对数据本身的判断。
渲染上有个小技巧:区划类数据按分类字段做分级配色,比统一填充色有用得多,一眼就能看出哪几个片区类型重复、哪个片区属性缺失。图层样式里选“分类”,字段选功能区类型,点“分类”按钮自动配色。
3.2 GeoPandas 做一次几何体检
QGIS 能看出“长得对不对”,看不出“几何合不合规”。自相交、空几何、多部件混入,这些在渲染时看着都正常,一算面积、一做叠加就炸。体检脚本我一般固定跑这几项:
import geopandas as gpd import shapely gdf = gpd.read_file("D:/gis/hf_zone/zones.shp", encoding="GBK") print("要素数:", len(gdf)) print("几何类型:", gdf.geom_type.value_counts().to_dict()) print("坐标系:", gdf.crs) # 空几何 / 无效几何,分开统计 empty = gdf[gdf.is_empty] invalid = gdf[~gdf.is_valid] print(f"空几何 {len(empty)} 条, 无效几何 {len(invalid)} 条") # 无效的原因要看清楚,别一律 make_valid 糊过去 for idx in invalid.index[:5]: print(idx, shapely.is_valid_reason(gdf.geometry[idx])) # 面积必须在投影坐标系下算,先转 UTM 50N(研究区经度带大致落在这) gdf_m = gdf.to_crs(epsg=32650) gdf_m["area_m2"] = gdf_m.geometry.area print(gdf_m["area_m2"].describe()) # 多部件要素要单独揪出来,有些流程不认 MultiPolygon multi = gdf[gdf.geom_type == "MultiPolygon"] print("多部件要素:", len(multi))geom_type.value_counts()是第一个要看的,如果混着Polygon和MultiPolygon,后面很多算子行为会不一致。is_valid_reason返回的字符串会告诉你到底是Self-intersection、Ring Self-intersection还是Hole lies outside shell,不同病因对应不同修法,一律buffer(0)会把带洞的面悄悄填平。
面积那行,to_crs(epsg=32650)之后再.area,单位才是平方米。要是直接在地理坐标系上.area,GeoPandas 会给你一个“度平方”的数,还附一条 warning,很多人 warning 一关就当没看见。
3.3 转成 GeoJSON 或灌进 PostGIS
下游要么是前端可视化要 GeoJSON,要么是分析库要进 PostGIS。两条路都走一遍,心里有数。
GeoJSON 方向,主要控制体积和精度:
# 精度降到 6 位(约 0.1 米),体积能砍掉一半以上,区划数据完全够用 gdf_out = gdf.to_crs(epsg=4326) gdf_out.to_file( "D:/gis/out/zones.geojson", driver="GeoJSON", encoding="utf-8", # GeoJSON 规范要求 UTF-8,别写 GBK COORDINATE_PRECISION=6 )to_crs(4326)是因为 GeoJSON 规范(RFC 7946)规定用 WGS84 经纬度,虽然不少工具能容忍别的,但跨平台传的时候容易出幺蛾子。COORDINATE_PRECISION这个参数是驱动层透传的,不是所有版本都认,如果报unexpected keyword,改用gdf_out.geometry = gdf_out.geometry.simplify(0)之类的方式处理精度。
PostGIS 方向,命令行工具比 Python 更稳:
# -s 指定源坐标系EPSG,-I 建 GiST 空间索引,-W 指定源文件编码,-D 用 COPY 模式(比 INSERT 快一个量级) shp2pgsql -s 4326:4326 -I -W GBK -D \ D:/gis/hf_zone/zones.shp \ public.dim_zone | psql -h 127.0.0.1 -U gis -d analysis-s 4326:4326两个数字分别是源 SRID 和目标 SRID,写成一样就是不做转换,直接按声明入库。真要转投影就写-s 4326:4490这种。-D用 COPY 而不是逐条 INSERT,几千条要素的差距是几秒对几分钟。入库后记得确认 SRID 有没有真的写进geometry_columns,有些老版本工具链会漏。
4. 进业务系统前的改造:裁剪、拓扑修复与空间索引
4.1 按研究范围裁剪,别把碎片一起带走
原始区划范围通常比你的研究区大。直接按属性筛是错的——功能区边界跟研究区边界往往不重合,筛出来要么多一块要么少一块。正确做法是做几何求交。
import geopandas as gpd zones = gpd.read_file("D:/gis/hf_zone/zones.shp", encoding="GBK").to_crs(32650) study = gpd.read_file("D:/gis/study_area.shp").to_crs(32650) # 求交,keep_geom_type=True 保证结果只留面,避免退化成一维线 clipped = gpd.overlay(zones, study, how="intersection", keep_geom_type=True) # 求交后会产生大量碎块,按面积过滤,阈值按研究区尺度定 clipped["area_m2"] = clipped.geometry.area print(clipped["area_m2"].describe()) clipped = clipped[clipped["area_m2"] > 10000] # 丢掉小于 1 公顷的碎屑 # 同一片区被切碎成多块的,按主键合并回来 clipped = clipped.dissolve(by="qhdm", aggfunc="first").reset_index()keep_geom_type=True这个参数值得单独说。求交时,两个面如果只是边界相切,结果可能是一条线或者一个点,这种退化几何混进面图层,后面做面积统计会莫名其妙多出 0 值记录。dissolve那步是为了把被研究区边界切开的片区重新合成整块——不合成的话,一个区在结果里出现三四条记录,汇总时业务方会来问你为什么数量对不上。
碎片阈值不能拍脑袋,先describe()看一眼面积分布,通常能明显看到一撮几十到几百平方米的碎块和一个数量级之上的正常片区,取中间值就行。
4.2 拓扑检查:重叠、自相交、缝隙
区划数据最要命的问题不是单要素无效,是要素之间的关系错误。三类问题依次查:
重叠。两块面有公共区域,做面积汇总时这块地被算两遍,总人口、总销售额全部虚高。查法是对每个要素取其与研究区内其他要素的交集:
import itertools, shapely # 用 STRtree 建空间索引,只对候选对做精确判断 tree = shapely.STRtree(clipped.geometry) overlaps = [] for i, geom in enumerate(clipped.geometry): for j in tree.query(geom): if j <= i: continue inter = geom.intersection(clipped.geometry.iloc[j]) if not inter.is_empty and inter.area > 1.0: # 1 平方米以下当舍入误差 overlaps.append((clipped.index[i], clipped.index[j], inter.area)) print("重叠对数:", len(overlaps)) print(sorted(overlaps, key=lambda x: -x[2])[:5])tree.query()返回的是候选索引,还得做一次精确求交,索引只负责缩小范围。容差设 1 平方米是有意的——浮点运算总会有微小重叠,不放容差会报出几百对“重叠”,全是噪声。
自相交。单个面自己跟自己交叉,成因多半是矢量化时点错了位置。is_valid已经能查出来,修的时候用shapely.make_valid()而不是buffer(0):前者会保留多部件并尽量维持原形状,后者可能把环合并成一个怪异的外壳。
缝隙。相邻面之间漏了一条缝,这条缝不属于任何一个片区,落在里面的点位在汇总时会被静默丢弃。查法是unary_union之后跟整体外包框对比,或者直接看落点统计的总和是不是小于原始点位数。缝隙小的话可以接受,大缝隙就得回去找源数据。
提示:拓扑修复一定要留一版原始数据不动。修完的图层另存,字段里加一列标记修复方式,出了问题能回溯。我见过修完面积整体少了 3% 却没人记得改过什么的现场。
4.3 空间索引与投影面积计算
数据量小的时候感觉不到索引的存在,几万个点做落面查询时差距就出来了。GeoPandas 的.sindex是 STRtree,PostGIS 的 GiST 索引原理类似,都是先按外包框粗筛再精确判断。
# 点位落面,用 sindex 加速,比逐个 contains 快两三个数量级 points = gpd.read_file("D:/gis/poi.shp").to_crs(32650) joined = gpd.sjoin(points, clipped, how="inner", predicate="within") print("落面成功:", len(joined), "/", len(points))PostGIS 那边,入库之后确认索引建上了:
-- 确认空间索引存在,没有就补一个 CREATE INDEX IF NOT EXISTS idx_dim_zone_geom ON public.dim_zone USING GIST (geom); -- 面积统计,务必用 ST_Transform 到等积或 UTM 投影再算 SELECT zone_type, count(*) AS n, round(sum(ST_Area(ST_Transform(geom, 32650)))::numeric, 2) AS area_m2 FROM public.dim_zone GROUP BY zone_type ORDER BY area_m2 DESC;ST_Area作用在地理坐标系上返回的是“度平方”,这个坑每年都有人踩。要么ST_Transform到投影,要么用ST_Area(geom::geography)——后者按椭球面算,精度更高但更慢,数据量不大时够用。
面积字段还有一层要核对:属性表里自带的面积,跟你自己算的,差多少。差在 1% 以内是正常的投影差异,差到 20% 以上就是口径不同——很可能原始字段用的是椭球面积或者某种地方坐标系下的值。这个差异得在交付说明里写清楚,不然业务方拿两份报表一对账,先怀疑的是你的代码。
5. 排查清单:五个把我拖到凌晨的现场问题
5.1 属性表打开全是方块和问号
现象:QGIS 里几何正常,属性表的中文名一片锟斤拷或者???,导出 CSV 后更乱。
原因:.dbf本身没有编码自描述能力,靠.cpg文件声明,缺了它就按系统默认猜。Windows 下默认 GBK、Linux 下默认 UTF-8,同一份数据在两个人机器上打开结果不一样。
解决:先确认有没有.cpg,内容是什么。没有的话直接建一个,写一行GBK或UTF-8保存。读的时候在代码里显式指定encoding,QGIS 里在图层属性的“数据源”页手动设编码。已经读乱了的,重新读一遍就能恢复,别在乱码基础上做任何修改再保存。
5.2 图层叠上去偏了几百米,或者掉进几内亚湾
现象:区划图层跟底图对不上,偏移量固定且方向一致;或者整层图形缩成一个小点飘在坐标原点附近。
原因:前者是坐标系贴错了(本该是 CGCS2000 被当成 WGS84,或者反过来的地方坐标系),后者是.prj缺失被默认当成EPSG:4326处理。
解决:先用ogrinfo -so看有没有 SRS WKT。缺失的情况下,用set_crs(allow_override=True)直接赋值,不要用to_crs转换。偏移量固定的情况,量一下偏移多少米,跟已知的坐标系差异对照,能快速定位是哪一套。实在说不清就找数据提供方要投影参数,别自己试。
5.3 区划编码匹配率是 0,字段值看着还正常
现象:拿编码字段跟业务表 join,明明两边都有值,一条也对不上。
原因:编码字段在.dbf里被定义成数值型,前导零被吃掉;或者数据集内部就混用了123456和123456.0两种表示,一个是字符串一个是浮点转出来的。
解决:读取时用dtype={"qhdm": str}强制成字符串。已经进库的数据,用 SQL 补齐:lpad(code::text, 6, '0')。更稳妥的是入库时就建一张维表,编码做标准化处理,主键固定成文本类型,后面所有关联都走这张表。
5.4 汇总面积比研究区总面积还大
现象:各片区面积加总,比研究区实测总面积多出百分之几到十几。
原因:三个来源叠在一起——面之间有重叠;面积在未投影的坐标系上算的;属性表自带的面积字段跟几何对不上却被当成真值直接用。
解决:按 4.2 的流程查重叠,容差 1 平方米;统一转到投影坐标系再算;把自带面积和计算面积做成两列,差值超过阈值(比如 5%)的记录单独筛出来人工看。这三件事按顺序做,基本能定位到根因。
5.5 保存时报错,或者属性表被锁住改不了
现象:编辑属性表提示“无法写入”,或者另存为新的 shp 时报“文件过大”。
原因:.dbf被 Excel、WPS 或者另开的 QGIS 进程占用了;或者要素数/文件体积撞上 shapefile 格式的 2GB 上限(本质是 32 位偏移量限制)。
解决:先关掉所有可能打开.dbf的程序,包括那些后台留着的预览窗口。体积问题目测要撞线的话,直接换格式——GeoPackage(.gpkg)或者 FlatGeobuf 都是单文件、无大小限制、支持 UTF-8 字段名,QGIS 和 GDAL 生态全兼容。区划数据这种规模不大但字段名长的,用 GeoPackage 比死守 shapefile 舒服得多。
5.6 一套五分钟的排查顺序
上面五条如果同时出现两三条,别一条条试,按这个顺序走:
ogrinfo -so -al看要素数、几何类型、Extent、SRS WKT —— 30 秒,定位坐标系和文件完整性。- 打开属性表看字段类型和编码 —— 1 分钟,定位编码和前导零。
- QGIS 里按分类字段渲染一遍 —— 1 分钟,定位属性缺失和几何异常。
- 体检脚本跑
is_valid、is_empty、面积分布 —— 1 分钟,定位拓扑问题。 - 落面查询看匹配率 —— 2 分钟,定位重叠和缝隙。
这五步跑完还没头绪的情况很少。我踩过的坑里,八成在前两步就能看见苗头,只是当时急着往下写代码,跳过了。
6. 从底图到指标:点位落面汇总与三层交叉验证
底图整理干净之后,真正产生价值的动作是把业务数据挂上去。最常见的是点位落面:一批带经纬度的 POI、门店、事件记录,按落在哪个功能区归集,输出每个片区的数量、密度、人均指标。
这里有个细节决定结果对不对:predicate用什么。用within,落在边界上的点会被丢掉;用intersects,落在两个片区公共边上的点会被算两次,汇总总数直接超过原始行数。我一般用within,然后对未匹配上的点做一次缓冲区回捞,把边界点归到最近的那个片区:
import geopandas as gpd zones = gpd.read_file("D:/gis/hf_zone/zones.shp", encoding="GBK").to_crs(32650) poi = gpd.read_file("D:/gis/poi.shp").to_crs(32650) hit = gpd.sjoin(poi, zones[["qhdm", "zone_type", "geometry"]], how="inner", predicate="within") miss = poi.loc[~poi.index.isin(hit.index)] # 未匹配的点按最近邻归入片区,容差 50 米,超出的算真的不在范围内 if len(miss): near = gpd.sjoin_nearest(miss, zones[["qhdm", "geometry"]], how="left", max_distance=50) hit = gpd.pd.concat([hit, near], ignore_index=True) summary = hit.groupby(["qhdm", "zone_type"]).size().rename("poi_cnt").reset_index() print(summary.head())max_distance=50这个容差要按点位精度定,GPS 采的点 10 米够了,地址解析出来的点可能要到 200 米。容差设太大,会把本该属于隔壁片区的点吸过来,指标就失真了。
出完报表别急着交付,跑三层校验:
| 校验项 | 判断标准 | 不通过说明什么 |
|---|---|---|
| 点位总数守恒 | 汇总数 = 原始表行数 − 明确排除数 | 有静默丢弃,通常是 within 边界点或坐标不在范围内 |
| 面积总和 | 各片区面积和 ≈ 研究区总面积,偏差 < 1% | 面有重叠或缝隙,回到 4.2 查拓扑 |
| 抽样人工核对 | 随机抽 20 个点,落面结果可解释 | 坐标系或边界数据本身有问题 |
三层里最容易被跳过的是第三层,也恰恰是最能救命的一层。数字守恒、面积对得上,但坐标系整体偏移 50 米的话,前两层照样通过——因为偏移是均匀的,所有点一起错。抽 20 个点,把经纬度贴到地图上看一眼落在哪条街上,五分钟的事。
从那以后我每次交付区划汇总类的结果,都强制先把这三层校验跑一遍再发出去,哪怕需求方只想要一个 Excel。图好看、数自洽、抽样能解释,这三条同时成立,才敢说这份数据能用。希望帮到你。
本文还有配套的精品资源,点击获取