Geopandas：如何还原未知平移的西班牙地图以实现图层叠加？

问题分析 #

你之前用整个图层的西南点计算偏移量，这会受伊比利亚半岛等其他区域的边界干扰，导致加那利群岛的平移偏移量计算错误。正确的做法是单独针对加那利群岛的要素计算偏移，而非用整个图层的边界。

具体解决方案 #

1. 精准提取加那利群岛的要素子集 #

首先从两个图层中单独分离出加那利群岛的部分，避免其他区域干扰：

import geopandas as gpd
import shapely.affinity
import fiona

# 读取并转换CRS（保留原有逻辑）
radon = gpd.read_file('./data/potencial_radon.gdb', layer=fiona.listlayers('./data/radon.gdb')[0]).to_crs(epsg=4326)
municipalities = gpd.read_file('./data/municipalities.gdb').to_crs(epsg=4326)

# 提取风险图中的加那利群岛要素（已确认是索引0-2）
radon_canarias = radon.loc[[0, 1, 2]].copy()

# 从市政图中提取加那利群岛的市政区域
# 方法1：通过属性筛选（如果有行政区划名称/代码字段，比如名称包含"Canarias"）
municipalities_canarias = municipalities[municipalities['nombre'].str.contains('Canarias', case=False)].copy()

# 方法2：通过空间范围筛选（加那利群岛在EPSG:4326下的大致范围：西经17.5-19.5，北纬27-29）
municipalities_canarias = municipalities.cx[-19.5:-17.5, 27:29].copy()

2. 计算加那利群岛专属的平移偏移量 #

基于加那利群岛自身的边界计算偏移，而非整个图层的边界：

# 获取市政图加那利区域的边界（返回值格式：minx, miny, maxx, maxy）
mun_canarias_bounds = municipalities_canarias.total_bounds
# 获取风险图加那利区域的边界
radon_canarias_bounds = radon_canarias.total_bounds

# 计算平移偏移量：用市政图的西南点减去风险图的西南点
xoff = mun_canarias_bounds[0] - radon_canarias_bounds[0]
yoff = mun_canarias_bounds[1] - radon_canarias_bounds[1]

3. 执行平移并合并图层 #

将加那利群岛的风险要素平移后，和原风险图的其他要素合并：

# 复制原风险图避免修改原始数据
moved_radon = radon.copy()

# 对加那利群岛的三个要素执行平移
for idx in [0, 1, 2]:
    moved_radon.at[idx, 'geometry'] = shapely.affinity.translate(
        moved_radon.loc[idx]['geometry'],
        xoff=xoff,
        yoff=yoff
    )

# 可视化验证匹配效果
import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(10,10))
municipalities_canarias.plot(ax=ax, edgecolor='black', facecolor='none')
moved_radon.loc[[0,1,2]].plot(ax=ax, alpha=0.5)
plt.show()

4. 更精准的控制点匹配（可选） #

如果上述方法仍有偏差，可以用已知地标点计算偏移：
比如找特内里费岛的端点或首府坐标，在市政图中获取该点的坐标，再在风险图中找到对应的点，直接用两点的坐标差作为偏移量：

# 示例：假设市政图中某地标点坐标为（-16.5785, 28.0932）
# 风险图中对应点的坐标为（x_r, y_r）
xoff = -16.5785 - x_r
yoff = 28.0932 - y_r

验证方法 #

平移完成后，可以用空间连接验证匹配度：

# 检查风险图加那利要素与市政图的空间匹配情况
join_result = gpd.sjoin(moved_radon.loc[[0,1,2]], municipalities_canarias, how='inner', predicate='intersects')
print(join_result)

内容的提问来源于stack exchange，提问作者Miquel Ferriol Galmés