upcoming sf will break tmap when adopting s2

[edzer]

See r-spatial/sf#1649 and consider commenting there in case this is not trivial to solve. @paleolimbot

> qtm(World, fill="#FFF8DC", projection=4326, inner.margins=0) +
+   tm_grid(x = seq(-180, 180, by=20), y=seq(-90,90,by=10), col = "gray70") +
+   tm_xlab("Longitude") +
+   tm_ylab("Latitude")
Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1
Error: Shape contains invalid polygons. Please fix it or set tmap_options(check.and.fix = TRUE) and rerun the plot
Execution halted

looks like another world map that is not valid on S2?

[mtennekes]

Apparently, S2 isn't very "natural'. World shapes from naturalearth fail:

library(rnaturalearth)
library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1

W = ne_countries(returnclass = "sf")

W_valid = st_is_valid(W)

which(!W_valid)
#> [1]   7  54 136

plot(W[!W_valid, 1])

W2 = st_make_valid(W)
plot(W2[!W_valid,1])

W2_valid = st_is_valid(W2)

which(!W2_valid)
#> integer(0)

^{Created on 2021-04-21 by the reprex package (v1.0.0)}

I guess that this is caused by polygons breaking the (-)180 longitude. Maybe it can be fixed by allowing longitudes to be outside the [-180, 180] range? E.g. by adding 360 to the longitudes of Russia that are < 0. Haven't tried it yet.

[mtennekes]

Forgot to mention that the issue has been fixed: tmap does the st_is_valid after st_transform. It will be made valid when tmap_options(check.and.fix = TRUE).

However, since the World is based on the shape from rnaturalearth, it has to be valid, also after re/unprojection.

[edzer]

What you consider "natural" is topological validity on R2 (plate carrée), not on S2. Since our Earth is round, I consider S2 the more natural interpretation of ellipsoidal coordinates, despite the 6 decades of GIS tradition ignoring this. Someone has to start.

Once valid on S2, there is no easy way back; you'd have to split at the dateline and break open the poles again.

I guess that this is caused by polygons breaking the (-)180 longitude.

That, and the singularity of poles (Antarctica making an excursion through +/-180,-90)

[edzer]

s2 1.0.5 (CRAN) and sf (github, branch sf_1_0) now have a st_make_valid that does a proper job for spherical coordinates. Still failing tmap and tmaptools:

Package: tmap
Check: examples
Old result: OK
New result: ERROR
  Running examples in ‘tmap-Ex.R’ failed
  The error most likely occurred in:
  
  > base::assign(".ptime", proc.time(), pos = "CheckExEnv")
  > ### Name: tm_xlab
  > ### Title: Axis labels
  > ### Aliases: tm_xlab tm_ylab
  > 
  > ### ** Examples
  > 
  > data(World)
  > 
  > qtm(World, fill="#FFF8DC", projection=4326, inner.margins=0) +
  + 	tm_grid(x = seq(-180, 180, by=20), y=seq(-90,90,by=10), col = "gray70") +
  + 	tm_xlab("Longitude") +
  + 	tm_ylab("Latitude")
  Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1
  Error: Shape contains invalid polygons. Please fix it or set tmap_options(check.and.fix = TRUE) and rerun the plot
  Execution halted

Package: tmaptools
Check: examples
Old result: OK
New result: ERROR
  Running examples in ‘tmaptools-Ex.R’ failed
  The error most likely occurred in:
  
  > base::assign(".ptime", proc.time(), pos = "CheckExEnv")
  > ### Name: crop_shape
  > ### Title: Crop shape object
  > ### Aliases: crop_shape
  > 
  > ### ** Examples
  > 
  > if (require(tmap) && packageVersion("tmap") >= "2.0") {
  +     data(World, NLD_muni, land, metro)
  + 
  +     #land_NLD <- crop_shape(land, NLD_muni)
  + 
  +     #qtm(land_NLD, raster="trees", style="natural")
  + 
  +     metro_Europe <- crop_shape(metro, World[World$continent == "Europe", ], polygon = TRUE)
  + 
  +     qtm(World) +
  +     tm_shape(metro_Europe) +
  +     	tm_bubbles("pop2010", col="red", title.size="European cities") +
  +     	tm_legend(frame=TRUE)
  + }
  Loading required package: tmap
  Error in s2_geography_from_wkb(x, oriented = oriented, check = check) : 
    Evaluation error: Found 1 feature with invalid spherical geometry.
  [34] Loop 5 edge 0 crosses loop 6 edge 4.
  Calls: crop_shape ... as_s2_geography.WKB -> new_s2_xptr -> s2_geography_from_wkb
  Execution halted

[Nowosad]

I think there are still some s2 problems in tmap. The below examples show an issue with tm_graticules():

library(sf)
#> Linking to GEOS 3.9.0, GDAL 3.2.2, PROJ 7.2.1
library(tmap)
data("World")
World2 = World
World2$geometry[World$iso_a3 == "USA"] = st_cast(World$geometry[World$iso_a3 == "USA"], "POLYGON")[6]
ortho_crs = st_crs("+proj=ortho +lon_0=11 +lat_0=46")
World_ortho = st_transform(World2, crs = ortho_crs)
World_ortho = st_make_valid(World_ortho)
World_ortho = st_cast(World_ortho, "MULTIPOLYGON")
World_ortho = World_ortho[!st_is_empty(World_ortho), ]

# works
tm_shape(World_ortho) +
  tm_polygons() +
  tm_layout()

# works
sf_use_s2(FALSE)
#> Spherical geometry (s2) switched off
tm_shape(World_ortho) +
  tm_polygons() +
  tm_graticules() +
  tm_layout()

# does not work
sf_use_s2(TRUE)
#> Spherical geometry (s2) switched on
tm_shape(World_ortho) +
  tm_polygons() +
  tm_graticules() +
  tm_layout()
#> Error in st_coordinates.sfc(st_geometry(lns_crop)[sum(lnsSel)]): not implemented for objects of class sfc_GEOMETRYCOLLECTION

^{Created on 2021-07-05 by the reprex package (v2.0.0)}

[mtennekes]

The problem was that the intersection between visual bounding box and graticules resulted in a geometrycollection. If should be solved, but no guarantees that it is s2-bug-free.

[paleolimbot]

If you use s2_intersects() you can specify the output type of the overlay op to only include lines using options = s2_options(dimensions = "polyline")!