Mapview in R and problem with projections/crs in Pacific area

Question

I am trying to plot data in R using mapView for the grid in the Pacific that crosses the longitude 180deg. The native crs is "+proj=merc +lon_0=150 +lat_ts=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0". The example of equivalent data is:

test.coords<-as.data.frame(cbind(c(runif(15,-180,-130),runif(5,160,180)),runif(20,40,60)))
test.sp <- SpatialPointsDataFrame(coords = cbind(test.coords$V1,test.coords$V2), data = test.coords,
                                 proj4string = CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"))
test.sp <- spTransform(test.sp, 
          "+proj=merc +lon_0=150 +lat_ts=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0")
plot(test.sp)
mapview(test.sp)

When using plot function, the points are centered, while using mapView they get split between two sides of the map (linked image). Can the view using mapView be centered on a different longitude?

Using native crs for this map does not help. I get an error.

mapview(plot, native.crs=TRUE)

Warning message: sf layer is not long-lat data

Thank you

mapView image Pacific

Can you please provide a reproducible example? We might need some checks in mapview to reproject correctly. What do you mean by 'not displaying correctly' when using native.crs? — TimSalabim, Dec 07 '19 at 10:22
Hi, I added additional sample data explaining the challenge. The issue is that I am working with longitudes 160 to 180 and -180 to -130. — Barbara Hutniczak, Dec 09 '19 at 17:12
See my answer below for a potential solution. Commenting just to clarify that `mapview(tes.sp, native.crs = TRUE)` (after reprojecting) does work as expected. You will not get any background map when you specify `native.crs = TRUE` and your data is not in `longlat`. The error you mention is a warning and can be ignored in this case. — TimSalabim, Dec 11 '19 at 18:21
What do you mean by "When using plot function, the points are centered"? At every step your points are split in the west and east hemispheres — mdsumner, Dec 11 '19 at 19:35

score 1 · Answer 1 · answered Dec 11 '19 at 18:10

Ok, so this is a semi-optimal solution to your question for several reasons:

this will only work for point data
it may not scale well for large point collections
projecting to another CRS may not work anymore (though that is not really related, as the issue at hand is about visualising)

For sets of lines or polygons, we would need another approach, but I am currently not aware of a sp/sf native solution to round-tripping coordinates of an object.

library(sp)
library(sf)
#> Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3
library(mapview)

test.coords<-as.data.frame(cbind(c(runif(15,-180,-130),runif(5,160,180)),runif(20,40,60)))
test.sp <- SpatialPointsDataFrame(coords = cbind(test.coords$V1,test.coords$V2), data = test.coords,
                                  proj4string = CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"))

test_sf = st_as_sf(test.sp)

shift = function(x) {
    geom = st_geometry(x)
    st_geometry(x) = st_sfc(
        lapply(seq_along(geom), function(i) {
            geom[[i]][1] = ifelse(geom[[i]][1] < 0, geom[[i]][1] + 360, geom[[i]][1])
            return(geom[[i]])
        })
        , crs = st_crs(geom)
    )
    return(x)
}

mapview(shift(test_sf)) +
    mapview(test_sf, col.regions = "orange")

^{Created on 2019-12-11 by the reprex package (v0.3.0)}

I chose to use sf rather than sp for this because mapview uses sf internally and was hoping to find a general approach to this problem which could potentially be integrated.

Hi Tim, thanks for looking at the problem at hand. This will indeed work fine for pints, but I also have associated polygons that will not align with the points. — Barbara Hutniczak, Dec 12 '19 at 20:05
So the challenge to build a mapview with various layers centered on Pacific remains. — Barbara Hutniczak, Dec 12 '19 at 20:11
I've just submitted a [pull request](https://github.com/r-spatial/sf/pull/12189) to **sf** which should take care of this also for polygons and lines too. As soon as it is merged, I will add a new argument to **mapview** to allow for recentering/shifting features to be centered on Pacific. — TimSalabim, Dec 12 '19 at 21:10

score 1 · Answer 2 · answered Mar 13 '20 at 23:42

Here is a workaround for polygons - it produces some warnings, but overall does the job.

# transform grid to standard latitude and longitude
# original grid in crs "+proj=merc +lon_0=150 +lat_ts=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"
grid<-spTransform(grid,"+init=epsg:4326")
summary(coordinates(grid))

#       V1               V2       
# Min.   :-179.8   Min.   :37.17  
# 1st Qu.:-168.8   1st Qu.:51.70  
# Median :-154.3   Median :55.02  
# Mean   :-118.0   Mean   :54.65  
# 3rd Qu.:-131.4   3rd Qu.:58.31  
# Max.   : 179.8   Max.   :65.56  

out <- lapply(1:length(grid), function(i) grid[i, ])

for (i in 1:length(grid)) {

  cds <- slot(slot(slot(out[[i]], "polygons")[[1]], "Polygons")[[1]], "coords")

  cds_polygon <- Polygons(list(Polygon(cds)), ID="out.x")
  out.x <- SpatialPolygons(list(cds_polygon))
  proj4string(out.x) <- CRS("+init=epsg:4326")

  if (cds[1,1]>0) { # depending to which side one want to move polygons that are on both sides of International Date Line
    out.y <- elide(out.x, shift=c(-360,0))
  } else {
    out.y<-out.x}

  if(i==1){grid.shifted<-out.y} else {
    grid.shifted<-raster::union(grid.shifted,out.y)
  }
}

# add data in the original grid
grid.shifted <- SpatialPolygonsDataFrame(grid.shifted, grid@data, match.ID = F)

mapview(grid.shifted)

Mapview in R and problem with projections/crs in Pacific area

2 Answers2