Raster image seems to be shifted using leaflet for

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.

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.