st_union reverses polygon orientation

[finnlindgren]

Describe the bug

When using st_union on two "upwards" facing polygons (i.e. in CCW order), the resulting polygon is changed so that it faces downwards (CW).

From section 6.1.11.1 of version 1.2.1 for the SF specification from https://www.ogc.org/standards/sfa

The exterior boundary LinearRing defines the “top” of the surface which is the
side of the surface from which the exterior boundary appears to traverse the
boundary in a counter clockwise direction. The interior LinearRings will
have the opposite orientation, and appear as clockwise when viewed from the “top”,

To Reproduce

# sf 1.0-9 and 1.0-10

# Plain polygon coordinates, CCW order
coord <- matrix(c(0, 1, 1, 0, 0, 0, 0, 1, 1, 0),5,2)
coord
#      [,1] [,2]
# [1,]    0    0
# [2,]    1    0
# [3,]    1    1
# [4,]    0    1
# [5,]    0    0

# Create polygon
pl <- sf::st_polygon(list(coord))

# ok, still CCW:
pl
# POLYGON ((0 0, 1 0, 1 1, 0 1, 0 0))

# ok, still CCW:
sf::st_union(pl)
# POLYGON ((0 0, 1 0, 1 1, 0 1, 0 0))

# incorrect, CW:
sf::st_union(pl, pl)
# POLYGON ((0 0, 0 1, 1 1, 1 0, 0 0))

Additional context

I've been working on conversion methods between sf classes and fmesher mesh storage classes in https://github.com/inlabru-org/inlabru/ (including INLA and the coming fmesher package that will take over the geometry handling from the INLA package), and was encountering problems when testing polygons with holes as well as conversion to and from sp, and traced at least one of the issues to this st_union inconsistency. The fmesher code handles general 3D triangles as well as a boundary segment specification that in 2D can store the complement of regular polygons (i.e. "free holes"), by using the a CCW traversion convention just as the quoted sf standard above (and sometimes without mandating closed sequences; we make use of non-closed "proto boundary" segments as well). In cases where we need to apply st_* methods in geospatial contexts we rely on ordering being handled consistently, so that we don't have to do all the special handling that was needed for sp objects with the ringDir issues etc.

