Raster over international date line

[tim-salabim]

This issue was filed in mapview (r-spatial/mapview#6)

As mapview uses leaflet and this is a leaflet issue I reproduce it here.

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))
leaflet() %>% addTiles() %>% addRasterImage(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))
leaflet() %>% addTiles() %>% addRasterImage(r)

Any idea how I can make this work?

Rotating the raster helps somewhat, but it would be nice for the raster to not be split into two parts at either end of the map (as is the case for polygons, see below).

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(rotate(r))

There is no problem with polygons (ideally this is the behaviour we would want for ratsers too):

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"))
leaflet() %>% addTiles() %>% addPolygons(data = p)

Thanks in advance,
Tim

[quantifish]

Hi all,

I came across this same issue as well and wanted to add that this only seems to be an issue with rasters. Adding circles or markers with points beyond 180 works fine as well:

p <- data.frame(lng = runif(100, 165, 185), lat = runif(100, -50, -31))
leaflet() %>% addTiles() %>% addMarkers(lat = p$lat, lng = p$lng)

While a raster can be rotated, it is far from ideal and looks terrible when the raster is split into two parts. Are there any ideas floating around on how to fix this up? I'd even take a quick fix for now.

Regards,
Darcy

[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
xr1 = extent(b1) %>% as.vector() %>% .[1:2] %>% diff()
xr2 = extent(b2) %>% as.vector() %>% .[1:2] %>% diff()
bb = c(extent(b1), extent(b2))[[which.max(c(xr1,xr2))]]

# map with bisected chunks, worldCopyJump = TRUE, fitBounds()
leaflet(options=leafletOptions(worldCopyJump = TRUE)) %>% 
  addTiles() %>%
  addRasterImage(b1, opacity=.5) %>%
  addRasterImage(b2, opacity=.5) %>%
  addGraticule() %>%
  fitBounds(bb@xmin, bb@ymin, bb@xmax, bb@ymax)

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.