st_intersection (with missing y) produces "degenerate" polygons

[jguelat]

Describe the bug
I'm using the st_intersection() function on a single sf object with overlapping polygons to extract the geometries with and without overlap. The polygons of the overlapping parts are OK, but the non-overlapping polygons almost always have some zero-width "part".

To Reproduce
Using the attached data, polygons 1 and 3 have a zero-width part while polygon 2 is OK.
testinter.zip

library(sf)
pols <- st_read("testinter.gpkg", "pols")
inter <- st_intersection(pols)
par(mfrow = c(1, 3))
plot(st_geometry(inter[1,]))
plot(st_geometry(inter[2,]))
plot(st_geometry(inter[3,]))

Additional context
Thanks a lot for all your work on the sf package!

sessionInfo()

R version 4.4.2 (2024-10-31 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 10 x64 (build 19045)

Matrix products: default

locale:
[1] LC_COLLATE=German_Switzerland.utf8 LC_CTYPE=German_Switzerland.utf8 LC_MONETARY=German_Switzerland.utf8
[4] LC_NUMERIC=C LC_TIME=German_Switzerland.utf8

time zone: Europe/Zurich
tz

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

other attached packages:
[1] sf_1.0-19

loaded via a namespace (and not attached):
[1] compiler_4.4.2 magrittr_2.0.3 class_7.3-23 tools_4.4.2 DBI_1.2.3 units_0.8-5
[7] proxy_0.4-27 Rcpp_1.0.14 KernSmooth_2.23-26 grid_4.4.2 e1071_1.7-16 classInt_0.4-11

sf::sf_extSoftVersion()

      GEOS           GDAL         proj.4 GDAL_with_GEOS     USE_PROJ_H           PROJ 
  "3.12.2"        "3.9.3"        "9.4.1"         "true"         "true"        "9.4.1" 

[edzer]

I get

library(sf)
# Linking to GEOS 3.12.2, GDAL 3.9.3, PROJ 9.4.1; sf_use_s2() is TRUE
pols <- st_read("testinter.gpkg", "pols")
# Reading layer `pols' from data source `/home/edzer/Downloads/testinter.gpkg' using driver `GPKG'
# Simple feature collection with 2 features and 1 field
# Geometry type: POLYGON
# Dimension:     XY
# Bounding box:  xmin: 2657047 ymin: 1217612 xmax: 2658977 ymax: 1219660
# Projected CRS: CH1903+ / LV95
inter <- st_intersection(pols)
par(mfrow = c(1, 3))
plot(st_geometry(inter[1,]))
plot(st_geometry(inter[2,]))
plot(st_geometry(inter[3,]))
inter
# Simple feature collection with 3 features and 3 fields
# Geometry type: POLYGON
# Dimension:     XY
# Bounding box:  xmin: 2657047 ymin: 1217612 xmax: 2658977 ymax: 1219660
# Projected CRS: CH1903+ / LV95
#     id n.overlaps origins                       geometry
# 1    1          1       1 POLYGON ((2657281 1218360, ...
# 1.1  1          2    1, 2 POLYGON ((2658787 1218444, ...
# 2    2          1       2 POLYGON ((2658977 1218282, ...

which looks OK?

[jguelat]

Thanks for your answer! The output has indeed 3 elements (one overlap and 2 non-overlapping ones), but if you look at the first and third polygons, you'll see that they both have something like a "line extension" (on the right side) that should not be here.

[edzer]

The geometry calculations are done at the level of the GEOS upstream library. It is of course possible that the implementation in sf for n-ary intersection is broken.

[jguelat]

Thanks! I was also suspecting a GEOS problem... I wonder how QGIS is solving it. When using the Union Processing algorithm (the name is confusing but it does the same as st_intersection), the output polygons don't have these "line extensions". Maybe there's some sort of post-processing cleaning happening in the background.

[edzer]

Thank you - that's something to look into!