I had previously discovered that the sp conversion method in st_as_sfc.SpatialPolygons is broken in sf 1.0-9 for multi-polygons (in a set of disjoint polygons all but the first polygons were removed by the conversion, as the "make_valid" code removed them, presumably for believing them to be holes, and reversing the ordering of the first polygon. That bug appears to have been fixed in 1.0-10 (at least it doesn't remove polygons) so I haven't filed an Issue for that problem. The related issues/discussion I recall seeing about sp/sf conversion issues focused on geometries read from existing files, but since fmesher is a meshing library, we are the ones creating the geometries, so need to know that we can rely on well-defined behaviour.

Session info

> sessionInfo()

R version 4.2.2 (2022-10-31)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 22.04.1 LTS

Matrix products: default
BLAS:   /home/flindgre/local/R-4.2.2/lib/R/lib/libRblas.so
LAPACK: /home/flindgre/local/R-4.2.2/lib/R/lib/libRlapack.so

locale:
 [1] LC_CTYPE=en_GB.UTF-8       LC_NUMERIC=C               LC_TIME=en_GB.UTF-8        LC_COLLATE=en_GB.UTF-8    
 [5] LC_MONETARY=en_GB.UTF-8    LC_MESSAGES=en_GB.UTF-8    LC_PAPER=en_GB.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C             LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] sf_1.0-10     reprex_2.0.2  usethis_2.1.6

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.10        rstudioapi_0.14    knitr_1.42         magrittr_2.0.3     units_0.8-1        tidyselect_1.2.0  
 [7] prompt_1.0.1       here_1.0.1         R6_2.5.1           rlang_1.0.6        fansi_1.0.4        dplyr_1.1.0       
[13] tools_4.2.2        grid_4.2.2         xfun_0.37          utf8_1.2.3         KernSmooth_2.23-20 cli_3.6.0         
[19] e1071_1.7-13       DBI_1.1.3          withr_2.5.0        class_7.3-21       rprojroot_2.0.3    tibble_3.1.8      
[25] lifecycle_1.0.3    purrr_1.0.1        vctrs_0.5.2        fs_1.6.0           glue_1.6.2         proxy_0.4-27      
[31] compiler_4.2.2     pillar_1.8.1       generics_0.1.3     classInt_0.4-8     pkgconfig_2.0.3 

> sf::sf_extSoftVersion()
          GEOS           GDAL         proj.4 GDAL_with_GEOS     USE_PROJ_H           PROJ 
      "3.10.2"        "3.4.1"        "8.2.1"         "true"         "true"        "8.2.1" 

[edzer]

Thanks for the clear issue! I guess this is caused by the upstream GEOS library (@dr-jts?) - when using s2geometry things are fine:

> p0 = st_sfc(pl)
> st_union(p0, p0) # no CRS, assume R2, use GEOS
Geometry set for 1 feature 
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 0 ymin: 0 xmax: 1 ymax: 1
CRS:           NA
POLYGON ((0 0, 0 1, 1 1, 1 0, 0 0))
> pg = st_sfc(pl, crs = 4326) # will use s2 by default
> st_union(pg, pg)
Geometry set for 1 feature 
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 0 ymin: 0 xmax: 1 ymax: 1
Geodetic CRS:  WGS 84
POLYGON ((0 0, 1 0, 1 1, 0 1, 0 0))

A way to "fix" this in current sf is pulling it through st_sfc:

> sf::st_union(p0, p0) |> st_sfc(check_ring_dir=TRUE)
Geometry set for 1 feature 
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 0 ymin: 0 xmax: 1 ymax: 1
CRS:           NA
POLYGON ((0 0, 1 0, 1 1, 0 1, 0 0))

[dr-jts]

GEOS uses CW as its canonical orientation for polygons. This was chosen historically before the SFS settled on CCW as its canonical orientation. If CCW orientation is desired then client code needs to enforce that.

And before the question is asked, it's a big job to change JTS/GEOS to adopt CCW orientation, and this is probably a breaking change. But note that JTS/GEOS does not require a particular orientation for input, so its flexible in that way.

[edzer]

Thanks, Martin, that clarifies! We'll put that somewhere in the documentation here.

[edzer]

Are you o.k., @finnlindgren if we leave this like it is and ask you to use st_sfc(check_ring_dir=TRUE) to make ring winding OGC compliant? I think it's a clash between de facto (GEOS) and de jure (OGC) standards.

[finnlindgren]

I added the st_sfc(check_ring_dir = TRUE) to the functions that require consistent orientations (assuming nobody needs/wants to have "downward facing" polygons, which will now be indistinguishable from upward facing polygons) and it seems like a workable workaround. If sf handled polygons in 3D euclidean space it would be a problem however, so we'll likely need to make our own classes to handle 3D triangles for meshing spheres and other manifolds, and only use the sf interface to transform the coordinates pointwise for the geospatial applications (we work in "+proj=geocent" a lot).

Related to this, will 1.0-10 go onto CRAN soon? The st_as_sfc.SpatialPolygons handling of multipolygon sp inputs is quite broken in 1.0-9 but I can see 1.0-10 has much expanded code for that method that appears to have solved the problem; in version 1.0-9 it sometimes remove all but the first polygon, again due to misinterpretation of the polygon orientations. I haven't filed a bug since 1.0-10 appears to have fixed it already.

[edzer]

Related to this, will 1.0-10 go onto CRAN soon?

Hopefully yes, will leave this open until it is accepted.

[edzer]

sf 1.0-10 is now on CRAN.

[idaflik]

Hi,

I continue having the same issue when trying to combine all countries into one polygon, st_sfc(check_ring_dir = TRUE) unfortunately is not fixing the issue for me. i'm trying to get a polygon that combines all landmasses of the earth and display it in d3.js.

library(needs)
needs(rnaturalearth, tidyverse, sf)

continents <- ne_countries(returnclass = "sf")%>%
st_set_precision(1000)%>%
st_make_valid()%>%
st_union()%>%
st_sfc(check_ring_dir=TRUE)

st_write(continents, "output/coastline.geojson")

when uploading it to https://observablehq.com/@bumbeishvili/rewind-geojson, it seems like the polygon is not simply wound the wrong way around but seems to be some kind of a mix-up, since the rewind tool (which in other cases works great) also is not able to produce a correct version of the polygon.

as i keep running into this issue and haven't managed to figure out a solution for months now, any help is much appreciated.