Rasters over international date line

[yvanrichard]

Mapview seems to not be supporting rasters going over the international date line.

This works fine:

library(raster)
r <- raster(xmn = 165, xmx = 180, ymn = -50, ymx = -31, crs = "+init=epsg:4326", nrows = 50, ncols = 50)
r[] <- rnorm(ncell(r))
mapview(r)

But not this, with the raster being cropped and the cells being stretched:

r <- raster(xmn = 165, xmx = 185, ymn = -50, ymx = -31, crs = "+init=epsg:4326", nrows = 50, ncols = 50)
r[] <- rnorm(ncell(r))
mapview(r)

Any idea how I can make this work? Thank you for your package. This is quite promising.

---

R version 3.2.3 (2015-12-10)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 14.04.3 LTS

[yvanrichard]

I forgot to mention that there is no problem with polygons crossing the international date line. For example, this works:

p <- SpatialPolygons(list(Polygons(list(Polygon(cbind(x=c(165, 165, 185, 185, 165),
                                                     y=c(-50, -31, -31, -50, -50)))), 1)),
                    proj4string=CRS("+init=epsg:4326"))
mapview(p)

[tim-salabim]

So far, I can only confirm that this is an issue of the underlying leaflet package (though it may well be a problem of leaflet.js).

r <- raster(xmn = 165, xmx = 185, ymn = -50, ymx = -31, crs = "+init=epsg:4326", nrows = 50, ncols = 50)
r[] <- rnorm(ncell(r))
leaflet() %>% addTiles() %>% addRasterImage(r)

The only way of getting around this in a half decent manner (with a warning) is to rotate the raster, however, this will draw the raster bit beyond the date line at the far western end of the map instead of simply extending it towards the east.

mapview(rotate(r))
leaflet() %>% addTiles() %>% addRasterImage(rotate(r))

All other attempts of cutting the raster at the date line and drawing it in 2 instances have failed so far.

I have reproduced your issue at https://github.com/rstudio/leaflet/issues/225. Let's see what they come up with.

[bbest]

I had a similar problem, and found using rotate() to cause the raster to be prohibitively huge since it fills out the global extent. Instead I bisected the raster into chunks using shift() and crop().

Here's the reproducible code and results...

define raster

library(raster)
library(leaflet)

# define raster that exceeds [-180,180]
r <- raster(xmn = 165, xmx = 185, ymn = -50, ymx = -31, crs = "+init=epsg:4326", nrows = 50, ncols = 50)
r[] <- rnorm(ncell(r))

map with > 180 chopped off

leaflet() %>% addTiles() %>% addRasterImage(r) %>% addGraticule()

option A: rotate()

# option 1: rotate raster
a = rotate(a)

# map rotated raster
leaflet() %>% addTiles() %>% addRasterImage(a) %>% addGraticule()

option B: shift() & crop() into 2 rasters

# option 2: shift and bisect into raster pieces left and right of dateline to minimize size
b  <- shift(r, -360)
b1 <- crop(b, extent(-180,180,-90,90), snap="in") %>% trim()
b2 <- crop(b, extent(-360,-180,-90,90), snap="in") %>% shift(360) %>% trim()

# map with bisected chunks, worldCopyJump = TRUE
leaflet(options=leafletOptions(worldCopyJump = TRUE)) %>% 
  addTiles() %>%
  addRasterImage(b1) %>%
  addRasterImage(b2) %>%
  addGraticule()

inspect difference in rasters between options A & B

The resulting maps of the two options above look the same, however you'll notice that the dimensions of the rotated raster a are more than 18x larger than b1 and b2!

# inspect:    (xmin,  xmax, ymin, ymax);           (nrow, ncol, ncell)
r  # extent:   165  ,  180  ,  -50, -31 ; dimensions : 50,  50,  2500
a  # extent : -179.8,  180.2,  -50, -31 ; dimensions : 50, 900, 45000
b  # extent : -195  , -180  ,  -50, -31 ; dimensions : 50,  50,  2500    
b1 # extent : -179.8, -175  ,  -50, -31 ; dimensions : 50,  12,   600
b2 # extent :  165  ,  179.8,  -50, -31 ; dimensions : 50,  37,  1850

For option B, I also used trim(), which doesn't do anything for this perfectly rectangular example, but if the non-NA latitudinal range of the 2 split rasters differed, than the extent (and memory size) would also shrink.

map based on extent of one of 2 rasters from option B

We can then map based on the extent of the larger raster from option B.

# get bounding box on bigger raster
bb = extent(b2) %>% as.matrix()

# map with bisected chunks, worldCopyJump = TRUE, fitBounds()
leaflet(options=leafletOptions(worldCopyJump = TRUE)) %>% 
  addTiles() %>%
  addRasterImage(b1, opacity=.5) %>%
  addRasterImage(b2, opacity=.5) %>%
  addGraticule() %>%
  fitBounds(bb["x","min"], bb["y","min"], bb["x","max"], bb["y","max"])

Although initial display cuts off at dateline of 180, we only need to pan over to the other side enough for the worldCopyJump = TRUE setting to take effect. The addGraticule() helps visualize this.

[tim-salabim]

Thanks!

[tim-salabim]

Sparked by https://github.com/r-spatial/stars/issues/256 I think I finally found a way to make this work. The approach is only slightly modifies option B by @bbest above. The issue with the original suggestion was that croping the raster into two parts of either side of the dateline sometimes resulted in a gap between the two rasters, depending on what snap ended up doing. We need to set snap = "in" to ensure we don't cross the dateline. But, this means we sometimes end up with extent values somewhere between 179 and 180 (e.g. 179.36 - see how this can change by changing xmx in the raster definition). Hence, I've added another shift call to slightly move the two sides of the raster towards the dateline.

library(raster)
library(leaflet)

# define raster that exceeds [-180,180]
r <- raster(xmn = 165, xmx = 211, ymn = -50, ymx = -31, crs = "+init=epsg:4326", nrows = 50, ncols = 50)
r[] <- rnorm(ncell(r))

leaflet() %>% addTiles() %>% addRasterImage(r) %>% addGraticule()

# option 2: shift and bisect into raster pieces left and right of dateline to minimize size
b  <- shift(r, -360)
b1 <- crop(b, extent(-180,180,-90,90), snap="in") %>% trim()
b1 = shift(b1, -180 - xmin(b1))
b2 <- crop(b, extent(-360,-180,-90,90), snap="in") %>% shift(360) %>% trim()
b2 = shift(b2, 180 - xmax(b2))

# map with bisected chunks, worldCopyJump = TRUE
leaflet(options=leafletOptions(worldCopyJump = TRUE)) %>% 
    addTiles() %>%
    addRasterImage(b1) %>%
    addRasterImage(b2) %>%
    addGraticule()

One thing that needs to happen with this approach though is that the color scales need to be defined globally so that coloring respects the value range of the entire image...

Maybe there is a better approach? At least this seems to work

[johnForne]

Kia ora Tim, much appreciated. I was super excited as you got rid of the gap over the 180 line - but, unfortunately, then released that you solution introduced a potentially more serious error in the data (not just the smoothness of how it shows on the map)...

The issue is that in shifting the cropped rasters it results in the coordinates of the pixels being shifted so that they may not long reflect the what was originally observed on the ground at a particular lat/long and recorded in a raster cell for that area.

So unless there's any other trick for getting leaflet to draw pixels that actually span the 180 line? - I'm a bit stumped.

Thanks again