企业官网建设流程全解析

简介：面向铁路地理信息分析、地图可视化与空间查询学习者，提供全国火车站站点位置的GIS矢量数据，可用于绘制站点分布、按区域检索车站、测算站间距离及规划线路等场景。压缩包共8个文件，约405KB，包含矢量主文件及配套索引、投影参考、属性表、字符编码和元数据等类型，分别承担站点坐标存储、快速定位、坐标系统定义、车站名称与所属线路等属性查询、跨语言读取支持，并兼顾投影与编码说明，便于在不同软件环境中稳定读取和复用。已有396人学习，适合需要轻量级全国站点底图的中级用户。数据可用于空间分析、专题制图与业务系统集成，也能作为教学演示和二次开发的基础数据，帮助读者快速搭建车站分布查询与线路规划原型，减少自行采集与整理坐标的成本。

1. 全国火车站站点位置 GIS 数据，本质是一张带几何列和属性列的站点表

拿到「全国火车站点位 GIS」这份资源，很多人的第一反应是双击打开看地图好不好看。但真正决定它能不能用的，是每个点要素上挂了几列属性、坐标落在哪个参考系、有没有一条记录能唯一锁定一个车站实体。GIS 里的全国火车站站点位置数据，剥开渲染层就是一张带几何列的表：站名、一对经纬度、若干归属和等级字段，一行一个站点。

它擅长回答「位置」这一类问题——物流网点到最近站点的距离、区域站点密度分布、换乘可达性粗算、给已有业务表补一列空间坐标做关联。适合做交通与区域分析、物流选址、城市规划出图的人，也适合想拿一份全国尺度真实点位练空间分析管线的从业者。它不擅长的事也很明确：线路走向、时刻表、客流，这份资源给不了，别指望拿它做排班。

真正拉开使用体验差距的，是坐标系和字段清洗这两步。这份笔记按「先定参考系、再洗数据、再做分析、最后避坑」的顺序拆，每一步都能直接抄。

2. 先对齐坐标系再画图：火车站点位的 CGCS2000、WGS84 与投影三选一

站点位置类的数据，最常见的翻车不是算法写错，而是坐标系没对齐。我见过一个很典型的场景：A 同学把点位图层直接叠到在线底图上，整层往一个方向偏了几百米，查了一下午以为是数据本身不准，最后发现问题出在读进来时 CRS 是空的，他随手 set_crs 成了一个不相干的投影。地图上看不出来，因为点还在城市大致范围内，一做距离计算就全错。

所以拿到数据的第一件事永远不是画图，是把参考系确认清楚并记录下来。火车站点位的坐标系无外乎三种情况：地理坐标（CGCS2000 或 WGS84，单位是度）、某个投影坐标（单位是米）、以及少数来源经过非线性偏移处理过的坐标。前两种可以互相转换，第三种只能近似纠正，这是后面避坑章节要重点讲的。

2.1 shapefile 是五件套，缺一件就可能读出错误数据

很多人拿到的是一个压缩包，解压后看到一堆同名不同后缀的文件，以为只需要 .shp。实际上 shapefile 是一个多文件格式，缺哪个都会出问题，而且报错方式各不相同，不看元数据很难定位。

这五个文件里，.prj 和 .cpg 最容易被忽略，也最容易在跨机器复制时丢失。我一般会把整个目录原样打包，不单独挑 .shp 传，血泪经验就来自某次只传了 shp 和 dbf，对方打开后中文全是问号。

2.2 三步判定参考系：读 .prj、看数值量级、量一段距离

判定的顺序不能乱。先读 .prj，这是最权威的来源；如果 .prj 缺失，再用数值量级推断；最后一定要用一个已知量做交叉验证，因为前两步都可能被骗。

第一步，直接看 WKT 原文。GEOGCS 开头说明是地理坐标系，PROJCS 开头说明是投影坐标系。用 ogrinfo 比肉眼读 WKT 快得多：

