sf::st_transform not honoring +lon_wrap

[jmlondon]

Hello

I'm working with a number of lines that cross back and forth of 180. When working in a projection (e.g. 3571), this isn't an issue. But, if I want to add those lines to a leaflet map, I need to transform to 4326. With sp::spTransform() I've been able to pass along the +lon_wrap=180 specification and this converts the longitude coordinates to 0-360 which leaflet respects

library(sp)
library(sf)
library(magrittr)
library(leaflet)

l <- "LINESTRING(-150.5 57.0,-160.5 57.0,-170.5 57.0,179.5 57.0,169.5 57.0,159.5 57.0)" %>% 
  sf::st_as_sfc(crs = "+init=epsg:4326") 

l_sp <- as(l,"Spatial")
l_sp <- sp::spTransform(l_sp,"+proj=longlat +datum=WGS84 +lon_wrap=180")
sp::proj4string(l_sp)
bbox(l_sp)

m <- l_sp %>% 
  leaflet() %>% 
  addProviderTiles("Esri.OceanBasemap") %>% 
  addPolylines(weight = 4, opacity = 0.35, color = 'black') 
m

However, when I try to accomplish the same thing with sf, the +lon_wrap appears to be ignored.

m <- l %>% 
  sf::st_transform("+proj=longlat +datum=WGS84 +lon_wrap=180") %>% 
  leaflet() %>% 
  addProviderTiles("Esri.OceanBasemap") %>% 
  addPolylines(weight = 4, opacity = 0.35, color = 'black') 
m

l %>% sf::st_transform("+proj=longlat +datum=WGS84 +lon_wrap=180") %>% 
  sf::st_crs()
l %>% sf::st_transform("+proj=longlat +datum=WGS84 +lon_wrap=180") %>% 
  sf::st_bbox()

> sessionInfo()
R version 3.3.3 (2017-03-06)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
Running under: macOS Sierra 10.12.3

leaflet_1.1.0.9000
sp_1.2-4
sf_0.3-4

Is this a known/expected difference between sp::spTransform and sf::st_transform? Is there a different approach I should be taking?

[edzer]

Yes: sp (rgdal) directly calls routines in proj.4, sf calls the interface through GDAL, and the latter seems to simply filter out the +lon_wrap part. Moving things the way sp does it would be a lot of work, and loose access to extras gdal offers, such as better handling of polygons that partly cannot be projected e.g. when they're behind the globe. But gdal has a wrap-dateline option, now in bd8d900 :

> st_wrap_dateline(st_sfc(st_linestring(rbind(c(-179,0),c(179,0))), crs = 4326))
Geometry set for 1 feature 
geometry type:  MULTILINESTRING
dimension:      XY
bbox:           xmin: -180 ymin: 0 xmax: 180 ymax: 0
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
MULTILINESTRING((-179 0, -180 0), (180 0, 179 0))

This option is also in line with the geojson rfc, which prescribes that longitudes should be in [-180,180], but geometries crossing the antimeridian should be cut there (see par. 3.1.9).

[jmlondon]

Thanks for the explanation.

I built the latest sf via devtools and started playing around with st_wrap_dateline. A few things I've noted:

l1 <- st_sfc(st_linestring(rbind(c(-179,0),
                                 c(179,0))), crs = 4326) %>% 
  sf::st_wrap_dateline()

l2 <- st_sfc(st_linestring(rbind(c(-150.5,57.0),
                                c(159.5,57.0))),crs=4326) %>% 
  sf::st_wrap_dateline()

If we examine l1 and l2, the results are different:

> l1
Geometry set for 1 feature 
geometry type:  MULTILINESTRING
dimension:      XY
bbox:           xmin: -180 ymin: 0 xmax: 180 ymax: 0
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
MULTILINESTRING((-179 0, -180 0), (180 0, 179 0))

> l2
Geometry set for 1 feature 
geometry type:  LINESTRING
dimension:      XY
bbox:           xmin: -150.5 ymin: 57 xmax: 159.5 ymax: 57
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
LINESTRING(-150.5 57, 159.5 57)

I also tried to replicate my previous example in leaflet using st_wrap_dateline instead of casting to sp. But, I'm still not getting the expected result. The linestring is split at 180, but leaflet doesn't recenter the extent. I've cc'd @bhaskarvk in case there's something on the leaflet side I should be doing.

m <- l1 %>% 
  leaflet() %>% 
  addProviderTiles("Esri.OceanBasemap") %>% 
  addPolylines(weight = 4, opacity = 0.35, color = 'black') 
m

[edzer]

There is another parameter, set by

l2 %>% st_wrap_dateline(c("WRAPDATELINE=YES", "DATELINEOFFSET=60"))

it has by default a value 10, and seems to determine the range around the dateline (default: +/- 5 degrees) between which geometries to split are expected to lie. The thing is, if you continue this experiment, at some stage the shortest line will no longer cross the dateline.

