简介本资源为厦漳泉地区矢量边界数据包面向从事GIS分析、城市规划、区域经济研究及地图制图的人员可用于行政区划可视化、空间统计与专题图制作。压缩包共20个文件约4.16MB以shp、shx、dbf、prj等Shapefile核心文件为主辅以sbx、sbn、cpg、xml、evf等索引与元数据文件覆盖厦门、漳州、泉州各区县行政区域边界可直接在ArcGIS中加载使用。资源包含Export_Output与ROI等多组图层便于按需提取与叠加分析。目前已有229人学习下载适合需要快速获取厦漳泉行政边界底图、开展区域空间研究的用户参考使用。1. 厦漳泉矢量边界.zip一份区域地理数据包到底能拿来做什么拿到「厦漳泉矢量边界.zip」这个文件名多数人的第一反应是解压看看里面有几个 shp。但真正做过区域空间分析的人会先问三个问题这份边界是行政区划面还是海岸线坐标系是 WGS84 还是 CGCS2000属性表里有没有区县编码这三个问题决定了它能不能直接进你的分析流程还是得先做一轮清洗。厦漳泉指的是福建南部厦门、漳州、泉州三市构成的经济协作区域做区域经济、物流可达性、人口热力、选址评估的人经常需要这一片的多边形底图。这份数据包的价值不在于「有边界」而在于它把三市及其下辖区县拼成一套可叠加、可裁剪、可做空间统计的矢量面。适合谁用做区域 GIS 分析的学生和工程师、需要底图做可视化的数据岗、以及要给业务系统补一层行政区划空间维度的后端开发。下面按「先搞清数据长什么样再动手跑通最后避开那些让人返工的坑」的顺序讲。2. 先看清数据矢量边界的格式、坐标系与属性结构2.1 解压后先做一次数据体检拿到压缩包别急着往 QGIS 里拖先看文件构成。常见的矢量边界包一般包含.shp、.shx、.dbf、.prj、.cpg这几个同名文件少一个都可能出问题。.shp存几何.shx是索引.dbf存属性表.prj定义坐标系.cpg声明字符编码。用命令行先列一遍# 解压到独立目录避免和别的数据混在一起 unzip 厦漳泉矢量边界.zip -d ./xzq_boundary # 查看文件清单和大小确认 shp 配套文件齐全 ls -lh ./xzq_boundary逻辑说明单独建目录是为了后面用脚本批量处理时不误伤其他数据。ls -lh看的是文件大小如果.shp只有几 KB 而.dbf却很大说明几何可能被简化过或者属性字段特别多这会影响后续叠加分析的精度。参数上没什么可调的重点是确认.prj存在——没有它坐标系就是黑匣子后面做面积计算和距离分析全是错的。2.2 坐标系判断为什么不能直接算面积矢量边界最常见的两种坐标系是地理坐标系如 WGS84、CGCS2000和投影坐标系如 UTM、高斯克吕格。地理坐标系的单位是度直接拿它算面积得到的是「平方度」没有任何物理意义。判断方法很简单看.prj文件内容或者用 GDAL 查# 用 ogrinfo 查看图层坐标系和字段结构 ogrinfo -so -al ./xzq_boundary/xzq_boundary.shp输出里Geometry:后面跟的是几何类型Feature Count:是要素数量Extent:是范围最关键的是Layer SRS WKT:那一段。如果看到GEOGCS开头说明是地理坐标系做面积统计前必须投影。我一般会先转成适合该区域的投影坐标系再算厦漳泉大致在东经 117° 到 119°、北纬 23° 到 25° 之间用 CGCS2000 3 度带高斯克吕格投影比较合适中央经线取 118° 或 119° 都行误差在可接受范围。提示如果.prj缺失或内容为空不要猜坐标系。先拿一个已知坐标点比如某市政府大楼的经纬度去比对几何范围确认后再补写.prj否则后面所有空间运算都是错的。2.3 属性表里该有什么区县编码与名称字段一份能用的行政区划边界属性表至少要有行政区代码和名称两个字段。代码通常是 6 位或 12 位前 6 位对应区县前 4 位对应地市。用 Python 快速看一眼字段import geopandas as gpd # 读取矢量边界注意指定编码中文属性容易乱码 gdf gpd.read_file(./xzq_boundary/xzq_boundary.shp, encodingutf-8) # 打印字段名、几何类型、坐标系和要素数量 print(字段:, list(gdf.columns)) print(几何类型:, gdf.geom_type.unique()) print(坐标系:, gdf.crs) print(要素数:, len(gdf)) # 查看前几行属性确认区县名称和代码是否完整 print(gdf.drop(columnsgeometry).head())逻辑说明encodingutf-8是为了防止中文区县名变成乱码如果读出来是问号或方块换成gbk再试。geom_type应该全是Polygon或MultiPolygon如果混进了LineString或Point说明数据包里的边界不纯需要按几何类型拆开处理。crs输出如果是EPSG:4326这类地理坐标系后面算面积前记得投影。参数上drop(columnsgeometry)只是为了让打印结果干净不影响数据本身。3. 把边界跑起来从裁剪到空间统计的最小工作流3.1 用边界裁剪栅格数据以夜间灯光或土地利用为例拿到边界后最常做的一件事是裁剪。比如你有一份全国范围的栅格数据只想保留厦漳泉区域。用rasterio配合geopandas可以一步到位import geopandas as gpd import rasterio from rasterio.mask import mask # 读取边界并统一到与栅格一致的坐标系 boundary gpd.read_file(./xzq_boundary/xzq_boundary.shp, encodingutf-8) with rasterio.open(./data/nightlight.tif) as src: # 如果边界和栅格坐标系不一致先投影边界 if boundary.crs ! src.crs: boundary boundary.to_crs(src.crs) # 用边界几何裁剪栅格cropTrue 会同时缩小输出范围 out_image, out_transform mask(src, boundary.geometry, cropTrue) out_meta src.meta.copy() out_meta.update({ height: out_image.shape[1], width: out_image.shape[2], transform: out_transform }) with rasterio.open(./output/xzq_nightlight.tif, w, **out_meta) as dest: dest.write(out_image)逻辑说明mask函数的核心参数是shapes和crop。shapes传的是几何列表这里直接给boundary.geometry它会自动处理多个多边形。cropTrue会把输出栅格的范围收紧到边界外接矩形减少无效像元。out_meta必须更新height、width和transform否则写出的文件坐标会错位。参数上如果边界是 MultiPolygon 且内部有空洞mask默认会按外环裁剪需要保留空洞的话得用all_touchedFalse并检查几何有效性。3.2 按区县做空间聚合统计每个单元的指标裁剪完栅格下一步往往是按区县统计均值或总和。思路是把栅格转成点或直接做分区统计。用rasterstats更省事from rasterstats import zonal_stats # 对每个区县多边形统计栅格均值、最大值和像元数 stats zonal_stats( ./xzq_boundary/xzq_boundary.shp, ./output/xzq_nightlight.tif, stats[mean, max, count], geojson_outFalse, nodata-9999 ) # 把统计结果写回属性表 import geopandas as gpd gdf gpd.read_file(./xzq_boundary/xzq_boundary.shp, encodingutf-8) gdf[mean_light] [s[mean] for s in stats] gdf[max_light] [s[max] for s in stats] gdf.to_file(./output/xzq_stats.shp, encodingutf-8)逻辑说明zonal_stats第一个参数是矢量路径第二个是栅格路径stats列表里mean是均值、max是最大值、count是有效像元数。nodata-9999要跟栅格本身的 nodata 值一致否则会把无效值算进去。参数上如果区县边界和栅格像元边界不完全对齐默认按像元中心点判断归属想要更精确可以加all_touchedTrue但会稍微慢一点。写回属性表时注意字段名不要超过 10 个字符这是 shapefile 的老限制。3.3 边界叠加与拓扑检查别让自相交毁掉分析矢量边界用久了难免遇到自相交、重叠、缝隙。做叠加分析前先跑一次拓扑检查import geopandas as gpd from shapely.validation import explain_validity gdf gpd.read_file(./xzq_boundary/xzq_boundary.shp, encodingutf-8) # 检查每个几何是否有效打印无效原因 for idx, geom in enumerate(gdf.geometry): if not geom.is_valid: print(f要素 {idx} 无效: {explain_validity(geom)}) # 用 buffer(0) 修复大部分自相交问题 gdf[geometry] gdf.geometry.apply(lambda g: g.buffer(0) if not g.is_valid else g) gdf.to_file(./output/xzq_fixed.shp, encodingutf-8)逻辑说明explain_validity会告诉你无效原因是「Self-intersection」还是「Ring Self-intersection」等。buffer(0)是修复自相交的常用技巧它会把多边形重新构造成有效几何但可能轻微改变边界形状所以修复后要对比面积变化。参数上如果buffer(0)后仍有无效几何说明问题比较严重需要手动编辑或换用make_valid方法。4. 避坑指南矢量边界处理中最容易翻车的五个地方4.1 中文属性乱码现象是区县名变问号原因是编码没对上现象用 QGIS 或 Python 打开后属性表里的区县名称显示为???或乱码方块。原因shapefile 的.dbf文件默认编码可能是GBK或GB2312而读取时用了UTF-8。解决先用encodinggbk读一遍确认能正常显示后再统一转成 UTF-8 保存。如果.cpg文件存在里面写的就是编码声明直接按它来。4.2 面积算出来是负数或极小值现象是统计结果离谱原因是坐标系没投影现象用gdf.area算面积结果要么是零点几要么是几万跟实际平方公里对不上。原因地理坐标系下area返回的是平方度不是平方米。解决先to_crs转到投影坐标系再算。厦漳泉区域可以用EPSG:4547或EPSG:4548这类 CGCS2000 3 度带投影转完再算面积就是平方米除以 10^6 得平方公里。4.3 裁剪后栅格全黑或全透明现象是输出没有数据原因是边界和栅格坐标系不一致现象用边界裁剪栅格后输出文件打开全是 nodata。原因边界是 WGS84栅格是投影坐标系mask时几何位置对不上。解决裁剪前先判断boundary.crs ! src.crs不一致就to_crs转换。另外检查边界范围是否真的和栅格有交集用gdf.total_bounds和src.bounds对比一下。4.4 叠加分析报拓扑异常现象是intersection报错原因是几何自相交现象做gdf.overlay或intersection时抛出TopologyException。原因边界多边形存在自相交或环方向错误。解决先跑is_valid检查用buffer(0)或make_valid修复。修复后重新检查面积变化如果变化超过 1%说明原始数据质量问题较大建议换数据源。4.5 属性表字段名被截断现象是字段名超过 10 字符后丢失原因是 shapefile 格式限制现象写回 shapefile 时字段名mean_light变成mean_ligh。原因shapefile 的.dbf字段名最多 10 个字符。解决要么改用 GeoPackage 或 GeoJSON 格式保存要么把字段名缩到 10 字符以内。如果必须用 shapefile建议用gdf.rename(columns{mean_light: mn_light})先改短。5. 进阶用法用边界做空间索引与可视化底图5.1 建立空间索引加速查询当边界要素多、查询频繁时空间索引能显著提速。GeoPandas 底层用rtree但手动建索引更可控import geopandas as gpd from shapely.geometry import Point gdf gpd.read_file(./xzq_boundary/xzq_boundary.shp, encodingutf-8) # 建立空间索引 sindex gdf.sindex # 查询某个点落在哪个区县 point Point(118.1, 24.5) # 先用索引粗筛再做精确判断 candidates list(sindex.intersection(point.bounds)) for idx in candidates: if gdf.geometry.iloc[idx].contains(point): print(gdf.iloc[idx][name]) break逻辑说明sindex.intersection传的是点的外接矩形返回可能相交的要素索引列表然后再用contains做精确判断。这样比遍历所有要素快一个数量级。参数上point.bounds返回(minx, miny, maxx, maxy)对于点来说四个值相同但接口要求这样传。5.2 用边界做分级设色底图可视化时按区县指标做分级设色很常见。用matplotlib配合geopandas可以快速出图import matplotlib.pyplot as plt gdf gpd.read_file(./output/xzq_stats.shp, encodingutf-8) # 按均值灯光指数分级设色cmap 用 YlOrRd 表示亮度 fig, ax plt.subplots(figsize(10, 8)) gdf.plot(columnmean_light, cmapYlOrRd, legendTrue, edgecolorgray, linewidth0.5, axax) ax.set_title(厦漳泉各区县夜间灯光均值) ax.axis(off) plt.savefig(./output/xzq_choropleth.png, dpi300, bbox_inchestight)逻辑说明column指定要映射的字段cmap是色带legendTrue会加图例。edgecolor和linewidth控制边界线样式太粗会盖住颜色。axis(off)去掉坐标轴出图更干净。参数上dpi300适合印刷屏幕展示 150 就够。如果字段有缺失值plot会默认留白可以用missing_kwds指定缺失颜色。5.3 一个我常犯的错忘了统一坐标系就做叠加早期做区域分析时我习惯性把边界和业务数据都读进来直接overlay结果要么报错要么结果偏移几百米。后来养成习惯任何两个图层做空间运算前先打印各自的crs不一致就统一到其中一个。这个习惯帮我省掉了大量排查时间。另一个教训是修复几何后一定要重新算面积buffer(0)有时会把小多边形吃掉面积变化超过 5% 就得警惕。希望帮到你。本文还有配套的精品资源点击获取