# -so 只输出摘要不遍历要素，-al 对所有图层生效，大文件下秒出结果 ogrinfo -so -al stations_all.shp # 如果 .prj 存在，直接把 WKT 打出来看前 600 字符 cat stations_all.prj | head -c 600

第二步，用 Python 看数值量级。这一步的判断标准很硬：经度在 -180 到 180、纬度在 -90 到 90 之间，才有可能是经纬度；如果 x 动辄几十万上百万，基本就是米制投影。

import geopandas as gpd # encoding 必须跟 .cpg 声明一致，声明 UTF-8 就用 utf-8，声明 GBK 就用 gbk gdf = gpd.read_file("stations_all.shp", encoding="utf-8") print("CRS:", gdf.crs) print("EPSG:", gdf.crs.to_epsg() if gdf.crs else None) print("要素数:", len(gdf)) print("几何类型:", gdf.geom_type.unique().tolist()) print("x 范围:", round(gdf.geometry.x.min(), 3), round(gdf.geometry.x.max(), 3)) print("y 范围:", round(gdf.geometry.y.min(), 3), round(gdf.geometry.y.max(), 3))

第三步，交叉验证。这一步最容易被跳过，但它能救你。做法是挑一对你大致知道直线距离的枢纽站，用测地线算一遍，再和常识对比。如果算出来差了十倍或差了一个数量级，说明坐标单位或参考系判断错了。

from pyproj import Geod # Geod 用椭球参数算测地线距离，输入必须是经纬度（度），不能喂投影坐标 geod = Geod(ellps="WGS84") a = gdf.iloc[0].geometry b = gdf.iloc[1].geometry _, _, dist_m = geod.inv(a.x, a.y, b.x, b.y) print(f"测地线距离: {dist_m / 1000:.1f} km")

参数说明：ellps 指定椭球，WGS84 和 CGCS2000 的椭球参数差异在厘米级，站点坐标本身的采集精度通常是几十米量级，用哪个都不影响结论。真正要警惕的是「以为它是经纬度，其实是投影坐标」，那种情况下 geod.inv 会返回荒谬的结果，反而成了一个免费的错误探测器。

关于 CGCS2000 和 WGS84 要不要转换，我的实际做法是不转。两者在经纬度层面差异小于点位采集精度，落库时把 EPSG 标成 4490 还是 4326 都可以，但必须在数据字典里写清楚标的是哪个，别让下游以为你做过七参数转换。

2.3 重投影：全国统计用等积，局部测距用等距，别碰 Web 墨卡托

一旦确认原始数据是经纬度，下一步就是选择做分析时用什么投影。这里的选择直接决定你的面积和距离算得对不对，而错误往往很隐蔽——数字看起来正常，只是偏了百分之十几。

我的选型原则是：全国尺度的密度、面积、覆盖统计，用等积投影；局部几百公里内的距离量算，用等距或局部 UTM 分带；任何情况下都不要用 Web 墨卡托（EPSG:3857）做面积统计，它在高纬度的面积拉伸可以到数倍，而点位分布在中高纬时这个误差会被直接带进结论。

import geopandas as gpd gdf = gpd.read_file("stations_all.shp", encoding="utf-8") # 只有在确认原始数值是经纬度、但 .prj 丢失时才这么做；否则不要覆盖已有 CRS if gdf.crs is None: gdf = gdf.set_crs(epsg=4326, allow_override=True) # 全国尺度统计：EPSG:6933 是等积圆柱投影，单位米，面积不失真 gdf_m = gdf.to_crs(epsg=6933) # 局部范围（比如单个城市群）：换成对应 UTM 分带，距离精度更好 # gdf_local = gdf.to_crs(epsg=32650) print(gdf_m.crs) print("重投影后要素数:", len(gdf_m))

参数说明：allow_override=True 表示允许覆盖已存在的 CRS 定义，这个参数只在明确知道原始坐标是什么的时候用，否则等于把错误合法化。to_crs 会重新计算每个点的坐标值，要素数不变，几何类型也不变，所以可以用 len 前后对比做一次完整性检查。

