Altair叠加边界与点图层未对齐问题求助

问题根源 #

你的点数据使用WGS84地理坐标系（EPSG:4326）（单位为度数，范围大致在±90/±180），而边界shapefile的geometry坐标量级大，说明它采用了投影坐标系（伦敦常用的是OSGB36，EPSG:27700，单位为米，坐标值通常在几万到几十万级别）。两种坐标系的单位和数值范围完全不匹配，直接叠加必然导致图层错位。

解决方案：统一两个图层的CRS #

必须将点数据和边界数据转换为同一坐标系，以下是两种可行方案：

方案1：将边界数据转换为WGS84（EPSG:4326） #

把边界数据转成和点数据一致的地理坐标系，这样可以直接用longitude/latitude编码点数据：

import geopandas as gpd
import altair as alt

# 读取边界Shapefile
shape_data = gpd.read_file("london_boundaries.shp")
# 转换CRS为WGS84
shape_data_wgs84 = shape_data.to_crs(epsg=4326)

# 你的点数据（假设为包含lat、long字段的DataFrame）
point_data = ...

# 绘制叠加地图
base_map = alt.Chart(shape_data_wgs84).mark_geoshape(
    fill='lightgray',
    stroke='white'
)

points_layer = alt.Chart(point_data).mark_circle(
    color='darkred',
    size=40
).encode(
    longitude='long:Q',
    latitude='lat:Q'
)

base_map + points_layer

方案2：将点数据转换为边界数据的投影坐标系 #

如果需要保留边界数据的投影坐标系（比如更适合局部地图的精度），可以将点数据转换为该坐标系，之后用投影坐标的x/y值编码：

import geopandas as gpd
import altair as alt
from shapely.geometry import Point

# 读取边界数据并获取其CRS
shape_data = gpd.read_file("london_boundaries.shp")
target_crs = shape_data.crs

# 处理点数据：转换为GeoDataFrame并转换CRS
point_data = ... # 包含lat、long字段的DataFrame
# 创建点几何对象
geometry = [Point(xy) for xy in zip(point_data['long'], point_data['lat'])]
# 初始CRS设为WGS84
point_gdf = gpd.GeoDataFrame(point_data, geometry=geometry, crs="EPSG:4326")
# 转换到边界数据的CRS
point_gdf_proj = point_gdf.to_crs(target_crs)

# 绘制叠加地图
base_map = alt.Chart(shape_data).mark_geoshape(
    fill='lightgray',
    stroke='white'
)

points_layer = alt.Chart(point_gdf_proj).mark_circle(
    color='darkred',
    size=40
).encode(
    x='geometry.x:Q',  # 投影坐标系的x坐标
    y='geometry.y:Q'   # 投影坐标系的y坐标
)

base_map + points_layer

注意事项 #

• 先确认边界数据的CRS：执行print(shape_data.crs)即可查看，伦敦区域的Shapefile大概率是EPSG:27700
• 优先用Geopandas做坐标转换：比Altair内置的transform_calculate地理转换更稳定可靠，后者仅支持部分坐标系转换，且容易出现精度问题
• 若尝试Altair内置转换（不推荐），可参考以下代码（仅适用于Vega支持的CRS）：

points_layer = alt.Chart(point_data).mark_circle(color='darkred', size=40).transform_calculate(
    # 将WGS84经纬度转换为EPSG:27700投影坐标
    projected='project([datum.long, datum.lat], "EPSG:4326", "EPSG:27700")'
).encode(
    x='projected[0]:Q',
    y='projected[1]:Q'
)

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