CODEWITHSUNDEEP
rdata.tabler-sf
user8795501
user8795501

Reputation: 87

st_union (sf) polygons with data.table

In near past I used a library sf with library data.table without special problems. But now I tried one of my code which finished with error. When I merge polygons by id within data.table the results I cannot ploted or write to shp. Workaround is transformation to sp and back to sf but with a lot polygons is this quite time consuming. Small example:

library(data.table)
library(sf)

a=structure(list(Id = c(1L, 2L, 3L, 3L, 1L, 1L, 2L), geometry = structure(list(
  structure(list(structure(c(455311.710673953, 455314.345317508, 
                             455315.716716433, 455316.87150159, 455312.03568616, 455311.710673953, 
                             376691.889842887, 376694.12745275, 376692.828258414, 376689.327648062, 
                             376689.255443473, 376691.889842887), .Dim = c(6L, 2L))), class = c("XY", 
                                                                                                "POLYGON", "sfg")), structure(list(structure(c(455316.87150159, 
                                                                                                                                               455315.716716433, 455319.072917605, 455321.382671023, 455316.87150159, 
                                                                                                                                               376689.327648062, 376692.828258414, 376693.694408316, 376690.518443961, 
                                                                                                                                               376689.327648062), .Dim = c(5L, 2L))), class = c("XY", "POLYGON", 
                                                                                                                                                                                                "sfg")), structure(list(structure(c(455312.93790784, 455316.87150159, 
                                                                                                                                                                                                                                    455317.485088015, 455313.370891238, 455312.93790784, 376686.224010367, 
                                                                                                                                                                                                                                    376685.646434684, 376682.57880773, 376682.542858023, 376686.224010367
                                                                                                                                                                                                ), .Dim = c(5L, 2L))), class = c("XY", "POLYGON", "sfg")), 
  structure(list(structure(c(455318.062480594, 455322.068278933, 
                             455321.815715457, 455314.634074832, 455318.062480594, 376684.52766027, 
                             376682.903819937, 376679.980419059, 376681.027049918, 376684.52766027
  ), .Dim = c(5L, 2L))), class = c("XY", "POLYGON", "sfg")), 
  structure(list(structure(c(455321.346477176, 455323.548076297, 
                             455324.414104129, 455322.104472781, 455321.346477176, 376688.028453727, 
                             376688.569652457, 376686.404430289, 376685.466014762, 376688.028453727
  ), .Dim = c(5L, 2L))), class = c("XY", "POLYGON", "sfg")), 
  structure(list(structure(c(455319.302470828, 455321.743510867, 
                             455323.90891614, 455321.382671023, 455319.302470828, 376693.379039664, 
                             376694.27186193, 376691.456859488, 376690.518443961, 376693.379039664
  ), .Dim = c(5L, 2L))), class = c("XY", "POLYGON", "sfg")), 
  structure(list(structure(c(455315.538127566, 455316.87150159, 
                             455318.105693484, 455318.856486941, 455315.35587659, 455315.538127566, 
                             376689.307811637, 376689.327648062, 376689.653453727, 376687.667430777, 
                             376687.667430777, 376689.307811637), .Dim = c(6L, 2L))), class = c("XY", 
                                                                                                "POLYGON", "sfg"))), n_empty = 0L, crs = structure(list(input = NA_character_, 
                                                                                                                                                        wkt = NA_character_), class = "crs"), class = c("sfc_POLYGON", 
                                                                                                                                                                                                        "sfc"), precision = 0, bbox = structure(c(xmin = 455311.710673953, 
                                                                                                                                                                                                                                                  ymin = 376679.980419059, xmax = 455324.414104129, ymax = 376694.27186193
                                                                                                                                                                                                        ), class = "bbox"))), row.names = c(NA, -7L), class = c("sf", 
                                                                                                                                                                                                                                                                "data.frame"), sf_column = "geometry", agr = structure(c(Id = NA_integer_), class = "factor", .Label = c("constant", 
                                                                                                                                                                                                                                                                                                                                                                         "aggregate", "identity")))