重投影做完，一定要回到地图上看一眼。不是看配色，是看海岸线附近的点是不是落在水里、边界附近的点有没有明显甩出去。投影转换本身不会产生这种问题，但如果原始数据里有脏点，投影后会放大偏差，早发现比后面分析完再发现省事得多。

3. 用 GeoPandas 把站点表读成干净表：字段修正、站名归一与脏点剔除

坐标系定好之后，第二层问题是属性。火车站站点数据的属性表看起来简单，实际脏得很有特点：站名有全角半角混用、有括号里塞别名、有同名不同站、有坐标落在水里的记录。这些问题不会让程序报错，只会让结果悄悄偏掉，属于典型的黑匣子型错误。

清洗的顺序我固定为三步：先做结构体检，确认字段类型和空值分布；再做站名归一，因为它是对外关联的唯一键；最后做几何层面的脏点剔除。顺序不能反，站名没归一就做去重，会误删真实站点；几何没体检就做空间连接，会得到笛卡尔积。

3.1 结构体检：字段、空值、几何有效性，一次看完

这一步的目标是把数据摸清楚，不急着改。重点看三样：字段名和类型是否符合预期、关键字段的空值比例、几何是否有效。

import geopandas as gpd import pandas as pd gdf = gpd.read_file("stations_all.shp", encoding="utf-8") print("字段:", gdf.columns.tolist()) print("类型:\n", gdf.dtypes) # 空值按比例排序，先处理影响最大的字段 na_ratio = gdf.isna().mean().sort_values(ascending=False) print("空值比例:\n", na_ratio[na_ratio > 0].round(3)) print("空几何:", int(gdf.geometry.is_empty.sum())) print("无效几何:", int((~gdf.geometry.is_valid).sum()))

字段类型这块最常见的坑是坐标被读成字符串。某些来源的数据把经纬度存在属性表的两列文本里，几何是由工具临时拼出来的，这时候 to_crs 会报错或者静默失败。判断方法很简单：如果 gdf.crs 是 None 但几何数值又像经纬度，就要怀疑是坐标列被当成普通字段了。

# 如果几何是空的，但属性表里有 lon / lat 两列，用它们重建点 from shapely.geometry import Point if gdf.geometry.is_empty.all() and {"lon", "lat"}.issubset(gdf.columns): gdf["geometry"] = [ Point(float(lon), float(lat)) for lon, lat in zip(gdf["lon"], gdf["lat"]) ] gdf = gdf.set_geometry("geometry").set_crs(epsg=4326, allow_override=True)

参数说明：Point 接收的是 (x, y)，也就是先经度后纬度，这个顺序反了是新手期最高频的错误，反了之后点会跑到地球另一侧，肉眼一看就知道，但程序不会报错。set_crs 之后一定要跟一次 to_crs 到米制，否则后面的 buffer 又会按度来算。

3.2 站名归一：全角半角、括号别名、后缀「站」字三板斧

站名是这份数据里最有价值的字段，因为它是唯一能跟外部业务表（比如物流网点表、客流表）关联的键。但它也是最脏的：同一个站在不同来源里可能写成「某某站」「某某」「某某（东）」「某某(东)」四种形态。不归一就直接 merge，你会发现匹配率低得离谱。

我的归一函数只做四件事，不做过度清洗，因为改得越多越容易出现不可逆的误伤：

import re import pandas as pd def norm_station_name(s): if pd.isna(s): return None s = str(s).strip() s = s.replace("（", "(").replace("）", ")") # 全角括号转半角 s = s.replace("\u3000", "") # 去掉全角空格 s = re.sub(r"\s+", "", s) # 去掉所有空白字符 s = re.sub(r"\(.*?\)", "", s) # 去掉括号内的别名或方位说明 if s.endswith("站") and len(s) > 2: s = s[:-1] # 去掉尾部「站」字，但保留短名 return s gdf["name_norm"] = gdf["name"].map(norm_station_name) print("归一后唯一站名数:", gdf["name_norm"].nunique()) print("原始唯一站名数:", gdf["name"].nunique())

