Raster over the 180th degree meridian

[johnForne]

Not sure if this is an issue per se, or more of a question / request...
As part of New Zealand's Environmental Reporting Programme, we regularly publish maps that span the 180th degree meridian (the anti-meridian). We're big fans of the sf package and able to display vector data in leaflet() maps. But I have an issue that I don't seem to be able to seamlessly display rasters that span the 180 line.

I've found that @tim-salabim (rstudio/leaflet#225) raised an issue with rasters over the 180 line that still isn't closed. I've also found that @jmlondon (r-spatial/sf#280) raised an issue relating to sf::st_transform not honoring +lon_wrap #280 (which showed me how to seamlessly display vector data over the 180 line - though the option of splitting the raster each side of the 180 line created a line down the middle). So I'm raising this issue in the stars space to check whether this issue can be addressed from the stars side of things...

I can create sf (vector) and stars objects with CRS of 4326 and that also have longitude values in 0-360 range (cf -180 to 180). Both display seamlessly in ggplot. As does the vector in Leaflet. But for some reason leaflet doesn't seem to know how to handle stars objects spanning the 180 line. Why can leaflet handle vector data spanning the 180 line - but not rasters? Can anyone please explain what the issue is and, better, offer a solution?

Thanks in advance,

John

The following code provides an example of a raster that I've trying to display seamlessly...

library(leaflet)
library(tidyverse)
library(sf)
library(stars)

r <- raster(xmn = 165, xmx = 195, ymn = -50, ymx = -31, crs = "+init=epsg:4326", nrows = 50, ncols = 50)
r[] <- rnorm(ncell(r))
leaflet() %>% addTiles() %>% addRasterImage(r)
s <- st_as_stars(r)
s <- s %>%
  st_wrap_dateline(c("WRAPDATELINE=YES", "DATELINEOFFSET=60"))
leaflet() %>% addTiles() %>% addStarsImage(s)

[mdsumner]

There's never any easy way out of this problem, the issue in this instance is that projecting to Mercator (for leaflet) truncatest at longitude = 180. You can create a copy and shift that westwards to the other side of the anti-meridian and it "works":

r <- raster(xmn = 165, xmx = 195, ymn = -50, ymx = -31, crs = "+proj=longlat", nrows = 50, ncols = 50)
r[] <- rnorm(ncell(r))
r2 <- raster::shift(r, dx = -360)
leaflet() %>% addTiles() %>% addRasterImage(r) %>% addRasterImage(r2)

The hassle is this kind of technique doesn't work for every case, so there's no good general answer. You could alternatively cut the raster in two parts, and so you only have as many cells as you started with. The only way to have a single raster is to have it extend from -180, 180 with all the empty cells in between. There's no tool in stars that can help built-in afaik, but the same copy and shift tech could work.

sf::st_wrap_dateline is for vector not raster

[johnForne]

Kia ora Michael, much appreciated :)
Following your suggest, I noticed that shifting the raster seems to distort it. If I run the code you provide it makes the raster have "lines" across it. However, I think that the raster should have more square pixels of randomly distributed colours. Try `plot(r2). Afaics rasters have this distorted look when they touch the 180 line... don't understand why this strange stretching occurs though. So I think the best we can do is crop the raster into two parts on either side of the 180 line.

Thanks heaps,
John
`

[johnForne]

p.s.
I've had a further dig into the issue of rasters over the 180 line and found that somehow when we used to use ArcGIS online we managed to publish rasters over the 180 line.

From what I can tell the ESRI World Ocean Base map is in EPSG 3857. And the webservice that we pull from ArcGIS online has a spatial reference of 4326.

The .mxd that we used to publish the webservice had a CRS of 4326. Yet the feature class for Sea Surface Temperature within the .mxd was a custom CRS (a kind of Lambert Conformal Conic). Though I think the fact that the feature class was not 4326 and was probably transformed on the fly is isn't significant...

The point is that I was intrigued what ArcGIS online is doing so that it is able to recognise that -180W actually joins next to 180E??? Is this type of thing that we can get added to leaflet in R?

Would really appreciate it if anyone is able to look into this, thanks.

[edzer]

@tim-salabim any ideas about this?

[tim-salabim]

@edzer I had another look at the (oldest open) mapview issues and might have a solution, though, admittedly not very elegant.

r-spatial/mapview#6 (comment)

Maybe stars can somehow avoid the issues that arise from snaping in raster?