st_polygonize returning empty geometry

[eisack00]

Hi,

I'm having an issue where st_polygonize is returning empty geometry. I've created two linestring grids which I union together which I know intersect and can polygonize using QGIS.

I am running

st_union(
  xgrid,
  ygrid) %>%
st_polygonize()

Where xgrid is:

> Simple feature collection with 7 features and 0 fields
Geometry type: LINESTRING
Dimension:     XY
Bounding box:  xmin: 1393427 ymin: 4916632 xmax: 1393459 ymax: 4916679
Projected CRS: NZGD2000 / New Zealand Transverse Mercator 2000
# A tibble: 7 x 1
                            geometry
                    <LINESTRING [m]>
1 (1393441 4916632, 1393441 4916679)
2 (1393427 4916632, 1393427 4916679)
3 (1393433 4916632, 1393433 4916679)
4 (1393440 4916632, 1393440 4916679)
5 (1393446 4916632, 1393446 4916679)
6 (1393459 4916632, 1393459 4916679)
7 (1393441 4916632, 1393441 4916679)

and ygrid is

Simple feature collection with 7 features and 0 fields
Geometry type: LINESTRING
Dimension:     XY
Bounding box:  xmin: 1393427 ymin: 4916632 xmax: 1393459 ymax: 4916679
Projected CRS: NZGD2000 / New Zealand Transverse Mercator 2000
# A tibble: 7 x 1
                            geometry
                    <LINESTRING [m]>
1 (1393427 4916679, 1393459 4916679)
2 (1393427 4916673, 1393459 4916673)
3 (1393427 4916659, 1393459 4916659)
4 (1393427 4916645, 1393459 4916645)
5 (1393427 4916632, 1393459 4916632)
6 (1393427 4916638, 1393459 4916638)
7 (1393427 4916679, 1393459 4916679)

What is returned is

Simple feature collection with 49 features and 2 fields (with 49 geometries empty)
Geometry type: GEOMETRYCOLLECTION
Dimension:     XY
Bounding box:  xmin: NA ymin: NA xmax: NA ymax: NA
Projected CRS: NZGD2000 / New Zealand Transverse Mercator 2000
# A tibble: 49 x 3
      id  id.1                 geometry
 * <int> <int> <GEOMETRYCOLLECTION [m]>
 1     1     1 GEOMETRYCOLLECTION EMPTY
 2     2     1 GEOMETRYCOLLECTION EMPTY
 3     3     1 GEOMETRYCOLLECTION EMPTY
 4     4     1 GEOMETRYCOLLECTION EMPTY
 5     5     1 GEOMETRYCOLLECTION EMPTY
 6     6     1 GEOMETRYCOLLECTION EMPTY
 7     7     1 GEOMETRYCOLLECTION EMPTY
 8     1     2 GEOMETRYCOLLECTION EMPTY
 9     2     2 GEOMETRYCOLLECTION EMPTY
10     3     2 GEOMETRYCOLLECTION EMPTY
# ... with 39 more rows

thanks

xygrid.zip

[rsbivand]

Please attach the RDS files for xgrid and ygrid (as a zip archive), otherwise it is very hhard to reproduce your case. Since the CRS is projected, this is not s2, but might be differences in GEOS precision settings between sf and QGIS. Does QGIS show its precision settings?

[rsbivand]

What does this do:

library(rgeos)
xgrid_sp <- as(xgrid, "Spatial")
ygrid_sp <- as(ygrid, "Spatial")
xy <- gUnion(xgrid_sp, ygrid_sp) # maybe needs byid=TRUE
xy_p <- gPolygonize(xy)
st_as_sfc(xy_p)

rgeos uses a different precision assumption from sf. See also https://github.com/r-spatial/sf/issues/1855.

[eisack00]

Using rgeos and converting to spatial provides the correct polygon output.

Not sure how to check the precision of QGIS. Overall this is all super helpful information for me thank you all!

• Added zip of rds files for x and y gird.

[rsbivand]

Thanks for the data files. @edzer unfortunately using the approach from #1855 doesn't seem to work, but converting to sp and using rgeos does work. So it is a precision issue, but not so easy to resolve. @eisack00 is the output feature count as expected?

[eisack00]

@rsbivand Output feature count is as expected and converting to sp using regos as a workaround works great. I've noticed a duplicate line for each grid which explains the unusual output feature count. However, I've tested the same operations with duplicates removed but the issue remains.