参数说明：len(s) > 2 这个判断是为了避免把两个字的站名砍成单字，导致大面积重名。re.sub(r"(.*?)", "", s) 里的非贪婪匹配很关键，用贪婪匹配会把两个括号之间的内容全吃掉。如果你的数据里括号内是正式的方位信息（比如「某某(东广场)」实际是两个不同实体），这一步就不要做，改成保留括号内容单独成一列。

归一之后立刻做同名检测。同名不等于重复，全国范围内同名不同地的站点是真实存在的，所以判断依据必须是坐标距离，不能是名字。

gdf_m = gdf.to_crs(epsg=6933) # 到米制才能算真实距离 rows = [] for name, sub in gdf_m.groupby("name_norm"): if len(sub) < 2: continue coords = list(zip(sub.geometry.x, sub.geometry.y)) max_d = max( ((x1 - x2) ** 2 + (y1 - y2) ** 2) ** 0.5 for i, (x1, y1) in enumerate(coords) for x2, y2 in coords[i + 1:] ) rows.append((name, len(sub), round(max_d))) dup = pd.DataFrame(rows, columns=["name_norm", "count", "max_dist_m"]) dup = dup.sort_values("max_dist_m", ascending=False) print(dup.head(20))

判定阈值我一般这样定：组内最大距离小于 200 米，视为同一实体重复录入，保留属性更全的一条；大于 5000 米，视为同名不同站，两条都留，但要在名字后面加区分后缀；介于两者之间的，人工看，这批通常不多，几十条以内。

3.3 脏点剔除：用稳健统计，别用三倍标准差

站点坐标出错的方式很有特点：不是随机噪声，而是整条记录用了行政区中心坐标、或者录入时小数点位置错了。这类点会严重拉偏后续的核密度和最近邻统计，但在地图缩略图上完全看不出来。

我用的方法是基于中位数和 MAD 的稳健离群检测，而不是三倍标准差。原因很实际：全国站点分布本身就不均匀，用均值和标准差做阈值，会把稀疏区域的合法站点大批误杀。

import numpy as np x = gdf.geometry.x.to_numpy() y = gdf.geometry.y.to_numpy() def mad_filter(v, k=6): med = np.median(v) mad = np.median(np.abs(v - med)) # 1.4826 是让 MAD 与标准差可比的缩放系数 return np.abs(v - med) <= k * 1.4826 * mad mask = mad_filter(x) & mad_filter(y) print("剔除点数:", int((~mask).sum()), "/", len(gdf)) gdf_clean = gdf[mask].copy().reset_index(drop=True)

参数说明：k 取 6 是宽松阈值，k 越小剔除越狠。不要用 k=3，那个尺度下稀疏区域的合法站点会被成片误杀。剔除完必须抽样看几条被剔除的记录，确认它们确实是脏点而不是边远地区的正常站点，这一步不能省，因为全国尺度的点分布本身就有长尾。

4. 火车站点位能做的三件事：核密度、最近邻与缓冲区抽样

数据洗干净之后，能做的分析其实就三大类，绝大多数业务需求都能归到其中一类：站点在哪儿扎堆（核密度）、某个位置离最近站点多远（最近邻）、站点覆盖了多大范围（缓冲区）。这三个方法本身都不复杂，难的是参数怎么给、结果怎么解释、边界在哪儿。

下面每个方法我都会给出可运行的代码和一个「什么时候不该用」的判断。因为在实际项目里，选错方法比写错代码代价更大。

4.1 核密度：看站点扎堆，但带宽决定结论

核密度的直觉很朴素：把每个点看成一个小山包，叠起来就是密度面。问题全在带宽上。带宽给大了，全国会糊成一整片红色，看不出任何结构；给小了，每个点都是一个孤立的尖峰，等于没做平滑。

