国家范围20km×20km网格构建、实体匹配及置换测试实现方案咨询

前置依赖 #

你可以基于Python生态完成全流程落地，需要用到的核心库如下：

• 地理数据处理：geopandas、shapely
• 数值计算：numpy、scipy
• 数据持久化：pickle、pandas

1. 20km等距网格生成（基于Haversine距离逻辑） #

无需依赖可视化工具，直接生成可调用的结构化网格数据：

1. 读取国家边界文件，提取全域 bounding box 极值：

import geopandas as gpd
import numpy as np

country_border = gpd.read_file("你的边界文件路径.shp/geojson")
min_lon, min_lat, max_lon, max_lat = country_border.total_bounds

2. 计算固定纬度步长：纬度每度对应地表距离约为111km，因此20km对应的纬度差为：
   lat_step = 20 / 111 ≈ 0.1802°
3. 生成网格结构：

grid_dict = {}
grid_id = 0
# 遍历所有纬度带
current_lat = min_lat
while current_lat < max_lat:
    # 计算当前纬度下的经度步长：经度1度距离=111km*cos(当前纬度弧度值)
    lon_step = 20 / (111 * np.cos(np.radians(current_lat)))
    current_lon = min_lon
    while current_lon < max_lon:
        # 生成网格多边形，过滤掉不在国家边界内的网格
        grid_poly = shapely.geometry.box(current_lon, current_lat, current_lon+lon_step, current_lat+lat_step)
        if grid_poly.intersects(country_border.geometry.iloc[0]):
            grid_dict[f"grid_{grid_id}"] = {
                "min_lon": current_lon,
                "max_lon": current_lon + lon_step,
                "min_lat": current_lat,
                "max_lat": current_lat + lat_step,
                "center": (current_lon + lon_step/2, current_lat + lat_step/2),
                "poly": grid_poly
            }
            grid_id += 1
        current_lon += lon_step
    current_lat += lat_step

生成的grid_dict就是可直接调用的结构化数据，可通过pickle.dump(grid_dict, open("grid_20km.pkl", "wb"))持久化存储，后续直接加载使用。

2. 关联实体-网格映射构建 #

1. 先批量匹配所有实体到所属网格，优先用空间索引提升匹配效率：

import pandas as pd
# 假设你的实体数据存在df里，字段为entity_id, lon, lat
entity_df = pd.read_csv("你的实体数据路径.csv")
# 转成geopandas格式
entity_gdf = gpd.GeoDataFrame(
    entity_df,
    geometry=gpd.points_from_xy(entity_df.lon, entity_df.lat),
    crs=country_border.crs
)
# 构建网格GeoDataFrame用于空间匹配
grid_gdf = gpd.GeoDataFrame(
    [{"grid_id": k, "geometry": v["poly"]} for k,v in grid_dict.items()],
    crs=country_border.crs
)
# 空间连接匹配，1次操作完成所有实体的网格归属判断
entity_to_grid = gpd.sjoin(entity_gdf, grid_gdf, predicate="within").set_index("entity_id")["grid_id"].to_dict()

2. 生成成对关联实体的映射：

# 假设你的关联对数据存在pair_df里，字段为x1_id, x2_id
pair_to_grids = {}
for _, row in pair_df.iterrows():
    pair_to_grids[(row["x1_id"], row["x2_id"])] = [
        entity_to_grid[row["x1_id"]],
        entity_to_grid[row["x2_id"]]
    ]

最终输出的pair_to_grids完全符合你要求的(X1,X2): [grid20,grid55]格式。

3. 置换检验（Permutation Test）实现 #

操作逻辑如下：

1. 先计算原始观测值：比如你要验证关联实体的空间距离更近，就先计算所有原始关联对的网格中心点Haversine距离的均值，记为obs_value
2. 执行置换操作：

import random
perm_times = 1000 # 置换次数可自行调整
perm_values = []
all_grid_ids = list(grid_dict.keys())
for _ in range(perm_times):
    # 随机打乱所有实体的网格归属，保持关联对的配对关系不变
    random.shuffle(all_grid_ids)
    perm_entity_grid = dict(zip(entity_to_grid.keys(), all_grid_ids[:len(entity_to_grid)]))
    # 计算本次置换后的特征值
    perm_dist = 0
    for x1, x2 in pair_to_grids.keys():
        g1, g2 = perm_entity_grid[x1], perm_entity_grid[x2]
        # 这里可以替换成你要验证的任意特征计算逻辑
        perm_dist += haversine(grid_dict[g1]["center"], grid_dict[g2]["center"])
    perm_values.append(perm_dist / len(pair_to_grids))

3. 显著性判断：统计obs_value在perm_values分布中的位置，得到p值，若p<0.05则可认为关联特征存在非随机的地理规律。

性能优化提示 #

• 实体数量超过10万时，优先用geopandas的sindex空间索引做匹配，比逐点判断快100倍以上
• 网格生成后可直接导出为pkl文件，无需每次重复计算
• 置换操作可通过多进程并行加速，10000次置换一般几分钟就能跑完

内容的提问来源于stack exchange，提问作者skynaive