Problem with st_transform_proj

[Nowosad]

I've tested the st_transform_proj() function and found out that it doesn't work properly with some of the projections, including Winkel Tripel, Boggs Eumorphic, or Natural Earth.
For these projections, the function creates a new object (without a warning or error) that can printed. However, it is impossible to plot the new object:

library(spData)
library(sf)
#> Linking to GEOS 3.5.0, GDAL 2.1.2, proj.4 4.9.2
library(lwgeom)
#> Linking to liblwgeom 2.5.0dev r16016
world_wintri = st_transform_proj(world, crs = "+proj=wintri")
head(world_wintri)
#> Simple feature collection with 6 features and 10 fields
#> geometry type:  MULTIPOLYGON
#> dimension:      XY
#> bbox:           xmin: -6013512 ymin: -6276076 xmax: 6211408 ymax: 4761920
#> epsg (SRID):    NA
#> proj4string:    +proj=wintri
#>   iso_a2            name_long     continent region_un       subregion
#> 1     AF          Afghanistan          Asia      Asia   Southern Asia
#> 2     AO               Angola        Africa    Africa   Middle Africa
#> 3     AL              Albania        Europe    Europe Southern Europe
#> 4     AE United Arab Emirates          Asia      Asia    Western Asia
#> 5     AR            Argentina South America  Americas   South America
#> 6     AM              Armenia          Asia      Asia    Western Asia
#>                type   area_km2      pop  lifeExp gdpPercap
#> 1 Sovereign country  652270.07 31627506 60.37446  1844.022
#> 2 Sovereign country 1245463.75 24227524 52.26688  6955.960
#> 3 Sovereign country   29694.80  2893654 77.83046 10698.525
#> 4 Sovereign country   79880.74  9086139 77.36817 63830.700
#> 5 Sovereign country 2784468.59 42980026 76.15861 18872.967
#> 6 Sovereign country   28656.60  3006154 74.67571  7706.133
#>                             geom
#> 1 MULTIPOLYGON (((5108567.494...
#> 2 MULTIPOLYGON (((1484047.791...
#> 3 MULTIPOLYGON (((1663361.831...
#> 4 MULTIPOLYGON (((4520545.763...
#> 5 MULTIPOLYGON (((-4723076.95...
#> 6 MULTIPOLYGON (((3531966.010...
plot(world_wintri)
#> Warning: plotting the first 9 out of 10 attributes; use max.plot = 10 to
#> plot all
#> OGR: Corrupt data
#> Error in CPL_gdal_dimension(st_geometry(x), NA_if_empty): OGR error

world_boggs = st_transform_proj(world, crs = "+proj=boggs")
plot(world_boggs$geom)
#> OGR: Corrupt data
#> Error in CPL_gdal_dimension(st_geometry(x), NA_if_empty): OGR error

world_natearth = st_transform_proj(world, crs = "+proj=natearth")
plot(world_natearth$geom)
#> OGR: Corrupt data
#> Error in CPL_gdal_dimension(st_geometry(x), NA_if_empty): OGR error

[edzer]

Yes, this is nasty -- before plotting we need to filter out empty geometries, which is currently done by !is.na(st_dimension(x)) and that takes the gdal path. There should be other ways to solve that without gdal, e.g. GEOSisEmpty_r in GEOS, but that would then call for another function.

The current workaround is to set

st_crs(world_wintri) = NA

before you plot.

[Nowosad]

In the new lwgeom version it is possible to plot a geometry in e.g. the Winkel tripel projection. However, the graticule argument doesn't work:

library(spData)
library(sf)
#> Linking to GEOS 3.5.0, GDAL 2.1.2, proj.4 4.9.2
library(lwgeom)
#> Linking to liblwgeom 2.5.0dev r16016

world_wintri = st_transform_proj(world, crs = "+proj=wintri")
plot(world_wintri$geom)

plot(world_wintri$geom, graticule = TRUE)
#> OGR: Corrupt data
#> Error in CPL_crs_equivalent(e1$proj4string, e2$proj4string): OGR error

[edzer]

See r-spatial/sf#509