st_transform +proj=wintri OGR error

[econandrew]

(Good chance this is from screwed up libraries, but I can't find any positive examples that this does work elsewhere so posting anyway.)

Example:

library(sf)

p1 = st_point(c(7,52))
geom.sf = st_sfc(p1, crs = 4326)
x <- st_transform(geom.sf, "+proj=aea")    # Works
x <- st_transform(geom.sf, "+proj=wintri") # Fails

# Equivalent in sp works
geom.sp = as(geom.sf, "Spatial")
wintri <- sp::CRS("+proj=wintri")
x <- sp::spTransform(geom.sp, wintri)

Output (homebrew, OS X El Cap):

> library(sf)
Linking to GEOS 3.6.2, GDAL 2.2.2, proj.4 4.9.3, lwgeom 2.4.0 r15853
> 
> p1 = st_point(c(7,52))
> geom.sf = st_sfc(p1, crs = 4326)
> x <- st_transform(geom.sf, "+proj=aea")    # Works
> x <- st_transform(geom.sf, "+proj=wintri") # Fails
OGR: Corrupt data
Error in CPL_transform(x, crs$proj4string) : OGR error
In addition: Warning message:
In CPL_crs_from_proj4string(x) :
  Cannot import crs from PROJ.4 string `+proj=wintri', missing crs returned
> 
> # Equivalent in sp works
> geom.sp = as(geom.sf, "Spatial")
> wintri <- sp::CRS("+proj=wintri")
> x <- sp::spTransform(geom.sp, wintri)

Any ideas?

[jsta]

Doesn't work for me either. Is it because of this?

While Winkel Tripel projection is defined in the proj library and can be called from the command line, it can't be used as a custom CRS in QGIS because there's no inverse transformation in the proj library.

The enhancement request to add this functionality has been closed since it seems that the inverse transformation isn't trivial

[econandrew]

Maybe, but inverse transformation shouldn't really be an issue here. Glad someone else can reproduce it though.

[edzer]

sf uses the GDAL api for coordinate transformations and conversions, rgdal (which sp calls) uses the proj api directly. My guess is that GDAL insists on inverse projections to be present, and hence does not interface this one.

[edzer]

Searching for "+proj=wintri gdal" gave a few hits suggesting this is a well-known problem, with gdal.

[rsbivand]

... and in rgdal, we have a kludge by @barryrowlingson to permit checking for missing inverses too, so it's even harder to drop something heavy on your foot.

[econandrew]

Ok that explains it. Still, Winkel Tripel may have messy maths, but since Nat Geo uses it for world maps it's pretty common in the wild. I guess in an ideal world GDAL would support numerical inverses (? - the ESRI tools seem to have no problem with wintri, including inverse-using operations).

In the meantime this is my +1 for the rgdal kludge to be applied here if possible.

[edzer]

So the feature request here is to enable st_transform, possibly optional, to directly use the PROJ.4 API rather than using the GDAL API, similar to how rgdal does this.

[edzer]

library(sf)

p1 = st_point(c(7,52))
geom.sf = st_sfc(p1, crs = 4326)
x <- st_transform(geom.sf, "+proj=wintri", use_gdal = FALSE) # Works!

It should be noted that downstream problems may be expected, as we now created an object with a proj.4 string that GDAL does not accept, but if this is the end of the script, there should be no problem.

Also, this only works if sf was built while linking against liblwgeom.

[edzer]

@econandrew can we close this issue?

[econandrew]

Just testing it - your point transformation works fine for me as does +over now, but wintri on map('world') still causes an error for me.

library(ggplot2)
library(maps)

world1 <- sf::st_as_sf(map('world', plot = FALSE, fill = TRUE))

# This (pre-this-issue) doesn't work as it wraps the polygons, which then invert
world2 <- sf::st_transform(world1, crs = "+proj=cea +datum=WGS84 +no_defs +over")
ggplot() + geom_sf(data = world2)

# This does work, excellent!
world3 <- sf::st_transform(world1, crs = "+proj=cea +datum=WGS84 +no_defs +over", use_gdal = FALSE)
ggplot() + geom_sf(data = world3)

# This still fails for me
world4 <- sf::st_transform(world1, crs = "+proj=wintri +datum=WGS84 +no_defs +over", use_gdal = FALSE)
ggplot() + geom_sf(data = world4)
#> OGR: Corrupt data
#> Error in CPL_crs_equivalent(e1$proj4string, e2$proj4string): OGR error

[edzer]

That is a downstream problem I predicted: the graticule asked for by ggplot wants to inverse transform the bounding box of world4, which is known to not work...

I think we found a clear limitation of sf here, compared to sp/rgdal.

[econandrew]

Ah, I see, that makes some sense. Knowing that as the cause, and borrowing from your own answer to tidyverse/ggplot2#2071, the following does work:

library(ggplot2)
library(maps)

world1 <- sf::st_as_sf(map('world', plot = FALSE, fill = TRUE))
world4 <- sf::st_transform(world1, crs = "+proj=wintri +datum=WGS84 +no_defs +over", use_gdal = FALSE)
ggplot() + geom_sf(data = world4) + theme_minimal() + coord_sf(datum = NA) + 
  theme(
    panel.grid = element_blank(), 
    line = element_blank(), 
    rect = element_blank(), 
    text = element_blank())

I think there's possibly still some work to be done on the ggplot side to make some of this smoother (e.g. could use_gdal be carried up to coord_sf), but as far as this package goes I think this is closed now. Thanks!

[edzer]

Cool! Here is the version with a manual graticule:

library(ggplot2)
library(maps)
world1 <- sf::st_as_sf(map('world', plot = FALSE, fill = TRUE))
world4 <- sf::st_transform(world1, crs = "+proj=wintri +datum=WGS84 +no_defs +over", use_gdal = FALSE)
gr = sf::st_graticule(lat = c(-89.9,seq(-80,80,20),89.9))
gr <- sf::st_transform(gr, crs = "+proj=wintri +datum=WGS84 +no_defs +over", use_gdal = FALSE)
ggplot() + geom_sf(data = gr, color = 'grey') + geom_sf(data = world4) + theme_minimal() +
  coord_sf(datum = NA) +
  theme(
    panel.grid = element_blank(),
    line = element_blank(),
    rect = element_blank(),
    text = element_blank())

[edzer]

As per future sf > 0.5-5, use_gdal will be deprecated, and package lwgeom has a function st_transform_proj that does projections entirely through Proj.4 (requiring liblwgeom).

[MartinBeal]

@edzer any idea for how to label these graticules (i.e. lat/lons)? I re-ran (with small edits) your code, and attempted to add labels using the scale_y_continuous() function as per #2857 but nothing happens. (FYI I've switched to Robinsons for simplicity).

Also tried coord_sf(label_graticule="NE") to no avail.

library(ggplot2)
library(maps)

robin <- "+proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"

world1 <- rnaturalearth::ne_countries(scale = "medium", returnclass = "sf") # load country dataset
world4 <- lwgeom::st_transform_proj(world1, crs = robin, use_gdal = FALSE)

gr <- sf::st_graticule(lat = seq(-80,80,40), lon = seq(-180, 180 ,by=40))
gr <- lwgeom::st_transform_proj(gr, crs = robin)

ggplot() + geom_sf(data = gr, color = 'grey') + geom_sf(data = world4) + theme_minimal() +
  coord_sf(datum=NA) +
  coord_sf() +
  scale_x_continuous(
    breaks = c(100, 180),
    labels = letters[1:2])