import numpy as np from scipy.stats import gaussian_kde gdf_m = gdf_clean.to_crs(epsg=6933) xy = np.vstack([gdf_m.geometry.x, gdf_m.geometry.y]) # bw_method 是相对带宽，按样本标准差的比例给，实测 0.05~0.15 是比较合理的区间 kde = gaussian_kde(xy, bw_method=0.08) # 在规则的格网中心上求密度值，而不是在全图逐像素算，速度快一个量级 gx = np.linspace(xy[0].min(), xy[0].max(), 300) gy = np.linspace(xy[1].min(), xy[1].max(), 300) gxx, gyy = np.meshgrid(gx, gy) positions = np.vstack([gxx.ravel(), gyy.ravel()]) z = kde(positions).reshape(gxx.shape) print("密度最大值:", z.max(), " 网格形状:", z.shape)

参数说明：bw_method 传浮点数时表示相对带宽，传字符串 "scott" 或 "silverman" 是按样本量自动算。实测在几千个点的规模下，自动带宽往往偏大，全国图上会糊成一片，手动调到 0.05 到 0.15 之间更好看。300x300 的格网是精度和速度的平衡点，要出图再提到 800 也行。

这里有个容易被忽略的点：gaussian_kde 的复杂度跟样本量成正比，样本上万之后会明显变慢。如果拿到的数据规模超过几万条，改成先做格网计数再上高斯滤波，速度快得多，结果在视觉上几乎一致。

4.2 最近邻：算每个业务点到最近站点的距离

这是业务侧用得最多的一个分析。核心就一句：目标点集到站点集的最近距离。别用双重循环，几万个点两两算会卡死，用 KD 树。

from scipy.spatial import cKDTree import numpy as np # 站点集和业务点集必须在同一个投影下，否则距离单位不一致 gdf_m = gdf_clean.to_crs(epsg=6933) biz_m = biz.to_crs(epsg=6933) station_xy = np.vstack([gdf_m.geometry.x, gdf_m.geometry.y]).T biz_xy = np.vstack([biz_m.geometry.x, biz_m.geometry.y]).T tree = cKDTree(station_xy) # k=1 只取最近的一个；workers=-1 开全部核心并行 dist, idx = tree.query(biz_xy, k=1, workers=-1) biz_m["nearest_station_m"] = dist biz_m["nearest_station"] = gdf_m["name_norm"].to_numpy()[idx] print(biz_m[["nearest_station", "nearest_station_m"]].head())

参数说明：query 返回的两个数组长度跟查询点数量一致，dist 单位是米（因为已经投影过），idx 是在站点数组里的下标，用它回填站名。workers 参数在 scipy 1.9 之后才有，老版本删掉即可。

一个常见的误用是拿站点集自己查自己。那种情况下要传 k=2，因为第一个结果永远是点自己，距离为 0。很多教程贴的代码在这里没写清楚，直接 query(xy, k=1) 会得到一堆 0，还以为是数据有问题。

4.3 缓冲区：算覆盖范围，注意分辨率和并集的开销

缓冲区最常见的用法是算「站点周边 N 公里覆盖了多少面积」。看起来简单，两个参数决定结果质量和运行时间：半径单位、圆的离散精度。

# resolution=16 时，一个圆被切成 64 段，视觉上已经够圆 # 半径单位是米，因为 gdf_m 已经投影过 buf = gdf_m.geometry.buffer(5000, resolution=16) # 求并集，得到去重后的总覆盖范围；geopandas 0.13 之前叫 unary_union cover = buf.union_all() print("覆盖面积 km2:", round(cover.area / 1e6, 1))

参数说明：resolution 越大越圆，但顶点数线性增长，几千个圆做并集时能明显感觉到差别。全国尺度下我一般用 8 到 16，局部出图才用 32。半径直接给数字就行，单位跟你当前 CRS 一致，投影到米制后就是米——这也是为什么前面反复强调要先转投影，如果在 4326 下写 buffer(0.05)，那是 0.05 度，不是 5 公里。

并集操作在点数多时是性能瓶颈。常见的优化是先用 sjoin 做粗筛，把互相重叠的圆分组成小簇，每个簇内部单独做并集，最后再合并。数据量在几千个点以内其实没必要优化，直接算就行。

5. 排查与避坑：坐标偏移、字段乱码、面积算错这四类翻车

