使用Cartopy add_wms绘制EPSG2056(LV95)瑞士地形图偏移170米求助

问题背景 #

使用Matplotlib+Cartopy调用瑞士联邦地形测量局（Swisstopo）的WMS服务时，基于LV95（EPSG:2056）投影绘制的地图，与已知LV95坐标的对象、官方下载的LV95格式数据叠加时出现约170米偏移。尝试通过wms_kwargs={'srs':'EPSG:2056'}显式指定坐标系时，触发重复参数错误。

问题原因 #

1. WMS版本参数差异：Swisstopo的WMS服务采用WMS 1.3.0标准，该版本使用CRS参数指定坐标系；而旧版WMS 1.1.1使用SRS参数。Cartopy的ax.add_wms自动处理参数时，可能因版本适配问题导致坐标系请求不准确。
2. 二次转换精度损失：ax.add_wms默认会先将LV95范围转换为WGS84（EPSG:4326）请求瓦片，再转换回LV95，这个过程会引入局部高精度坐标系的转换误差。

解决方案 #

方案1：使用cartopy.io.img_tiles.WMS直接请求LV95瓦片

绕过自动转换流程，直接指定WMS 1.3.0版本和LV95坐标系，确保请求的瓦片原生匹配目标投影，彻底避免转换误差。

修改后的完整代码：

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from cartopy.io.img_tiles import WMS

LV95 = ccrs.epsg(2056)
x0 = 2682956; y0 = 1246234; extent = 500; 
xmin=x0-extent/2; xmax=x0+1.5*extent; ymin=y0-extent; ymax=y0+extent

# 初始化Swisstopo WMS服务，指定版本和坐标系
swisstopo_wms = WMS(
    url='https://wms.geo.admin.ch/',
    layers='ch.swisstopo.swisssurface3d-reliefschattierung-multidirektional',
    wms_version='1.3.0',
    crs=LV95
)

fig = plt.figure(figsize=(8,8), dpi=96)
ax = fig.subplots(subplot_kw={'projection':LV95})

# 添加WMS瓦片，13为缩放级别（可根据范围大小调整）
ax.add_image(swisstopo_wms, 13)

ax.set_xticks([xmin, x0, xmax]); ax.set_yticks([ymin, y0, ymax]) 
ax.set_xlim(xmin,xmax); ax.set_ylim(ymin,ymax)
ax.ticklabel_format(axis='x', useOffset=x0)
ax.ticklabel_format(axis='y', useOffset=y0)
ax.set_xlabel("(LV95) X [m]"); ax.set_ylabel("(LV95) Y [m]")

plt.show()

方案2：强制add_wms使用正确的CRS参数

若坚持使用ax.add_wms，需传递CRS参数（而非SRS）适配WMS 1.3.0标准，避免重复参数冲突：

# 修改原代码中的add_wms部分
ax.add_wms(
    wms='https://wms.geo.admin.ch/',
    layers='ch.swisstopo.swisssurface3d-reliefschattierung-multidirektional',
    wms_kwargs={'CRS': 'EPSG:2056'}
)

注：若仍出现重复参数报错，优先使用方案1。

关键注意事项 #

• 缩放级别：方案1中的13是Swisstopo WMS的固定缩放级别，范围越小可选择越高的级别（如14、15）以获得更清晰的瓦片。
• 数据一致性：所有叠加的本地数据必须严格使用EPSG:2056投影，避免跨投影转换带来的额外误差。

内容的提问来源于stack exchange，提问作者Jörg Rychen