Raster image seems to be shifted using leaflet for R

Question

I want to plot some spatial data using the leaflet package in R, however the generated raster image seems to be shifted compared to a reference grid. I suspect map projection issue, but I am not expert on the topic, so any help would be appreciated.

Here is a minimal code to plot the map:

library(leaflet)
library(sp)
library(raster)

set.seed(111)

# create dummy data -rectangular grid with random values
m = 10
n = 10
x = seq(45,48,length.out = m)
y = seq(15,18,length.out = n)
X = matrix(rep(x, each = n), nrow = n)
Y = matrix(rep(y, m), nrow = n)

# collector dataframe
points = data.frame(value = rnorm(n*m), lng = c(Y), lat = c(X))

## create raster grid
s = SpatialPixelsDataFrame(points[,c('lng', 'lat')], data = points)
# set WGS84 projection
crs(s) = sp::CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")

r = raster(s)

# add coloring
pal = colorNumeric(c("#0C2C84", "#f7f7f7", "#F98009"), points$value,
                    na.color = "transparent")

## plot map
leaflet() %>% addProviderTiles("CartoDB.Positron") %>%
   addRasterImage(r, colors = pal, opacity = 0.6)

This produces this map, which is ok at first sight:

If a grid is added to this map:

## grid
dx = diff(x)[1]/2
dy = diff(y)[1]/2

rect_lng = sapply(points$lng, function(t) c(t-dx, t+dx))
rect_lat = sapply(points$lat, function(t) c(t-dy, t+dy))

leaflet() %>% addProviderTiles("CartoDB.Positron") %>% 
  addRectangles(
    lng1=rect_lng[1,], lat1=rect_lat[1,],
    lng2=rect_lng[2,], lat2=rect_lat[2,],
    fillColor = "transparent",
    weight = 1
  ) %>%
  addRasterImage(r, colors = pal, opacity = 0.6)

The map looks like this:

Here we can see that the grids are not matching.

What is the reason of this mismatch? How could it be eliminated? I tried various projections in vain. The only thing that worked was to use addRectangle instead of addRasterImage, however that requires much more computation and slows down the process, thus I want to avoid. Note that in the above example addRectangle is used only for having a reference, in the final code I do not want to use it.

For maps with more cells(grids) the mismatch is quite large, can be larger than the size of a single cell.

Thanks in advance for any help.

EDIT

The problem might be related to the projection issue between ellipsoid and sphere projections, see the last question here:

to convert between WGS84 and mercator on the sphere there will be substantial shifts in the Y mercator coordinates. This is because internally cs2cs is having to adjust the lat/long coordinates from being on the sphere to being on the WGS84 datum which has a quite differently shaped ellipsoid.

However, I was not able to solve the problem with the recommended 'trick': +nadgrids=@null.

Uncomfortable

Answer 1

Author of the leaflet R package here. It looks to me like the raster layer renderer I wrote starts to drift when the source raster has very few pixels relative to the number of pixels that are rendered onscreen. You can see this by making the following modification to the raster:

r1 <- r
nrow(r1) <- 600
ncol(r1) <- 600
r <- resample(r, r1, method = "ngb")

I'll see if I can improve matters on the rendering side, but in the meantime a resample like this may be the easiest workaround, though admittedly it's inelegant.

Archipelago

I think the problem is that the projection operation I do (leaflet::projectRasterForLeaflet) creates an output raster with the same pixel/cell dimensions as the input raster (in this case, 10x10), whereas normally you want to project using the output dimensions you are actually rendering to. But in this case there aren't really fixed output dimensions, since the user can zoom in and out. I'm making a small projected raster and then scaling it up linearly to match the zoom, and that's likely what's introducing the drift. – Archipelago

Thank you for the workaround and for the possible explanation! The resampled raster r is behaving strangely in colorNumeric; I changed the domain to values(r) but I receive a warning that some values are outside the color scale and indeed some pixels are transparent (holes in the raster plot). For the above example domain = c(1.01*min(values(r)), 1.01*max(values(r))) solved the problem, but a more elegant solution would be better. If you have any suggestion in this matter that would be great. Thanks! – Uncomfortable

Hmmm. Does domain = c(minValue(r), maxValue(r)) not work (where r is the already-resampled raster)? Why is 1.01 needed? – Archipelago

No, domain = c(minValue(r), maxValue(r)) does not work, that is why I applied the 1.01 multiplier (tried with the resampled r). – Uncomfortable

@JoeCheng Has there been a solution for this in the package? I have the same issue, and I tried the resampling you suggested. It improved things slightly but the raster and grids are still slightly off. – Nansen