ST_Intersection returning weird results

[latot]

Hi again @edzer has been a while :)

Today happened something weird, I was trying to use st_intersection from one geometry to a sf object, but it start returning nothing, but if only one feature intersects, all the results becomes only one value! From what I understand and I checked the docs, it should be a pairwise operation.

sample.zip

# All repeated just to do a pairwise intersection
a <- sf::st_read(a, quiet = TRUE)
b <- sf::st_read(b, quiet = TRUE)

# All of them has the same polygon, while is only one in b who intersects
sf::st_intersection(b, a)
Simple feature collection with 6 features and 2 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 670008.4 ymin: 5930778 xmax: 670050.8 ymax: 5930827
Projected CRS: WGS 84 / UTM zone 18S
     b  a                           geom
5   20 10 POLYGON ((670012.4 5930815,...
5.1 20 10 POLYGON ((670012.4 5930815,...
5.2 20 10 POLYGON ((670012.4 5930815,...
5.3 20 10 POLYGON ((670012.4 5930815,...
5.4 20 10 POLYGON ((670012.4 5930815,...
5.5 20 10 POLYGON ((670012.4 5930815,...

#Empty...
# Same result as sf::st_intersection(b[1:3,], a[1:3,])
sf::st_intersection(b[1:3,], a)
Simple feature collection with 0 features and 2 fields
Bounding box:  xmin: NA ymin: NA xmax: NA ymax: NA
Projected CRS: WGS 84 / UTM zone 18S
[1] b    a    geom
<0 rows> (o 0- extensión row.names)

# The element 5 in b now intersects, and fill all with "a" instead intersect it
# Similar result as sf::st_intersection(b[1:5,], a[1:5,]) (changes the number of rows that shows)
sf::st_intersection(b[1:5,], a)
Simple feature collection with 6 features and 2 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 670008.4 ymin: 5930778 xmax: 670050.8 ymax: 5930827
Projected CRS: WGS 84 / UTM zone 18S
     b  a                           geom
5   20 10 POLYGON ((670012.4 5930815,...
5.1 20 10 POLYGON ((670012.4 5930815,...
5.2 20 10 POLYGON ((670012.4 5930815,...
5.3 20 10 POLYGON ((670012.4 5930815,...
5.4 20 10 POLYGON ((670012.4 5930815,...
5.5 20 10 POLYGON ((670012.4 5930815,...

Here is two issues:

• Is not retuning the intersection, just one particular geometry all the time
• With some... sets of data, it returns all empty, and if intersects return all with one geometry

I'm very confused, no idea why this happens, I tested it on CRAN and Git versions.

sessionInfo()
R version 4.4.1 (2024-06-14)
Platform: x86_64-pc-linux-gnu
Running under: Gentoo Linux

Matrix products: default
BLAS:   /usr/lib64/libblas.so.3.12.0 
LAPACK: /usr/lib64/liblapack.so.3.12.0

locale:
 [1] LC_CTYPE=es_CL.utf8       LC_NUMERIC=C             
 [3] LC_TIME=es_CL.utf8        LC_COLLATE=es_CL.utf8    
 [5] LC_MONETARY=es_CL.utf8    LC_MESSAGES=es_CL.utf8   
 [7] LC_PAPER=es_CL.utf8       LC_NAME=C                
 [9] LC_ADDRESS=C              LC_TELEPHONE=C           
[11] LC_MEASUREMENT=es_CL.utf8 LC_IDENTIFICATION=C 

Thx!

[Dory111]

Hi,

I'm not the user you tagged, but I came across this post when I was dealing with the same issue. For me it was that my geometry set was a series of linestrings that i cast to polygons. This created overlapping points that st_intersection did not like, but also didn't throw an error for instead returning empty geometries. I fixed this by throwing st_make_valid() around my geometry. Maybe this will work for you as well.

If you found another solution I'd be interested in hearing it as well.

[latot]

Hi! thx for the suggestion, I was using a lot of data, so I moved from R to Postgis and simplify all the workflow..... I did not found any other solution for this......