ogr2ogr describes option -datelineoffset as " offset from dateline in degrees (default long. = +/- 10deg, geometries within 170deg to -170deg will be split)". That seems to be wrong (10 -> +/- 5 degrees, although I didn't actually try with ogr2ogr).

[edzer]

Please report back whether this is still an open issue to you, or whether we can close.

[jmlondon]

I think from an sf perspective, this can be closed. It might be a nice feature to make the arguments to sf::st_wrap_dateline for setting the offset more user friendly. If I find some time, I can take a look at the function and submit a pull request.

I still, however, can't find a solution for adding an sf object that wraps 180 to a leaflet map. The most reliable option is still to cast to sp and thenspTransform with +lon_wrap=180 before adding to leaflet. That's probably more of a leaflet issue and I'll try to find some time to raise the issue over there. The sf + leaflet workflow for quickly visualizing spatial data is a great improvement ... if your work doesn't straddle 180 ;).

[edzer]

See the example in Ops.sfg:

library(sf)
library(maps) 
w = st_as_sf(map('world', plot = FALSE, fill = TRUE))
w2 = (st_geometry(w) + c(360,90)) %% c(360) - c(0,90)
plot(w2, axes = TRUE)

(which doesn't look nice, but does the -180,180 -> 0,360, afaics)

[johnForne]

Kia ora

I'm working with a sf object polygon that crosses the 180 line. I used to be able to run processes, such as intersecting with another polygon also straddling the 180 line but converting to sp:: object and then using sp::spTransform("+proj=longlat +datum=WGS84 +lon_wrap=180"). However, I've just discovered that if I re-run sp::spTransform("+proj=longlat +datum=WGS84 +lon_wrap=180") it no longer seems to transform that the extent from being described in terms of -180W to 180 E to 0 - 360.

I thought it might have something to do with the fact that I was no longer using the same version of sp:: So I rolled back from 2.1-4 to 1.6-1 and then 2.0-0. However, this rolling back didn't seem to resolve the issue.

I then noticed a note that sp package is now running under evolution status 2 (evolution status 2 uses the sf package instead of rgdal). Which lead me (via Google) to sp evolution status: examples of migration from retiring packages. I might be missing something but this report seems to suggest that versions prior to 2.0-0 should still be using rgda which directly calls routines in proj.4.

library(tidyverse)
library(sf)
library(sp)
library(leaflet)
library(dggridR)

# create tibble with extent of the bbox
e <- tibble(xmin = 156.9167, xmax = 193.0833, ymin = -57.58333, ymax = -23.91667)

# use e to create a sf object 
nz_bbox <-
  rbind(
    c(e$xmin, e$ymax),
    c(e$xmin, e$ymin),
    c(e$xmax, e$ymin),
    c(e$xmax, e$ymax),
    c(e$xmin, e$ymax)
  ) %>%
  list() %>%
  st_polygon() %>%
  st_sfc(., crs = "+proj=longlat +datum=WGS84 +no_defs") %>%
  st_sf()

nz_bbox %>% 
  leaflet() %>% 
  addTiles() %>% 
  addPolygons(color = "black")  # this plots nicely over the 180 line because leaflet:: respects 0-360 notation see map below

# write to temp folder
nz_bbox %>% 
  st_write(paste0(tempdir(), "/poly.shp"))

# create example grid
dggs3000 <- dgconstruct(area = 3000, metric = TRUE, resround='nearest')

# now subset the dggs to the area of the supplied rasters----
nz_grid_3000 <- dgshptogrid(dggs3000, "/tmp/RtmptHGBme/poly.shp", cellsize = 0.01) # this sf object seems to retain the 0-360 notation from the 'poly.shp' object.

# Reading layer `poly' from data source `/tmp/RtmptHGBme' using driver `ESRI Shapefile'
# Simple feature collection with 1 feature and 1 field
# Geometry type: POLYGON
# Dimension:     XY
# Bounding box:  xmin: 156.9167 ymin: -57.58333 xmax: 193.0833 ymax: -23.91667
# CRS:           4326

nz_grid_3000 %>% 
  leaflet() %>% 
  addTiles() %>% 
  addPolygons(color = "black")  # for some reason this does not plot well over the 180 line despite initially appearing to use the 0-360 notation (see above)

# Interestingly having plotted 'nz_grid_3000' in a leaflet map for some reason this object no longer has the 0-360 notation??? (see below)

# > nz_grid_3000
# Simple feature collection with 4498 features and 1 field
# Geometry type: POLYGON
# Dimension:     XY
# Bounding box:  xmin: -179.9935 ymin: -58.13773 xmax: 179.9937 ymax: -23.42187
# CRS:           EPSG:4326

# If I then try to intersect the grid with the bbox polygon to "smooth the jagged edge of the grid then the output 'nz_grid_0300_inter' object has -180 to 180 notation... even if I try setting back to 0-360 notation...
nz_grid_3000_inter <- nz_grid_300 %>% 
  as("Spatial") %>%
  sp::spTransform("+proj=longlat +datum=WGS84 +lon_wrap=180") %>%
  st_as_sf() %>% 
  st_intersection(nz_bbox)

# Simple feature collection with 38473 features and 1 field
# Geometry type: POLYGON
# Dimension:     XY
# Bounding box:  xmin: -179.9997 ymin: -57.78281 xmax: 179.9991 ymax: -23.91667
# CRS:           +proj=longlat +datum=WGS84 +no_defs 

# And if I plot it using leaflet:: then it isn't wrapped over the 180 line...
nz_grid_3000_inter %>% 
  leaflet() %>% 
  addTiles() %>% 
  addPolygons(color = "black") 

nz_grid_3000_inter %>% 
  leaflet() %>% 
  addTiles() %>% 
  addPolygons(color = "black") %>% 
  addPolygons(data = nz_bbox)

In sum, can you please clarify why converting to sp and then spTransforming using lon_wrap=180 not longer seems to work. And why rolling back to an earlier version (<2.0-0) doesn't seem to get it working? Given my old trick no longer seems to be working I'd really appreciate any help with how to run processes on polygons that straddle the 180 line? E.g. what can I do to st_intersection(nz_grid_300, nz_bbox)???

Thanks heaps in advance,

John

[edzer]

# Interestingly having plotted 'nz_grid_3000' in a leaflet map for some reason this object no longer has the 0-360 notation??? (see below)

I don't think nz_grid_3000 ever had the 0-360 range, e.g.

> dgshptogrid(dggs3000, "/tmp/RtmpT2lIEP/poly.shp", cellsize = 0.1)
Reading layer `poly' from data source `/tmp/RtmpT2lIEP' using driver `ESRI Shapefile'
Simple feature collection with 1 feature and 1 field
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 156.9167 ymin: -57.58333 xmax: 193.0833 ymax: -23.91667
Geodetic CRS:  +proj=longlat +datum=WGS84 +no_defs
Simple feature collection with 4500 features and 1 field
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -179.9935 ymin: -58.20392 xmax: 179.9937 ymax: -23.4396
Geodetic CRS:  WGS 84
# A tibble: 4,500 × 2
   seqnum                                                               geometry
    <dbl>                                                          <POLYGON [°]>
 1 117614 ((-166.5671 -57.1074, -166.8732 -56.85598, -167.4115 -56.79074, -167.…
 2 117776 ((-166.5123 -57.56449, -166.8261 -57.31386, -167.396 -57.26755, -167.…
 3 117857 ((-167.352 -57.72444, -167.6614 -57.47217, -168.2067 -57.40134, -168.…
 4 171724 ((156.6049 -57.11113, 156.2565 -57.34227, 156.3855 -57.62739, 156.868…
 5 171805 ((157.2158 -57.44764, 156.8685 -57.68068, 157.0037 -57.96542, 157.491…
 6 171806 ((156.819 -56.59428, 156.4764 -56.82633, 156.6049 -57.11113, 157.0813…
 7 171807 ((156.4413 -55.74068, 156.1033 -55.97183, 156.2256 -56.25668, 156.691…
 8 171886 ((157.4228 -56.92929, 157.0813 -57.16322, 157.2158 -57.44764, 157.697…
 9 171887 ((157.028 -56.07676, 156.6911 -56.30977, 156.819 -56.59428, 157.289 -…
10 171888 ((156.6519 -55.22392, 156.3194 -55.45609, 156.4413 -55.74068, 156.900…
# ℹ 4,490 more rows
# ℹ Use `print(n = ...)` to see more rows

the first bounding box you see is reported by st_read-ing the shapfile, the second from the object returned by dgshptogrid. Maybe this is an issue with dggridR?

[johnForne]

Thanks Edzer

Good to clarify that the two different bounding boxes come from different routes.

Given your suggestion that this issue could have something to do with 'dggridR::', I have tried with polygons that were not created using 'dggridR::' and I still seem to run into the issue...

If I use the 'nz_bbox' polygon from above I still don't seem to be able to process this polygon in a way that respects the 180 line. For example, if I use a negative buffer to shrink the polygon, it returns a polygon that wraps the wrong way wrong the 180 line.

So I'm not sure that this is just a 'dggridR' issue?

I'm interested whether it is still possible to even use the ' as("Spatial") %>%
sp::spTransform("+proj=longlat +datum=WGS84 +lon_wrap=180") %>%
st_as_sf()' type trick that @jmlondon mentions above?

nz_bbox_sub <- nz_bbox %>% 
  st_transform(2193) %>% 
  st_buffer(-1) %>% 
  st_transform(4326)

nz_bbox %>% 
  as("Spatial") %>%
  sp::spTransform("+proj=longlat +datum=WGS84 +lon_wrap=180") %>%
  st_as_sf() %>% 
  st_wrap_dateline(c("WRAPDATELINE=YES")) %>% # this doesn't seem to make any difference either
  leaflet() %>% 
  addTiles() %>% 
  addPolygons(data = nz_bbox_sub, color = "red")

[johnForne]

p.s. just discovered that ' st_shift_longitude()' seems to work and provides a nicer way of keeping the objects as sf (rather than having to convert to sp).