上面几步做完，正常的分析流程就能跑通了。但实际项目里真正耗时间的从来不是主流程，是那些不报错、只给你一个错误结果的环节。下面四条是我在火车站点位数据上反复遇到的，每条都按「现象、原因、解决」写清楚，遇到问题时可以对着查。

5.1 现象：整层点位在底图上单向偏移几百米

原因基本可以锁定在坐标来源上。一部分在线地图服务产出的坐标经过非线性偏移处理，跟 CGCS2000 或 WGS84 不是同一个参考系，直接叠加就会整体错位。这种偏移不是简单的平移，不同位置的偏移量还不一样，所以用平移量去凑是凑不准的。

解决分两步。先确认是不是偏移坐标：拿两三个位于城市中心的站点，跟公开的经纬度对照，如果偏移方向和量级在城市间不一致，就是非线性偏移。确认之后，要么找一份未偏移的原始数据，要么用公开发布的开源纠偏算法做近似还原。我的习惯是优先换数据源，纠偏只能做到米级近似，做物流选址这类需要精确距离的场景不够用。如果只是出示意图，偏移几百米其实可以接受，别为这个耗一整天。

5.2 现象：中文站名读出来是乱码或问号

原因在 .dbf 的编码声明。shapefile 本身没有强制编码标准，OGR 在读不到 .cpg 时会按 ISO-8859-1 猜，中文自然全乱。这时候数据结构是好的，只是字符被错误解码。

解决方式是先看 .cpg 内容，再用对应的 encoding 参数读：

# 看编码声明，输出通常是 UTF-8、936 或 GBK cat stations_all.cpg

# 声明是 UTF-8 就用 utf-8，是老数据常见 GBK 用 gbk 或 gb18030 gdf = gpd.read_file("stations_all.shp", encoding="utf-8") # 如果已经读错了，可以转成 GeoJSON 时指定编码再重读 # ogr2ogr -f GeoJSON -lco ENCODING=UTF-8 stations.geojson stations_all.shp

用 gb18030 比 gbk 更保险，它是超集，能覆盖更多生僻字。如果两种都试过还是乱码，那就不是编码问题，是 .dbf 本身就写坏了，只能换数据。

5.3 现象：缓冲区和面积算出来的数字离谱

原因几乎总是同一个：在经纬度坐标系下直接调 buffer 或 area。度不是长度单位，1 度在不同纬度对应的实际距离差好几倍，算出来的「面积」没有任何物理意义。

解决就一句：任何涉及距离和面积的操作之前，先 to_crs 到米制投影。验证方法也很简单，算完看一眼量级——单个站点的 5 公里缓冲区面积应该在 78 平方公里左右，如果算出来是 0.0002 或者 7 亿，那就是单位搞错了。

5.4 现象：空间连接后结果条数暴增几十倍

原因是几何无效。重复点、自相交、坐标完全相同的多个记录，在 sjoin 时会产生笛卡尔积，行数直接炸开。这个错误最坑的地方是它不报错，你只会觉得结果「怎么这么多」。

解决方式是在所有空间连接之前加一道体检：

# 先修复无效几何，再按几何去重 gdf["geometry"] = gdf.geometry.make_valid() gdf = gdf.drop_duplicates(subset="geometry").reset_index(drop=True) print("去重后:", len(gdf))

make_valid 在 shapely 2.x 里是几何方法，老版本用 buffer(0) 代替也能达到类似效果。去重时用几何本身做键，不要用站名，同名不同站是真实存在的，按名字去重会误删。

5.5 一份可以一次性跑完的自检清单

上面四条都是事后排查，我更推荐事前自检。把下面这段存成脚本，每次拿到新数据先跑一遍，能挡掉八成问题：

def sanity_check(gdf, name="stations"): report = { "要素数": len(gdf), "CRS已定义": gdf.crs is not None, "无空几何": bool((~gdf.geometry.is_empty).all()), "几何均有效": bool(gdf.geometry.is_valid.all()), "几何全为点": bool((gdf.geom_type == "Point").all()), "站名无空值": bool(gdf["name"].notna().all()), "坐标量级正常": bool( abs(gdf.geometry.x).max() < 180 and abs(gdf.geometry.y).max() < 90 ), } print(f"=== {name} ===") for k, v in report.items(): print(f"{k}: {v}") return report sanity_check(gdf_clean)

