-
Notifications
You must be signed in to change notification settings - Fork 299
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
sf 1.0 will switch to using S2 geometry for ellipsoidal coordinates #1649
Comments
Dear package author, while checking reverse dependency's against the sf 1.0 candidate which uses s2 as a geometry back-end when using ellipsoidal coordinates, we found a regression. Please verify that this is not a false positive, and update your package on CRAN in a timely fashion before See #1649 for details. |
For those interested in a low-res world map valid on S2, created by @paleolimbot : library(sf)
# Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1
library(s2)
st_as_sf(s2_data_tbl_countries)
# Simple feature collection with 177 features and 2 fields
# Geometry type: MULTIPOLYGON
# Dimension: XY
# Bounding box: xmin: -180 ymin: -85.60904 xmax: 180 ymax: 83.64513
# CRS: NA
# First 10 features:
# name continent
# 1 Afghanistan Asia
# 2 Angola Africa
# 3 Albania Europe
# 4 United Arab Emirates Asia
# 5 Argentina South America
# 6 Armenia Asia
# 7 Antarctica Antarctica
# 8 French Southern and Antarctic Lands Seven seas (open ocean)
# 9 Australia Oceania
# 10 Austria Europe
# geometry
# 1 MULTIPOLYGON (((61.21082 35...
# 2 MULTIPOLYGON (((16.32653 -5...
# 3 MULTIPOLYGON (((20.59025 41...
# 4 MULTIPOLYGON (((51.57952 24...
# 5 MULTIPOLYGON (((-65.5 -55.2...
# 6 MULTIPOLYGON (((43.58275 41...
# 7 MULTIPOLYGON (((-59.57209 -...
# 8 MULTIPOLYGON (((68.935 -48....
# 9 MULTIPOLYGON (((145.398 -40...
# 10 MULTIPOLYGON (((16.97967 48...
|
For those getting errors involving redundant, invalid, or crossing edges, it is likely because s2 assumes the shortest distance on the sphere between two points (not along a constant latitude or longitude as would appear in a standard cartesian plot with points connected). A few things worth trying are:
I'm certainly happy to PR some of this in (but @edzer should decide how/if it's done since he's spent more time in sf!) |
We can't call library(tmap)
library(sf)
# Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1
library(s2)
data(World)
x = st_transform(World, 4326)
# qtm(x) # -->> fails
w = s2_geog_from_wkb(st_as_binary(st_geometry(x)))
# Error in s2_geography_from_wkb(wkb_bytes, oriented = oriented, check = check) :
# Loop 7 is not valid: Edge 171 crosses edge 174
# Calls: s2_geog_from_wkb -> new_s2_xptr -> s2_geography_from_wkb
# Execution halted |
And the following also didn't help: w = s2_geog_from_wkb(st_as_binary(st_segmentize(st_geometry(x), units::set_units(10, km))))
w = s2_geog_from_wkb(st_as_binary(st_segmentize(st_geometry(x), units::set_units(1, km))))
w = s2_geog_from_wkb(st_as_binary(st_segmentize(st_geometry(x), units::set_units(.1, km))))
g = st_geometry(x)
st_crs(g) = NA # so we segmentize in "degrees"
w = s2_geog_from_wkb(st_as_binary(st_segmentize(g, .01)))
w = s2_geog_from_wkb(st_as_binary(st_segmentize(g, .1)))
w = s2_geog_from_wkb(st_as_binary(st_segmentize(g, .001)))
w = s2_geog_from_wkb(st_as_binary(st_segmentize(g, .0001))) |
I think the key is library(tmap)
#> Warning: package 'tmap' was built under R version 4.0.3
library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
# Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1
library(s2)
#> Warning: package 's2' was built under R version 4.0.5
data(World)
x = st_transform(World, 4326)
# qtm(x) # -->> fails
w = s2_geog_from_wkb(st_as_binary(st_geometry(x)), check = FALSE)
library(s2plot) # remotes::install_github("paleolimbot/s2plot")
s2plot(s2_rebuild(w)) Created on 2021-04-19 by the reprex package (v2.0.0) |
Bingo! This makes |
I think so - it's something like a Union. I think it might be worth exposing As a note, when I was testing with BigQuery, I did notice that they do this operation when taking any foreign input (wkt, wkb). I think they also round coordinates on input to snap level 30. |
That might be a bit hard: you'd have to pass it e.g. in every Are there any risks on the way back: S2 geometries converted to wkt, then to sf, then e.g. projected to whatever? |
just to confirm – 6de3420 fixes all errors with polygons from |
The "risk" is mostly that it's slow and shouldn't have to be. I think that any operation that returns a gemoetry (intersection, etc.) will probably be unaffected but I forget the details and I'm worried that it could do things like change the direction of a line (I guess that a GEOS overlay operation might do this too though). The overlay ops might be much faster than GEOS though so users might not notice? Two ways to provide an escape hatch here without running all input through
I think these are along the same lines as how users currently have to handle invalid input. In the meantime s2 should provide an equivalent to |
Re: |
FWIW, here are the failing polygons from the most recent #> Simple feature collection with 6 features and 6 fields
#> Geometry type: MULTIPOLYGON
#> Dimension: XY
#> Bounding box: xmin: -180 ymin: 8.229188 xmax: 180 ymax: 81.2504
#> CRS: +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
#> dataset name st_valid fixed_prec fixed_prep_script fixed_lwgeom
#> 1 ne_countries_110 Sudan TRUE FALSE FALSE FALSE
#> 2 ne_countries_110 Russia TRUE FALSE TRUE FALSE
#> 3 ne_map_units_110 Sudan TRUE FALSE FALSE FALSE
#> 4 ne_map_units_110 Russia TRUE FALSE TRUE FALSE
#> 5 ne_sovereignty_110 Sudan TRUE FALSE FALSE FALSE
#> 6 ne_sovereignty_110 Russia TRUE FALSE TRUE FALSE For |
Note that Centos 7 runs very old PROJ/GDAL/GEOS, so you need also to handle those. The CRAN Solaris instance runs moderately old PROJ (5.2.0), GDAL (2.2.4) and GEOS I think 3.7, but may be 3.6. It is very useful in detecting packages that do not condition for PROJ versions. |
I have good experience with world dataset from It comes from the Eurostat / via GISCO and seems to be on a NC license, don't know if that is an issue... |
Thanks @jlacko That is a good source indeed. I work in official statistics, so Eurostat has my sympathy vote. I'll continue with sharing details on the 'new world' dataset in |
Glad to be of service! Eurostat are good guys, even though the giscoR package is not officially endorsed by them and just kind of sort of hijacks their API, but more importantly (against current sf / s2 versions):
|
Clever: > giscoR::gisco_get_countries() |> st_bbox()
xmin ymin xmax ymax
-180.00000 -89.90000 179.99999 83.65187 they created a small gap at the antimeridean. |
I have pondered how to update the library(sf)
#> Linking to GEOS 3.9.0, GDAL 3.2.2, PROJ 7.2.1
library(spData)
library(giscoR)
sf_use_s2(FALSE)
#> Spherical geometry (s2) switched off
bb = st_bbox(giscoR::gisco_get_countries())
world2 = st_intersection(world, st_as_sfc(bb))
#> although coordinates are longitude/latitude, st_intersection assumes that they are planar
#> Warning: attribute variables are assumed to be spatially constant throughout all
#> geometries
sf_use_s2(TRUE)
#> Spherical geometry (s2) switched on
isv = sf::st_is_valid(world2)
world2$name_long[!isv]
#> [1] "Sudan"
world3 = st_make_valid(world2)
plot(world3[,1]) Created on 2021-06-13 by the reprex package (v2.0.0) @edzer @paleolimbot is that make any sense? Can I update the spData package with |
Almost there! I get that Sudan is not valid in library(sf)
#> Linking to GEOS 3.8.1, GDAL 3.2.1, PROJ 7.2.1
library(spData)
#> To access larger datasets in this package, install the spDataLarge
#> package with: `install.packages('spDataLarge',
#> repos='https://nowosad.github.io/drat/', type='source')`
library(giscoR)
sf_use_s2(FALSE)
#> Spherical geometry (s2) switched off
bb = st_bbox(giscoR::gisco_get_countries())
world2 = st_intersection(world, st_as_sfc(bb))
#> although coordinates are longitude/latitude, st_intersection assumes that they are planar
#> Warning: attribute variables are assumed to be spatially constant throughout all
#> geometries
sf_use_s2(TRUE)
#> Spherical geometry (s2) switched on
isv = sf::st_is_valid(world2)
world2$name_long[!isv]
#> [1] "Sudan"
world3 <- world2
world3$geom <- world3$geom %>%
s2::as_s2_geography(check = F) %>%
s2::s2_union() %>%
s2::s2_rebuild(s2::s2_options(split_crossing_edges = TRUE)) %>%
st_as_sfc()
plot(world3[,1]) isv = sf::st_is_valid(world3)
world2$name_long[!isv]
#> character(0)
# experimental pkg, I just use it to check s2 objects
s2plot::s2plot(world3, col = "grey90", projection = s2plot::s2plot_projection_orthographic("POINT (0 0)")) s2plot::s2plot(world3, col = "grey90", projection = s2plot::s2plot_projection_orthographic("POINT (180 0)")) s2plot::s2plot(world3, col = "grey90", projection = s2plot::s2plot_projection_orthographic("POINT (0 90)")) s2plot::s2plot(world3, col = "grey90", projection = s2plot::s2plot_projection_orthographic("POINT (0 -90)")) Created on 2021-06-13 by the reprex package (v0.3.0) |
Thank you a lot @paleolimbot ! I have already added the changes to the dev version of spData (Nowosad/spData@957c1db) and plan to update it on CRAN sometimes this week. |
Thanks @paleolimbot and @Nowosad! I've healed tmap's library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
library(spData)
#> To access larger datasets in this package, install the spDataLarge
#> package with: `install.packages('spDataLarge',
#> repos='https://nowosad.github.io/drat/', type='source')`
library(giscoR)
library(tmap)
sf_use_s2(FALSE)
#> Spherical geometry (s2) switched off
# load 'broken' World dataset from tmap
data(World)
bb = st_bbox(giscoR::gisco_get_countries())
W2 = World |>
st_transform(crs = 4326) |>
st_intersection(st_as_sfc(bb))
#> although coordinates are longitude/latitude, st_intersection assumes that they are planar
#> Warning: attribute variables are assumed to be spatially constant throughout all
#> geometries
sf_use_s2(TRUE)
#> Spherical geometry (s2) switched on
isv = sf::st_is_valid(W2)
W2$name[!isv]
#> factor(0)
#> 177 Levels: Afghanistan Albania Algeria Angola Antarctica Argentina ... Zimbabwe
tm_shape(W2) + tm_polygons() tm_shape(W2, projection = "+proj=eck4") + tm_polygons() # old default CRS for World # round earth tests:
tm_shape(W2, projection = "+proj=ortho +lon_0=0 +lat_0=0") + tm_polygons() tm_shape(W2, projection = "+proj=ortho +lon_0=0 +lat_0=90") + tm_polygons() tm_shape(W2, projection = "+proj=ortho +lon_0=0 +lat_0=45") + tm_polygons()
#> Error: Shape contains invalid polygons. Please fix it or set tmap_options(check.and.fix = TRUE) and rerun the plot
# valid?
W2_ortho_0_0 = st_transform(W2, "+proj=ortho +lon_0=0 +lat_0=0")
st_is_valid(W2_ortho_0_0)
#> [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
#> [13] TRUE TRUE TRUE FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
#> [25] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
#> [37] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
#> [49] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
#> [61] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
#> [73] TRUE FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
#> [85] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
#> [97] TRUE TRUE TRUE TRUE TRUE TRUE FALSE TRUE TRUE TRUE TRUE TRUE
#> [109] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
#> [121] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
#> [133] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
#> [145] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
#> [157] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
#> [169] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
W2_ortho_0_90 = st_transform(W2, "+proj=ortho +lon_0=0 +lat_0=90")
st_is_valid(W2_ortho_0_90)
#> [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
#> [13] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE TRUE
#> [25] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE FALSE FALSE
#> [37] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE TRUE
#> [49] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
#> [61] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
#> [73] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
#> [85] FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
#> [97] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
#> [109] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
#> [121] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
#> [133] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
#> [145] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
#> [157] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE TRUE TRUE
#> [169] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
W2_ortho_0_45 = st_transform(W2, "+proj=ortho +lon_0=0 +lat_0=45")
st_is_valid(W2_ortho_0_45)
#> [1] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
#> [26] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
#> [51] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
#> [76] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
#> [101] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
#> [126] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
#> [151] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
#> [176] NA NA
# s2plot works for 0 45:
library(s2plot)
s2plot::s2plot(W2, col = "grey90", projection = s2plot::s2plot_projection_orthographic("POINT (0 45)")) Created on 2021-06-14 by the reprex package (v2.0.0) |
Polygons that are partly out of sight (behind the globe) may come back as invalid: you need to clip them to the part(s) of the polygon that is visible before projecting. I think that is what |
Any update on the invalid polygon error from |
…h remained in spite of st_make_valid(). Related to r-spatial/sf#1649
I love tmap and have been using it constantly for the 30 Day Map Challenge. It has been working great until yesterday. On code that ran great 2 days ago, I now get: Error: Shape contains invalid polygons. Please fix it or set tmap_options(check.and.fix = TRUE) and rerun the plot constantly. After playing with marmap and having it hose my R environmant, especially tidyverse. I had to delete my entire R environment along with libraries. I worked with Jenny Bryan and she helped with the tidyverse issue tidyverse/tidyverse#316 (comment) We agreed about the fresh start with R. But, now, I am getting errors with tmap. I am getting the same errors off your sample vignettes (tmap: version changes). Can you assist? Thanks. Jim Here is my sys info.
R version 4.2.2 (2022-10-31 ucrt) Matrix products: default locale: attached base packages: other attached packages: loaded via a namespace (and not attached): |
Well, a reprex is missing and your issue sounds like a |
As announced here in June last year, and also in the NEWS file around the same time, it is planned that where possible
sf
will use the spherical geometry operators of packages2
instead of theGEOS
routines that assume projected coordinates. I've now run reverse dependency checks and found regressions in the following packages ("check" means: pkg maintainer informed):I have opened issues / written emails to all packages affected (issues mentioned below).
Testing your package while
sf
usess2
can be done by settingor by specifying the environment variable
_SF_USE_S2
to the valuetrue
(case sensitive); the latter can be done from within R withbut that needs to be done before
sf
is loaded.The text was updated successfully, but these errors were encountered: