5

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...
Yollanda Beetroot
  • 323
  • 6
  • 15
  • Could you please tell what your expected output is? – Quinten Mar 19 '23 at 09:26
  • The expected output is a raster which shall have spatial points overlaid over. In other word, display things in something else than lat-long. To be integrated in a shiny app that let's you click on pixels and display some characteristics of the said pixel. Interpolation of the raster creates local artefacts, and I want to be able to work in planar (like plot() allows you to). I'm doing it already with ggplot(), which conveniently does not force background reprojection on your spatial layers, but I wanted to take advantage of the fluidity of leaflet for zooming and exploring the raster. – Yollanda Beetroot Mar 20 '23 at 07:53

1 Answers1

0

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)
I_O
  • 4,983
  • 2
  • 2
  • 15
  • Thanks for your time @I_O. Your solution brings me to the same point I'm stuck at: I have to let leaflet reproject my layers (which causes raster interpolation) so that they overlay perfectly. This is not good. I need to find a way to have leaflet letting me work with my planar layers which natively overlay, as with the function plot() in my example. Leaflet is doing some hidden reprojection work in the background, and I can't find how to deactivate it. Do you see what I mean? – Yollanda Beetroot Mar 06 '23 at 14:26
  • 1
    I see. I tried supplying leaflet with EPSG:27700 versions of raster and vector, set the leaflet CRS to 27700 (via leaflet options) and prevent reprojection from side of Leaflet. Only that the boundaries of raster and vector mismatch by up to 80 km now. – I_O Mar 06 '23 at 18:02
  • Thanks for trying. I think the solution is in working with those layers as non-geographical layer who do not need to be projected. I guess the "L.CRS.Simple" is made for this purpose but I can't manage to have it work either... – Yollanda Beetroot Mar 07 '23 at 13:12
  • 1
    @Yollanda: I tried the "L.CRS.Simple", following this post: https://stackoverflow.com/questions/70096385/projected-coordinate-reference-systems-in-leaflet-and-r . The curious effect was to have the vector shape zoomable along a raster of constant size. I'm lost :-) – I_O Mar 07 '23 at 13:36