setDT(a)


################original data
#############################

#plot work ok
plot(a$geometry)

##write to shp work ok
st_write(a,"ab.shp", delete_layer=T)


################union polygons by id
####################################
aa<-a[,.(geometry=st_union(geometry)), by=Id] %>% sf::st_as_sf()

##plot error:  "Error in CPL_geos_is_empty(st_geometry(x)) : Not a matrix."
plot(aa$geometry)

## write error: "Error in CPL_write_ogr(obj, dsn, layer, driver, as.character(dataset_options),  :Not a matrix."
st_write(aa,"abc.shp", delete_layer=T)


#####################transform from data.table to sf
####################################################
aa=st_as_sf(aa)

##plot error  "Error in CPL_geos_is_empty(st_geometry(x)) : Not a matrix."
plot(aa)

## write to shp with error: "Error in CPL_write_ogr(obj, dsn, layer, driver, as.character(dataset_options),  :Not a matrix."
st_write(aa,"abc.shp", delete_layer=T)


############################transform sf to sp and back
#######################################################
aa=st_as_sf(as(st_as_sf(aa),"Spatial"))

##plot work ok
plot(aa)

##write to shp work ok
st_write(a,"ab.shp", delete_layer=T)

EDIT:

I found part solution base by https://github.com/Rdatatable/data.table/issues/4681 . If after st_union I add line with code like aa=aa[1:nrow(aa),] ploting and writing working as expected.

setDT(a)
aa<-a[,.(geometry=st_union(geometry)), by=Id] %>% sf::st_as_sf()
aa=aa[1:nrow(aa),]
plot(aa)
st_write(aa,"aa.shp")

Upvotes: 2

Views: 3011

Answers (1)

at80
at80

Reputation: 780

It's because your id 1 has no overlap and consists of multiple polygons, so it's coercing it to a multipolygon instead of a polygon. The others only have two overlapping polygons, so its dissolving them to a single polygon. I tested this using the dplyr method and received a similar error. You can see this in the data.table below.

aa<- a[,.(geometry=st_union(geometry)), by=Id]
#   Id geometry
#1:  1  <XY[3]> <--- three polygons
#2:  2  <XY[1]>
#3:  3  <XY[1]>

# Geometry set for 3 features 
# geometry type:  MULTIPOLYGON <--- Multipolygon conflicting with polygon geometries
# dimension:      XY
# bbox:           xmin: 455311.7 ymin: 376685.5 xmax: 455324.4 ymax: 376694.3
# CRS:            NA
# MULTIPOLYGON (((455321.3 376688, 455323.5 37668...
# POLYGON ((455315.5 376689.3, 455316.9 376689.3,...
# POLYGON ((455318.1 376684.5, 455322.1 376682.9,...

This may be a bug. An alternative would be to use st_combine instead of st_union if that works for your needs. Or you could try getting the bad geometries back to the correct geometry.

Edit: Looking at it more closely you probably don't want just the st_combine method. It seems that st_combine will combine but not union the geometries. One work around was to use st_combine then st_union. This works because it combines all of the features by group, then it unions them by feature.

# first way with st_combine
# EDIT this didn't work with just st_combine added the st_union at the end
aa<-a[,.(geometry=st_combine(geometry)), by=Id] %>% sf::st_as_sf() %>% st_union(by_feature=TRUE)
plot(aa$geometry)

# The hackier way
aa<- a[,.(geometry=st_union(geometry)), by=Id] %>% sf::st_as_sf() 
newgeom <- aa[which(st_geometry_type(aa) == 'POLYGON'), ] %>% st_cast("MULTIPOLYGON")
fixedaa <- rbind(aa[which(st_geometry_type(aa) != 'POLYGON'), ], newgeom)

Upvotes: 1

Related Questions