简介这份资源提供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。三种形态各有取舍选哪种取决于你后面要干什么形态优点局限适合的场景ShapefileGIS 软件通用性最好多文件、字段名限 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, encodinggbk) # 打印完整坐标系定义 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, encodinggbk, # 乱码时优先试 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}, # 不加这句前导零会被吃掉 encodingutf-8-sig # 带 BOM 的 CSV 用这个避免第一列名带不可见字符 ) # 从经纬度列构建几何并声明源坐标系 gdf gpd.GeoDataFrame( df, geometrygpd.points_from_xy(df[lon], df[lat]), crsEPSG: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 projaea lat_125 lat_247 lat_00 lon_0105 x_00 y_00 datumWGS84 unitsm 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, encodinggbk).to_crs(albers) # 村点也投影过去两边坐标系必须一致才能做空间关系判断 pts gdf.to_crs(albers) # 左连接保留所有村点未匹配的用 NaN 标出 joined gpd.sjoin( pts, boundary[[adcode, geometry]].rename(columns{adcode: adcode_bnd}), howleft, predicatewithin # within 表示点落在面内边界上的点用 intersects ) # 统计未匹配数量这批点要么在边界外要么边界数据有空洞 unmatched joined[joined[adcode_bnd].isna()] print(未匹配数量:, len(unmatched), 占比:, len(unmatched) / len(joined))predicate参数有两个常用值within严格判断点是否在面内部落在边界线上的点会被判为不匹配intersects则包含边界线。全国尺度的数据用intersects更稳妥未匹配率能低一些。howleft保证村点一条不丢这样你才能看到未匹配的那部分长什么样。未匹配率超过 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, layervillage, driverGPKG, spatial_indexTrue )spatial_indexTrue这个参数在小数据上看不出差别数据量上到十万级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.10105E11或者前面的零没了变成 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, indexFalse, encodingutf-8-sig)程序读取的那份则始终保证dtype{adcode: str}。两套文件分开互不干扰。4.3 现象按村名关联一条业务数据匹配出几十条现象是把业务表和村点按village_name做 join结果行数从几万膨胀到几十万一条记录对应十几个村。原因是村名重名率远超直觉。同一个村名在不同省份、不同县、甚至同一个县的不同乡镇都会重复出现。用村名当键等于放弃唯一性约束。解决方式是建立复合键。首选adcode12 位全码唯一性最好。如果业务表只有村名没有代码退而求其次用「县级代码前 6 位 乡镇名 村名」三元组匹配率会降到可接受范围。实在什么代码都没有就只能靠空间位置反查把业务地址地理编码成点再和村点做最近邻匹配取最近的那个村。这条路误差最大只适合做粗略归属不能用于统计口径。4.4 现象属性表里中文村名全是问号现象是打开属性表village_name一列全是???其他数字字段正常。原因是.dbf的编码和读取时指定的编码不一致。.cpg文件如果写的是UTF-8而实际内容是 GBK或者.cpg压根不存在读取端就会猜错。解决有两条路。第一条在读取时显式指定encodinggbk大多数国内分发的 Shapefile 都是 GBK。第二条用 GDAL 的命令行工具直接转一次编码转完再读# 把 Shapefile 转成 GeoPackage同时强制按 UTF-8 写入属性表 ogr2ogr -f GPKG village_utf8.gpkg village_2021.shp \ -lco ENCODINGUTF-8 \ --config SHAPE_ENCODING GBK--config SHAPE_ENCODING GBK告诉 GDAL 源文件按 GBK 解析-lco ENCODINGUTF-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], binsbins, labelslabels) 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, geometrygpd.points_from_xy(agg[lon], agg[lat]), crsEPSG: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_state42)[ [village_name, adcode, lon, lat] ] sample.to_csv(spot_check.csv, indexFalse, encodingutf-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))第三道分布合理性。按县级汇总村点数量画个直方图再按省份排个序。正常的分布应该是长尾的少数县几百个村多数县几十到一百多个。如果某个县算出几千个村或者某个县只有个位数那大概率是代码截取或者空间连接出了问题。三道关走完剩下的就是使用层面的经验了。从那以后我每次拿到村点数据不管对方说得多确定都强制走一遍这三道关再动分析——血泪教训换来的一条数据源的口头承诺永远不如自己抽三十个点看一眼。另外提醒一句村点数据本身是静态快照行政区划却在动用之前先确认手上这份的年份别拿旧点配新界也别拿新点去对旧年鉴。希望帮到你。本文还有配套的精品资源点击获取 
长河君 · 2026/10/26 10:33:20 · 阅读 1.2k 2021全国村名点数据：坐标系纠偏与空间连接实战
简介这份资源提供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。三种形态各有取舍选哪种取决于你后面要干什么形态优点局限适合的场景ShapefileGIS 软件通用性最好多文件、字段名限 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, encodinggbk) # 打印完整坐标系定义 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, encodinggbk, # 乱码时优先试 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}, # 不加这句前导零会被吃掉 encodingutf-8-sig # 带 BOM 的 CSV 用这个避免第一列名带不可见字符 ) # 从经纬度列构建几何并声明源坐标系 gdf gpd.GeoDataFrame( df, geometrygpd.points_from_xy(df[lon], df[lat]), crsEPSG: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 projaea lat_125 lat_247 lat_00 lon_0105 x_00 y_00 datumWGS84 unitsm 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, encodinggbk).to_crs(albers) # 村点也投影过去两边坐标系必须一致才能做空间关系判断 pts gdf.to_crs(albers) # 左连接保留所有村点未匹配的用 NaN 标出 joined gpd.sjoin( pts, boundary[[adcode, geometry]].rename(columns{adcode: adcode_bnd}), howleft, predicatewithin # within 表示点落在面内边界上的点用 intersects ) # 统计未匹配数量这批点要么在边界外要么边界数据有空洞 unmatched joined[joined[adcode_bnd].isna()] print(未匹配数量:, len(unmatched), 占比:, len(unmatched) / len(joined))predicate参数有两个常用值within严格判断点是否在面内部落在边界线上的点会被判为不匹配intersects则包含边界线。全国尺度的数据用intersects更稳妥未匹配率能低一些。howleft保证村点一条不丢这样你才能看到未匹配的那部分长什么样。未匹配率超过 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, layervillage, driverGPKG, spatial_indexTrue )spatial_indexTrue这个参数在小数据上看不出差别数据量上到十万级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.10105E11或者前面的零没了变成 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, indexFalse, encodingutf-8-sig)程序读取的那份则始终保证dtype{adcode: str}。两套文件分开互不干扰。4.3 现象按村名关联一条业务数据匹配出几十条现象是把业务表和村点按village_name做 join结果行数从几万膨胀到几十万一条记录对应十几个村。原因是村名重名率远超直觉。同一个村名在不同省份、不同县、甚至同一个县的不同乡镇都会重复出现。用村名当键等于放弃唯一性约束。解决方式是建立复合键。首选adcode12 位全码唯一性最好。如果业务表只有村名没有代码退而求其次用「县级代码前 6 位 乡镇名 村名」三元组匹配率会降到可接受范围。实在什么代码都没有就只能靠空间位置反查把业务地址地理编码成点再和村点做最近邻匹配取最近的那个村。这条路误差最大只适合做粗略归属不能用于统计口径。4.4 现象属性表里中文村名全是问号现象是打开属性表village_name一列全是???其他数字字段正常。原因是.dbf的编码和读取时指定的编码不一致。.cpg文件如果写的是UTF-8而实际内容是 GBK或者.cpg压根不存在读取端就会猜错。解决有两条路。第一条在读取时显式指定encodinggbk大多数国内分发的 Shapefile 都是 GBK。第二条用 GDAL 的命令行工具直接转一次编码转完再读# 把 Shapefile 转成 GeoPackage同时强制按 UTF-8 写入属性表 ogr2ogr -f GPKG village_utf8.gpkg village_2021.shp \ -lco ENCODINGUTF-8 \ --config SHAPE_ENCODING GBK--config SHAPE_ENCODING GBK告诉 GDAL 源文件按 GBK 解析-lco ENCODINGUTF-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], binsbins, labelslabels) 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, geometrygpd.points_from_xy(agg[lon], agg[lat]), crsEPSG: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_state42)[ [village_name, adcode, lon, lat] ] sample.to_csv(spot_check.csv, indexFalse, encodingutf-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))第三道分布合理性。按县级汇总村点数量画个直方图再按省份排个序。正常的分布应该是长尾的少数县几百个村多数县几十到一百多个。如果某个县算出几千个村或者某个县只有个位数那大概率是代码截取或者空间连接出了问题。三道关走完剩下的就是使用层面的经验了。从那以后我每次拿到村点数据不管对方说得多确定都强制走一遍这三道关再动分析——血泪教训换来的一条数据源的口头承诺永远不如自己抽三十个点看一眼。另外提醒一句村点数据本身是静态快照行政区划却在动用之前先确认手上这份的年份别拿旧点配新界也别拿新点去对旧年鉴。希望帮到你。本文还有配套的精品资源点击获取