Leaflet中EPSG:27700坐标系下无插值栅格与多边形叠加问题

我需要在Leaflet（或mapview）中显示多个图层，其中包含一个EPSG:27700坐标系的栅格。目前仅能通过默认经纬度投影实现图层叠加，但这会导致栅格重投影并产生插值，而项目要求禁止插值，因此必须在EPSG:27700坐标系下操作。尝试使用CRS.Simple在笛卡尔平面显示，但无法让同坐标系的多边形（或sp对象）与未插值的栅格正确叠加，以下是最小可复现代码：

library("raster")
library("leaflet")
library("eurostat")
library("sf")

## 获取投影到英国网格EPSG27700的UK空间数据框
europe <- get_eurostat_geospatial(resolution = 10, nuts_level = 1,  year = 2021)
UK_spdf <- as_Spatial(europe[grepl("UK", europe$id),])
UK_spdf <- spTransform(UK_spdf, crs("+init=epsg:27700 +units=km +datum=WGS84"))

## 创建一个投影为EPSG:27700的虚拟栅格
r <- rasterize(UK_spdf, raster(UK_spdf, ncols = 100, nrows = 200))

## 在默认绘图中两个图层叠加正常
plot(r) ; plot(UK_spdf, add=TRUE)

## 栅格可以正常加载
leaflet() %>% 
  addRasterImage(r, project = FALSE) ## project=FALSE用于避免插值

## 但多边形无法正确显示？
leaflet() %>% 
  addPolygons(data = UK_spdf)
## 无法正常工作...

## 必须转换成经纬度才能显示：
leaflet() %>% 
  addTiles() %>%
  addPolygons(data = spTransform(UK_spdf, crs("+proj=longlat"))) %>%
  addRasterImage(r)
## 但我们不想这么做，因为这会导致栅格被重投影并产生插值


## 那么如何在简单平面坐标系下让两者叠加？
crs <- leafletCRS(crsClass = "L.CRS.Simple") ## 或许简单投影能解决？
leaflet(options = leafletOptions(crs = crs)) %>% 
  addPolygons(data = UK_spdf) %>%
  addRasterImage(r, project = FALSE)
## 还是无法正常工作...

核心问题 #

Leaflet默认依赖WGS84经纬度投影，直接使用非Web墨卡托/经纬度的投影（如EPSG:27700）时，需要手动对齐图层的坐标原点和缩放逻辑——L.CRS.Simple不会自动适配你的投影坐标系范围。

步骤1：对齐多边形与栅格的坐标系统 #

当用project=FALSE加载栅格时，Leaflet会把栅格左上角作为原点(0,0)，像素直接对应地图坐标。因此需要将多边形坐标转换为以栅格左上角为原点的相对坐标，同时修正y轴方向差异（Leaflet的y轴向下，R栅格的y轴向上）：

# 获取栅格的范围
r_extent <- extent(r)
# 复制多边形数据避免修改原对象
UK_spdf_aligned <- UK_spdf
# 计算x轴偏移：减去栅格左上角x值
UK_spdf_aligned@coords[,1] <- UK_spdf@coords[,1] - r_extent@xmin
# 计算y轴偏移：栅格右上角y值减去多边形y值，反转y轴方向
UK_spdf_aligned@coords[,2] <- r_extent@ymax - UK_spdf@coords[,2]

步骤2：配置Leaflet视图与CRS #

设置Simple CRS，并将地图初始视图适配到栅格尺寸：

# 创建Simple CRS
crs <- leafletCRS(crsClass = "L.CRS.Simple")

# 计算栅格的宽度和高度
r_width <- r_extent@xmax - r_extent@xmin
r_height <- r_extent@ymax - r_extent@ymin

# 加载并对齐图层
leaflet(options = leafletOptions(crs = crs, minZoom = -5, maxZoom = 10)) %>%
  addRasterImage(r, project = FALSE) %>%
  addPolygons(data = UK_spdf_aligned, color = "#ff0000", weight = 2) %>%
  setView(lng = r_width/2, lat = r_height/2, zoom = 0) # 居中显示

关键说明 #

• 调整minZoom和maxZoom参数可以适配大尺寸栅格范围，避免缩放异常。
• 该方法完全避免了栅格重投影插值，所有图层都基于原始EPSG:27700坐标对齐显示。

可选：用mapview简化操作 #

如果使用mapview，可以直接指定CRS并禁用重投影，无需手动调整坐标：

library(mapview)
mapviewOptions(crs = 27700)
mapview(r, project = FALSE) + mapview(UK_spdf)

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