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
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
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
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

TimSalabim TimSalabim · Accepted Answer · 2019-12-11T18:10:06

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.

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

2 Answers