st_union (or st_erase) takes ages

[raff-k]

I would like to simply union two polygons. One of them is created by a raster to polygon conversion using RSAGA (I assume this polygon is problematic). However, the st_union-procedure takes ages, while, unfortunately, ArcGIS finishs in seconds (which annoys me a bit!).
In ArcGIS one can set the xy-tolerance (default is 0.001), maybe this speeds up the process. I tried to reduce the precision to 0.001, as well, but the results looks really strange :-/

Here, is the code with the shapefiles:

# load packages
if(!require("pacman")){install.packages("pacman")}
pacman::p_load(sf, lwgeom, curl)

# get data from github
URLshape1 <- ("https://raw.githubusercontent.com/raff-k/DataUpload/master/shape1.rds")
URLshape2 <- ("https://raw.githubusercontent.com/raff-k/DataUpload/master/shape2.rds")

# set pathes
path.shape1 <- file.path(tempdir(), "shape1.rds")
path.shape2 <- file.path(tempdir(), "shape2.rds")

# download data
download.file(url = URLshape1, destfile = path.shape1, method = "curl")
download.file(url = URLshape2, destfile = path.shape2, method = "curl")

# load data
shape1 <- readRDS(file = path.shape1)
shape2 <- readRDS(file = path.shape2) # the "problematic" sf



# 1 SIMPLE UNION 
union.simpl <- sf::st_union(x = shape1, y = shape2) # ... takes ages (same for st_erase)


# 2 UNION: with reduced precision (0.001)
shape1.prec <- shape1 %>% sf::st_set_precision(x = ., precision = 0.001)
shape2.prec <- shape2 %>% sf::st_set_precision(x = ., precision = 0.001)

# ... check geometry
if(!all(sf::st_is_valid(shape1.prec))){shape1.prec <- lwgeom::st_make_valid(shape1.prec)}
if(!all(sf::st_is_valid(shape2.prec))){shape2.prec <- lwgeom::st_make_valid(shape2.prec)}


# ... perform union
union.prec <- sf::st_union(x = shape1.prec, y = shape2.prec) # ... really fast
plot(union.prec) # looks really strange

[edzer]

If your coordinates are in units metre, and you want to round to mm, then you need to use a precision value of 1000, not 0.001. A value of 0.001 rounds to km.

[edzer]

Weird, large, complex polygon, and performance -- this may be relevant: http://blog.cleverelephant.ca/2018/09/postgis-external-storage.html

[edzer]

I liked your suggestion of setting precision expressed as distance measures. You can do this now with sf too:

> st_set_precision(shape1, units::set_units(1, mm))
Simple feature collection with 2 features and 2 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 551076.2 ymin: 5190411 xmax: 575786.2 ymax: 5210671
epsg (SRID):    32633
proj4string:    +proj=utm +zone=33 +datum=WGS84 +units=m +no_defs
precision:      1000 
  Focus     A_ha                       geometry
1     1  3257.03 POLYGON ((560246.2 5190411,...
2     2 13661.57 POLYGON ((560356.2 5199901,...

[raff-k]

thank you for the reply and the hints :-)
however, after setting the precision in the new way, it still does not speed up the process (I stopped after 25 minutes).
Is there a way to alter the geom column by "SET STORAGE EXTERNAL" through the sf-package (or maybe through a workaround by restoring the sf in another file format)?

shape1.prec <- st_set_precision(shape1, units::set_units(1, mm))
shape2.prec <- st_set_precision(shape2, units::set_units(1, mm)) # the "problematic" sf

union.prec <- sf::st_union(x = shape1.prec, y = shape2.prec) 

[jannes-m]

If you have already used RSAGA to create one of the shapefiles, why don't you use it for the unioning as well. SAGA is really fast with these kind of operations 😃 .

[MxNl]

Weird, large, complex polygon, and performance -- this may be relevant: http://blog.cleverelephant.ca/2018/09/postgis-external-storage.html

Thanks @edzer for this hint. It helped me a lot. I had to union quite detailed polygons of Europea landmasses (Scandinavian coastline really is an issue here :-)). I stopped sf::st_union() after 4 days. With that trick it was done in less than 30 min.

[LDalby]

@MxNl Just to clarify - you ended up doing your union in PostGIS?

[MxNl]

@MxNl Just to clarify - you ended up doing your union in PostGIS?

@LDalby: That's correct, I used PostGIS with the described trick from R with the DBI package. The polygons that I have for Europe are quite large. So I guess that the poor performance really has to do with the compression of the geometry as described in the link. And again: I stopped sf::st_union after 4 days. PostGIS + trick did it in less than 30 min.