简介：这份资源提供2021年全国村名点矢量数据，面向地理信息、城乡规划、区域经济与空间分析学习者，以及需要村庄点位底图开展地图制图、分布统计和空间核验的从业者。压缩包共6个文件、约88.73MB，主体为shp矢量数据，并配有prj坐标系文件、dbf属性表、shx索引、xml元数据及txt说明，便于在ArcGIS、QGIS等软件中直接加载、查看字段与坐标信息。已有397人学习下载。读者可用于全国村庄点位分布展示、专题地图制作、行政区或研究区筛选、空间分析与数据比对，也能借助属性表进行名称检索、分类统计和可视化表达，适合作为GIS课程练习、论文制图与项目前期底图素材；使用前建议结合说明文件核对坐标系、字段含义与数据年份口径，按研究范围裁剪后再开展分析，并按需转换投影。
1. 全国村名点数据不是一张表：先搞清它是点位库还是地名库
不少人第一次搜「2021年全国村名点数据」，脑子里浮现的是一份大 Excel：一列省、一列市、一列村名，翻到底就完事。真拿到手才会发现，它更像一套点位库——每条记录对应一个村的地理位置，外加行政区划归属和层级代码。这两者差别不小：名录只能做检索和对照，点位才能落到地图上算距离、算覆盖、做空间连接。做末端网点规划、县域商业摸底、乡村服务设施可达性分析的人，卡点几乎都卡在同一处——县以上有边界，村一级没坐标。这份资源补的正是村这一级的空间落点，适合要在 QGIS 里叠边界出图的、要把村名和业务表按代码关联的、以及要给几十万个点做聚合统计的从业者。它不适合指望拿它当村级行政边界用的人，点不是面，边界得另找。
2. 拆包先看三件事：文件构成、字段语义、坐标系
拿到一份村点数据，先别急着打开属性表看村名。顺序应该是：先看目录里有哪些文件，再看表头有哪些列，最后确认坐标系。这个顺序反了，后面每一步都在返工。我见过太多次直接双击 .shp 打开、结果中文全是问号、点位整体飘到海里的情况，问题无一例外都出在这三件事上。这一章把这三件事拆开讲清楚，后面两章的代码才有落脚点。
2.1 目录结构与文件格式：shp 从来不是单个文件
Shapefile 是一组文件共用同一个主名的集合，少一个都可能出问题。常见的配套关系是这样的：
| 扩展名 | 作用 | 缺失后果 |
|---|---|---|
| .shp | 存几何（点/线/面） | 直接打不开 |
| .dbf | 存属性表 | 只剩图形，没有村名和代码 |
| .shx | 几何索引 | 部分软件能打开但定位异常 |
| .prj | 坐标系声明 | 坐标变成"裸数字"，无法判断投影 |
| .cpg | 属性表编码声明 | 中文村名乱码 |
先看一眼目录，确认文件齐不齐：
# 列出目录，重点看 .prj 和 .cpg 在不在，以及总大小是否合理 ls -lh village_2021/ # 直接读 .prj 的文本内容，能看出是地理坐标系还是投影坐标系 cat village_2021/*.prj # 看属性表前几行，确认编码和分隔符 file village_2021/*.dbf.prj里如果出现GEOGCS，说明是经纬度地理坐标系；出现PROJCS，说明是投影坐标系（单位通常是米）。这一步只用三十秒，但能省掉后面半小时的排查。
除了 Shapefile，这套数据还可能是 GeoJSON 或 CSV。三种形态各有取舍，选哪种取决于你后面要干什么：
| 形态 | 优点 | 局限 | 适合的场景 |
|---|---|---|---|
| Shapefile | GIS 软件通用性最好 | 多文件、字段名限 10 字符、编码易乱 | QGIS/ArcGIS 直接打开 |
| GeoJSON | 单文件、UTF-8、Web 友好 | 体积大、无字段类型约束 | 前端地图、接口传输 |
| CSV + 经纬度列 | 最轻、易 diff、易入库 | 无几何、无坐标系声明 | 批量清洗、进数据库 |
如果拿到的是 CSV，一定要确认里面有没有单独的经度列和纬度列。有些包只给了「中心点坐标」一个合并字段，那就得先切分，别硬塞。
2.2 字段表：哪几列能直接拿去 join
村点数据的字段不算多，但每一列的可靠性差别很大。下面这张表是我实际处理时会给同事的说明，重点看最后一列：
| 字段 | 类型 | 语义 | 能不能当关联键 |
|---|---|---|---|
| prov_name | 字符 | 省级名称 | 不建议，名称有别名和简称 |
| city_name | 字符 | 地级名称 | 不建议 |
| county_name | 字符 | 县级名称 | 不建议 |
| town_name | 字符 | 乡镇名称 | 可做辅助校验 |
| village_name | 字符 | 村名 | 只做展示，重名率极高 |
| adcode | 字符(12) | 行政区划代码 | 首选关联键 |
| urban_rural | 字符(3) | 城乡分类代码 | 做筛选维度 |
| lon / lat | 数值 | 经纬度 | 用来建点 |
adcode是这份数据里最值钱的一列。它稳定、无歧义、可排序、可截断，前 6 位就是县级代码，拿去和任何一份统计年鉴、行政区划表都能对上。村名则完全不同，「王庄」「李村」「张家湾」这类名字在全国范围内的重名数量能到四位数，用村名做 join 键等于给自己埋雷。
城乡分类代码值得单独说一句。这类代码通常是三位，首位为 1 表示城镇属性，首位为 2 表示乡村属性，后面两位细分到主城区、城乡结合区、镇中心区、乡中心区、村庄等。做网点规划时，用这个字段一刀切出「乡村点位」，比按村名猜要靠谱得多。
2.3 坐标系不是玄学：.prj 和 crs 属性里写着答案
坐标系这件事，说玄学也玄学，说清楚也就三行代码。用 geopandas 读进来直接打印：
import geopandas as gpd # encoding 先试 gbk，中文乱码再换 utf-8 gdf = gpd.read_file("village_2021.shp", encoding="gbk") # 打印完整坐标系定义 print(gdf.crs) # 打印 EPSG 编码，None 表示没有 .prj 或无法识别 print(gdf.crs.to_epsg() if gdf.crs else None)如果返回EPSG:4326，说明是 WGS84 经纬度；返回EPSG:4490，说明是 CGCS2000 经纬度；这两个在村点这个精度上差别只有厘米级，可以视为等价。如果返回None，那就是.prj丢了，这时候只能靠经验判断：把点直接丢到在线底图上，看整体偏移方向和幅度。
偏移这件事有规律可循。整体往东北方向偏几百米、且偏移量随纬度变化，基本可以判定原始坐标是 GCJ-02（国内互联网地图常用的加密偏移坐标系）；偏移量更大且方向不规则，可能是 BD-09。这两种都不能靠to_crs解决，to_crs只做基准面之间的数学变换，处理不了人为加偏。这一点后面第 3 章会给出具体打法。
2.4 行政区划代码的 12 位怎么读
12 位代码不是一串随机数字，它是分层级的：
1 1 0 0 0 0 0 0 0 0 0 0 └─┬─┘ └─┬─┘ └─┬─┘ └─┬─┘ 省级 市级 县级 乡村级 2位 2位 2位 6位前 6 位是县级单位代码，第 7 到 9 位是乡镇级单位代码，第 10 到 12 位是村级单位代码。实际使用中，最常用的是截取前 6 位做县级汇总、截取前 9 位做乡镇汇总。
2021 年这个时间点需要特别注意。这一年前后有一轮较大规模的区划调整，撤县设区、合村并居都有发生，部分村级代码被废止或合并。如果你手上的边界数据是近两年的，而村点是 2021 年的，直接按代码 join 会出现一批匹配不上的记录，这不是数据坏了，是版本错位。第 4 章会讲怎么兜底。
3. 从点表到可用底图：读入、纠偏、空间连接的完整链路
数据拆明白之后，接下来是一套固定的处理链路：读入并把编码理顺，判断坐标系并做必要转换，和行政区划边界做空间连接补全归属，最后导出一份同事能直接打开的成果。这条链路我一般写成四个函数，每步都能单独跑、单独验证，出错时不用从头再来。
3.1 读入与编码：先把中文理顺再谈别的
第一步永远是编码。Shapefile 的属性表编码由.cpg声明，但很多分发出来的包压根没带这个文件，pandas 和 geopandas 就会按默认编码解析，中文村名直接变成问号或者一串乱码。
import pandas as pd import geopandas as gpd from shapely.geometry import Point # 路线一：从 Shapefile 读，指定编码 gdf = gpd.read_file( "village_2021.shp", encoding="gbk", # 乱码时优先试 gbk，其次 utf-8、utf-8-sig columns=[ # 只读需要的列，大文件能省一半内存 "prov_name", "city_name", "county_name", "town_name", "village_name", "adcode", "urban_rural", "lon", "lat" ] ) # 路线二：从 CSV 读，关键是把 adcode 强制读成字符串 df = pd.read_csv( "village_2021.csv", dtype={"adcode": str}, # 不加这句，前导零会被吃掉 encoding="utf-8-sig" # 带 BOM 的 CSV 用这个，避免第一列名带不可见字符 ) # 从经纬度列构建几何，并声明源坐标系 gdf = gpd.GeoDataFrame( df, geometry=gpd.points_from_xy(df["lon"], df["lat"]), crs="EPSG:4326" )dtype={"adcode": str}这一行是这套流程里最容易被忽略、后果最严重的一行。不加它，pandas 会把 12 位代码当整数读，前导零丢掉，后面的所有关联全错。columns参数则是性能优化，全国村点动辄几十万行，少读几列能明显缩短时间。points_from_xy的顺序是先 x 后 y，也就是先经度后纬度，写反了整个数据集会跑到南极附近，这是新手最常见的翻车点之一。
读完之后马上做一次完整性检查，别急着往下走：
# 三个必查项：空几何、空代码、坐标范围 print("空几何数量:", gdf.geometry.isna().sum()) print("代码为空数量:", gdf["adcode"].isna().sum()) print("经度范围:", gdf["lon"].min(), gdf["lon"].max()) print("纬度范围:", gdf["lat"].min(), gdf["lat"].max())经纬度范围能帮你快速识别坐标写反的情况。正常的全国范围，经度大致在 70 到 140 之间，纬度在 3 到 55 之间。如果经度出现 20 到 50 这样的数值，那基本可以确定 lon 和 lat 两列被搞反了，交换一下重跑就行。
3.2 坐标转换与纠偏：三种情况三套打法
坐标系处理要分情况，不能一套代码打天下。
情况 A，源数据本来就是 WGS84 或 CGCS2000。这种情况最省事，直接投影到你要用的坐标系即可。做距离计算时不要在 4326 上直接算，单位是度，算出来的数没有物理意义。全国范围我一般用 Albers 等积投影：
# 全国范围常用的 Albers 等积投影，标准纬线 25/47，中央经线 105 albers = "+proj=aea +lat_1=25 +lat_2=47 +lat_0=0 +lon_0=105 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs" gdf_m = gdf.to_crs(albers) # 此时 gdf_m.geometry.distance() 的单位是米情况 B，源数据是 GCJ-02。这种数据的特点是：直接叠在 WGS84 底图上，整体往东北偏，偏移量在几百米量级。to_crs处理不了，必须用专门的纠偏算法。常见做法是找一个开源的转换脚本，核心是一个分段函数加椭球参数：
import math # 椭球参数，长半轴与第一偏心率平方 A = 6378245.0 EE = 0.00669342162296594323 def _transform_lat(x, y): ret = -100.0 + 2.0 * x + 3.0 * y + 0.2 * y * y + 0.1 * x * y + 0.2 * math.sqrt(abs(x)) ret += (20.0 * math.sin(6.0 * x * math.pi) + 20.0 * math.sin(2.0 * x * math.pi)) * 2.0 / 3.0 ret += (20.0 * math.sin(y * math.pi) + 40.0 * math.sin(y / 3.0 * math.pi)) * 2.0 / 3.0 ret += (160.0 * math.sin(y / 12.0 * math.pi) + 320 * math.sin(y * math.pi / 30.0)) * 2.0 / 3.0 return ret def _transform_lng(x, y): ret = 300.0 + x + 2.0 * y + 0.1 * x * x + 0.1 * x * y + 0.1 * math.sqrt(abs(x)) ret += (20.0 * math.sin(6.0 * x * math.pi) + 20.0 * math.sin(2.0 * x * math.pi)) * 2.0 / 3.0 ret += (20.0 * math.sin(x * math.pi) + 40.0 * math.sin(x / 3.0 * math.pi)) * 2.0 / 3.0 ret += (150.0 * math.sin(x / 12.0 * math.pi) + 300.0 * math.sin(x / 30.0 * math.pi)) * 2.0 / 3.0 return ret def gcj02_to_wgs84(lng, lat): """输入 GCJ-02 经纬度，输出 WGS84 经纬度，误差在米级""" dlat = _transform_lat(lng - 105.0, lat - 35.0) dlng = _transform_lng(lng - 105.0, lat - 35.0) radlat = lat / 180.0 * math.pi magic = math.sin(radlat) magic = 1 - EE * magic * magic sqrtmagic = math.sqrt(magic) dlat = (dlat * 180.0) / ((A * (1 - EE)) / (magic * sqrtmagic) * math.pi) dlng = (dlng * 180.0) / (A / sqrtmagic * math.cos(radlat) * math.pi) return lng - dlng, lat - dlat这个函数对全国村点逐条跑会比较慢，几十万行大概需要几分钟。实际处理时可以用向量化版本，或者先按县级分区、抽样验证纠偏效果再批量跑。判断要不要纠偏有个土办法：随便抽五个点，把坐标贴到带卫星影像的地图上，看落在村子中间还是落在村外的田里。偏几百米和落在村中心，肉眼一眼能分出来。
情况 C，源数据是 BD-09。需要先转到 GCJ-02 再转到 WGS84，两步走，顺序不能反。
3.3 与行政区划边界做空间连接
村点自带省市区县名称，但这些名称来自属性表，不一定和你的边界数据完全一致。做统计汇总之前，用空间关系补一遍归属更稳：
# 边界数据读入后先统一到同一坐标系 boundary = gpd.read_file("county_boundary.shp", encoding="gbk").to_crs(albers) # 村点也投影过去，两边坐标系必须一致才能做空间关系判断 pts = gdf.to_crs(albers) # 左连接，保留所有村点，未匹配的用 NaN 标出 joined = gpd.sjoin( pts, boundary[["adcode", "geometry"]].rename(columns={"adcode": "adcode_bnd"}), how="left", predicate="within" # within 表示点落在面内；边界上的点用 intersects ) # 统计未匹配数量，这批点要么在边界外，要么边界数据有空洞 unmatched = joined[joined["adcode_bnd"].isna()] print("未匹配数量:", len(unmatched), "占比:", len(unmatched) / len(joined))predicate参数有两个常用值：within严格判断点是否在面内部，落在边界线上的点会被判为不匹配；intersects则包含边界线。全国尺度的数据，用intersects更稳妥，未匹配率能低一些。how="left"保证村点一条不丢，这样你才能看到未匹配的那部分长什么样。
未匹配率超过 1% 就要停下来查。常见原因有三个：边界数据本身是多边形拼接有缝隙；村点坐标系还没纠偏，整体偏出去了；边界数据的行政区划版本和村点不一致。
3.4 导出一份同事能直接打开的成果
导出格式我一般选 GeoPackage，不用 Shapefile。原因很实际：Shapefile 字段名限 10 个字符，adcode_bnd这种名字会被截断成adcode_bn，交付给不熟悉的人就是个黑匣子；编码问题也反复出现。GeoPackage 是单文件、UTF-8、字段名不截断，QGIS 和 ArcGIS 都能直接打开。
# 只保留交付需要的列，避免带一堆中间字段出去 out = joined[[ "prov_name", "city_name", "county_name", "town_name", "village_name", "adcode", "adcode_bnd", "urban_rural", "geometry" ]].rename(columns={"adcode_bnd": "bnd_code"}) # 写 GeoPackage，同时建空间索引，QGIS 打开会快很多 out.to_file( "village_2021_clean.gpkg", layer="village", driver="GPKG", spatial_index=True )spatial_index=True这个参数在小数据上看不出差别，数据量上到十万级，QGIS 里拖动地图的流畅度完全是两个体验。交付前顺手在 QGIS 里加载一次，按urban_rural做个分类符号化，确认点的分布形态符合常识——东部密集、西部稀疏，这个宏观格局不对，说明前面某一步出了系统性问题。
4. 五个真会翻车的点：坐标偏移、代码丢零、重名、乱码、版本漂移
这一章是我自己在处理村点数据时踩过的坑，按出现频率从高到低排。每一条都按「现象 → 原因 → 解决」写，遇到对应症状可以直接对号入座。
4.1 现象：所有点整齐地往一个方向偏几百米
现象很典型：把村点叠到卫星影像上，所有点不是随机散落，而是整体朝同一个方向平移，偏移量在不同纬度略有差异。
原因只有一个：源数据是 GCJ-02 或 BD-09，被当成 WGS84 用了。这类坐标系是国内互联网地图服务为了合规做的非线性加偏，偏移量不是常数，随经纬度变化，所以不能靠加个固定偏移量修正。
解决分三步。第一，确认来源：如果这份数据是从某个在线地图接口导出的，那几乎肯定是加密坐标系；如果是测绘部门或公开的区划数据，一般是 WGS84 或 CGCS2000。第二，跑 3.2 节的纠偏函数。第三，纠偏后必须抽检，随机挑 5 到 10 个点贴到影像上，确认落在村委会、村口或者聚落中心，而不是落在村外的农田里。这一步不能省，纠偏函数写错一个符号，结果照样是偏的，只是偏的方向变了。
4.2 现象：Excel 双击打开，12 位代码变成 4 位或者科学计数法
现象是adcode列显示成1.10105E+11，或者前面的零没了，变成 10 位。
原因很明确：Excel 按数值类型解析 CSV，超过 11 位的整数会转科学计数法，带前导零的会直接丢零。
解决的关键是「别用 Excel 中转」。整个流程里 CSV 只应该被 pandas 读，不应该被人打开再另存。如果非要给人看，导出时把代码列加个前缀：
# 导出给人工查看的版本，代码列加制表符前缀强制成文本 view = gdf[["prov_name", "county_name", "village_name", "adcode"]].copy() view["adcode"] = "\t" + view["adcode"].astype(str) # 制表符让 Excel 按文本处理 view.to_csv("village_view.csv", index=False, encoding="utf-8-sig")程序读取的那份则始终保证dtype={"adcode": str}。两套文件分开，互不干扰。
4.3 现象：按村名关联，一条业务数据匹配出几十条
现象是把业务表和村点按village_name做 join，结果行数从几万膨胀到几十万，一条记录对应十几个村。
原因是村名重名率远超直觉。同一个村名在不同省份、不同县、甚至同一个县的不同乡镇都会重复出现。用村名当键，等于放弃唯一性约束。
解决方式是建立复合键。首选adcode12 位全码，唯一性最好。如果业务表只有村名没有代码，退而求其次用「县级代码前 6 位 + 乡镇名 + 村名」三元组，匹配率会降到可接受范围。实在什么代码都没有，就只能靠空间位置反查：把业务地址地理编码成点，再和村点做最近邻匹配，取最近的那个村。这条路误差最大，只适合做粗略归属，不能用于统计口径。
4.4 现象：属性表里中文村名全是问号
现象是打开属性表，village_name一列全是???，其他数字字段正常。
原因是.dbf的编码和读取时指定的编码不一致。.cpg文件如果写的是UTF-8而实际内容是 GBK，或者.cpg压根不存在，读取端就会猜错。
解决有两条路。第一条，在读取时显式指定encoding="gbk"，大多数国内分发的 Shapefile 都是 GBK。第二条，用 GDAL 的命令行工具直接转一次编码，转完再读：
# 把 Shapefile 转成 GeoPackage，同时强制按 UTF-8 写入属性表 ogr2ogr -f "GPKG" village_utf8.gpkg village_2021.shp \ -lco ENCODING=UTF-8 \ --config SHAPE_ENCODING GBK--config SHAPE_ENCODING GBK告诉 GDAL 源文件按 GBK 解析，-lco ENCODING=UTF-8指定输出按 UTF-8 写。这一条命令能同时解决乱码和格式转换两个问题，比在代码里反复试编码省事。
4.5 现象：用近两年的边界数据接 2021 年的村点，大批匹配不上
现象是空间连接之后未匹配率偏高，或者按代码 join 时有一批记录对不上，但村名看起来又是对的。
原因是行政区划在这几年有调整，撤县设区会改变县级代码，合村并居会废止部分村级代码。2021 年的村点用近两年的边界去套，版本错位是必然的。
解决办法是准备一张新旧代码对照表。匹配顺序分三级：第一级，12 位全码精确匹配；第二级，前 9 位乡镇码匹配，允许村级码不同；第三级，按村名加空间位置兜底，取距离最近且名称一致的记录。三级匹配都失败的，单独拉一张表出来人工看，通常数量不大，几十到几百条，逐个确认比写复杂规则划算。
提示：版本漂移这件事没有一劳永逸的解法，只能在流程里留一个「未匹配清单」的出口。每次处理完都看一眼这个清单的规模和分布，规模突然变大就是上游数据换了版本。
5. 把村点接进分析：覆盖半径、网格聚合、叠加路网
数据洗干净只是起点，真正体现价值的是把它接进具体的分析任务。这一章给三个我实际用过、复现成本不高的场景：算服务半径覆盖、按网格做聚合、和路网叠加修正直线距离。三个场景的共同前提是，距离计算必须在投影坐标系下做，不能在经纬度上直接算。
5.1 服务半径覆盖：每个村到最近设施有多远
这是村点数据最直接的用法。手上有若干设施点（卫生室、快递驿站、农资网点都行），要算每个村到最近设施的距离，输出一张覆盖情况表。
import geopandas as gpd import numpy as np from shapely.ops import nearest_points # 两边都投影到同一等积坐标系，单位统一为米 villages = gpd.read_file("village_2021_clean.gpkg").to_crs(albers) facilities = gpd.read_file("facilities.gpkg").to_crs(albers) # 建空间索引，几十万点对几千设施，不建索引会跑到天荒地老 fac_sindex = facilities.sindex def nearest_info(point): """返回最近设施的距离（米）和设施索引""" idx = list(fac_sindex.nearest(point, 1))[0] if len(idx) == 0: return np.nan, None fidx = idx[0] dist = point.distance(facilities.geometry.iloc[fidx]) return dist, fidx results = villages.geometry.apply(nearest_info) villages["dist_m"] = [r[0] for r in results] villages["fac_id"] = [r[1] for r in results] # 按 5 公里一档做覆盖分级，这个阈值按业务场景改 bins = [0, 1000, 3000, 5000, 10000, np.inf] labels = ["1km内", "1-3km", "3-5km", "5-10km", "10km外"] villages["cover_level"] = pd.cut(villages["dist_m"], bins=bins, labels=labels) print(villages["cover_level"].value_counts())sindex.nearest返回的是基于外包框的候选集，严格来说还需要再算一次精确距离，代码里的point.distance做的就是这件事。分档阈值不要照搬，1 公里和 5 公里在不同业务里含义完全不同，末端配送看 3 公里，农资服务看 10 公里，按实际服务能力定。
5.2 网格聚合：把散点变成可统计的面
村点是散点，做热力图或者按区域统计时，需要先聚合成网格。两种做法：规则渔网和六边形网格。渔网实现简单，六边形在邻接关系上更均匀，不会出现方形网格那种对角距离不一致的问题。
import h3 # 用 H3 做六边形聚合，分辨率 7 约 5 平方公里一格 # 参数顺序是 lat, lng，写反了整个网格会跑到地球另一侧 villages["h3_r7"] = [ h3.geo_to_h3(lat, lng, 7) for lat, lng in zip(villages.geometry.y, villages.geometry.x) ] # 按网格统计村点数量和城乡属性构成 agg = villages.groupby("h3_r7").agg( village_cnt=("adcode", "count"), rural_cnt=("urban_rural", lambda s: (s.str.startswith("2")).sum()) ).reset_index() # 把网格中心点还原成几何，方便落图 agg["center"] = agg["h3_r7"].apply(lambda c: h3.h3_to_geo(c)) agg["lat"] = agg["center"].apply(lambda p: p[0]) agg["lon"] = agg["center"].apply(lambda p: p[1]) agg_gdf = gpd.GeoDataFrame( agg, geometry=gpd.points_from_xy(agg["lon"], agg["lat"]), crs="EPSG:4326" ).to_crs(albers)分辨率 7 对应的平均面积在 5 平方公里左右，这个粒度做县域分析刚好，再细会碎，再粗会糊。geo_to_h3的参数顺序是纬度在前、经度在后，这一点和points_from_xy正好相反，写代码时特别容易搞混，建议在函数封装里就把顺序固定死。
5.3 与路网叠加：直线距离不等于实际可达
5.1 算出来的是直线距离。做初步筛选够用，但真要评估可达性，直线距离会系统性低估——山区一个村直线距离 2 公里，实际绕山路可能 8 公里。
常见做法是接网络分析。手上有路网数据的话，用最短路径算法从设施点出发做一次多源 Dijkstra，得到每个村点的路网距离。没有路网数据，可以用等时圈近似：以设施点为中心，按不同出行方式的速度做缓冲区，再和村点做空间关系判断。这种方法精度差一档，但胜在不需要路网，适合做早期摸底。
另一个容易忽略的点是村委会坐标和村域中心不是一回事。村点通常落在村委会所在地，一个行政村下辖几个自然村，村委会在最东边的自然村，那这个点的位置对整个村来说就是偏的。做村级覆盖统计时，这个偏差会累积。解决办法是如果手上有自然村级别的点位，按行政村聚合取质心；没有的话，就在结论里注明「以村委会所在地代表全村」，别让下游误以为是村域质心。
5.4 出图与交付：符号化比配色重要
最后一步是出图。QGIS 里加载 GeoPackage，按urban_rural的首位字符做分类，城镇类和乡村类用两种颜色区分，比按村名随机着色有意义得多。点太密的时候打开聚合渲染，缩放到县级再展开成单点。
交付时我一般给两个文件：一个带完整属性的 GeoPackage 给做分析的人，一个简化到只剩村名、代码、经纬度的 CSV 给做表的人。两份文件同源同版本，避免两边对不上账又回来问。
6. 交付前的三道核验关，以及一个降低返工率的习惯
村点数据这种东西，错得不会很明显，但会错得很均匀。坐标系偏了，所有点一起偏；代码丢零了，所有关联一起错。这类错误不会在单条记录上暴露，只会在汇总结果上体现成一个说不通的数字。所以我每次处理完，都要走三道核验关再往下游交。
第一道，坐标系抽检。随机抽 30 个点，导出成 CSV，把经纬度贴到带影像的地图上逐个看落点。这一步纯人工，慢，但能抓住所有系统性的坐标问题。
# 固定随机种子，保证每次抽的是同一批，方便复现和对比 sample = villages.sample(30, random_state=42)[ ["village_name", "adcode", "lon", "lat"] ] sample.to_csv("spot_check.csv", index=False, encoding="utf-8-sig")第二道，代码对账。把adcode前 6 位截出来去重，和边界数据的县级代码集合做差集，两边都应该接近空。差集规模超过几十个，说明版本错位，得回去查对照表。
# 两边县级代码集合做对称差，看差在哪些代码上 village_codes = set(villages["adcode"].str[:6]) boundary_codes = set(boundary["adcode"].astype(str).str[:6]) print("村点有边界无:", len(village_codes - boundary_codes)) print("边界有村点无:", len(boundary_codes - village_codes))第三道，分布合理性。按县级汇总村点数量，画个直方图，再按省份排个序。正常的分布应该是长尾的，少数县几百个村，多数县几十到一百多个。如果某个县算出几千个村，或者某个县只有个位数，那大概率是代码截取或者空间连接出了问题。
三道关走完，剩下的就是使用层面的经验了。从那以后我每次拿到村点数据，不管对方说得多确定，都强制走一遍这三道关再动分析——血泪教训换来的一条：数据源的口头承诺，永远不如自己抽三十个点看一眼。另外提醒一句，村点数据本身是静态快照，行政区划却在动，用之前先确认手上这份的年份，别拿旧点配新界，也别拿新点去对旧年鉴。
希望帮到你。
本文还有配套的精品资源，点击获取