最后一项「坐标量级正常」是关键，它把「经纬度还是投影坐标」这个最容易搞错的判断固化成了自动检查。如果它返回 False 而 CRS 又显示已定义，八成是 CRS 被错误地 set 过，这时候别继续往下做，先把参考系捋清楚。

6. 把它接进真实分析：用站点做一次网格化可达性粗算

站点数据最终多数会落到「可达性」这类问题上：某片区域到最近火车站有多远、多快。这里有个容易被跳过的判断——用直线距离还是路网距离。直线距离快、稳、不需要额外数据，但会系统性低估山区和跨水系区域的时间成本；路网距离准，但要依赖路网数据和路径规划，成本高一个量级。

我的做法是先跑直线距离做一版粗算，把量级和分布摸清楚，再决定值不值得上路网。粗算的流程是：把研究区切成规则格网，算每个格网中心到最近站点的距离，按阈值分级。

import numpy as np import pandas as pd from scipy.spatial import cKDTree gdf_m = gdf_clean.to_crs(epsg=6933) station_xy = np.vstack([gdf_m.geometry.x, gdf_m.geometry.y]).T tree = cKDTree(station_xy) # 格网边长 20 km，可以在全国尺度上把计算量控制在百万级以内 step = 20000 minx, miny, maxx, maxy = gdf_m.total_bounds xs = np.arange(minx, maxx, step) ys = np.arange(miny, maxy, step) gx, gy = np.meshgrid(xs, ys) centers = np.column_stack([gx.ravel(), gy.ravel()]) dist, idx = tree.query(centers, k=1, workers=-1) grid = pd.DataFrame({ "x": centers[:, 0], "y": centers[:, 1], "nearest_m": dist, "nearest_station": gdf_m["name_norm"].to_numpy()[idx], }) # 按 10/30/60 km 分三级，阈值按业务场景调，不要照搬 bins = [0, 10000, 30000, 60000, np.inf] labels = ["10km内", "10-30km", "30-60km", "60km外"] grid["level"] = pd.cut(grid["nearest_m"], bins=bins, labels=labels) print(grid["level"].value_counts(normalize=True).round(3))

参数说明：step 决定格网密度，20 公里在全国尺度上大约得到几十万个格网单元，KD 树查询几秒内完成；如果只做一个省，step 可以降到 2 到 5 公里。bins 的三个阈值是我按客运场景随手定的，做货运或做应急响应时阈值完全不同，一定要按场景重设，不要拿别人的阈值直接用。

这版结果能回答的问题是「哪些格网连 60 公里内都没有站点」，以及「站点密度高的区域跟人口或业务分布是否重合」。它回答不了「开车要多久」，因为直线距离没有考虑地形和路网，山区误差可能到两三倍。所以我一般会把这个结果当成一张筛子，先圈出可疑区域，再对那几百个格网单独跑路网路径规划，而不是全量做。

分级的解释也要留个心眼：60 公里外不等于没有服务，可能只是站点布局本身稀疏，这在人口稀少区域是正常现象，不构成问题。反过来，10 公里内有站也不等于方便，中间可能隔着一条没有桥的河。这些边界写进结论里，比把图做得好看更重要。

从那次把一个没定义 CRS 的图层直接拿去做缓冲区、结果面积差了六个数量级之后，我养成了一个习惯：任何站点数据进流程之前，先跑一遍第 5 章那段自检，CRS、几何有效性、坐标量级三项必须全绿才往下走。这份全国火车站点位数据本身质量不错，字段结构清楚、点位密度也够，真正决定产出质量的，是你在这三道工序上花的功夫——参考系对齐、站名归一、脏点剔除，每一步都不复杂，但跳过任何一步都会在后面某个环节以你想不到的方式找回来。希望帮到你。

本文还有配套的精品资源，点击获取