raster and polygons in leaflet, without raster interpolation

312 views Asked by At

I'm trying to display several layers in leaflet (or mapview), one of which is a raster in EPSG:27700. The only way I manage to adequately overlay those layers is through the default latlong projection, which implies a reprojection of the raster and therefore its interpolation. I cannot have interpolation in this project, so I need to work at the EPSG:27700.

How can I display additional layers to an unprojected raster? I have tried using CRS.Simple, as I'd like to display every thing in a simple cartesian plan, but without success. I don't mind loosing the beautiful background tiles. But whatever I try, I cannot have my polygon (also EPSG27700) layer (or any sp object) to correctly display with my not interpolated raster. I hope the MWE below illustrates efficiently my problem:

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

## get UKK spdf projected on british grid EPSG27700
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"))

## build a dummy raster projected on EPSG:27700
r <- rasterize(UK_spdf, raster(UK_spdf, ncols = 100, nrows = 200))

## the two layers overlay well in default plots
plot(r) ; plot(UK_spdf, add=TRUE)

## raster can be loaded 
leaflet() %>% 
  addRasterImage(r, project = FALSE) ## project=FALSE to prevent interpolation

## how to get the polygons right?
leaflet() %>% 
  addPolygons(data = UK_spdf)
## does not work...

## you need to have it in lat long:
leaflet() %>% 
  addTiles() %>%
  addPolygons(data = spTransform(UK_spdf, crs("+proj=longlat"))) %>%
  addRasterImage(r)
## but we don't want that, as it implies that our raster will have to be reprojected and therefore interpolated


## so how to have them together on a simple planar coordinate system?
crs <- leafletCRS(crsClass = "L.CRS.Simple") ## maybe simple projection can help?
leaflet(options = leafletOptions(crs = crs)) %>% 
  addPolygons(data = UK_spdf) %>%
  addRasterImage(r, project = FALSE)
## does not work...
1

There are 1 answers

4
I_O On

If you don't mind using other spatial packages, you can try this:

library(eurostat)
library(sf)
library(dplyr)
library(stars) ## provides `st_rasterize`
library(leafem) ## provides `addStarsImage`
library(leaflet)

europe <- get_eurostat_geospatial(resolution = 10, 
                                  nuts_level = 1,
                                  year = 2021)

UK <- europe |> filter(CNTR_CODE == 'UK')

r_27700 <- UK |>
  select(geometry) |>
  st_transform(27700) |> 
  st_rasterize(nx = 100, ny = 200)


crs_option <- leafletCRS(code = "EPSG:27700")


leaflet(options = leafletOptions(crs = crs_option)) |>
  addStarsImage(r_27700, project = TRUE) |> ## let Leaflet handle (re)projection
  addPolygons(data = UK)