#r #data.table #sf
#r #данные.таблица #sf
Вопрос:
В недалеком прошлом я без особых проблем использовал библиотеку sf с библиотекой data.table. Но теперь я попробовал один из своих кодов, который завершился с ошибкой. Когда я объединяю полигоны по идентификатору в data.table, результаты, которые я не могу отобразить или записать в shp. Обходным путем является преобразование в sp и обратно в sf, но с большим количеством полигонов это довольно трудоемко. Небольшой пример:
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)
Редактировать:
Я нашел часть базы решений с помощью https://github.com/Rdatatable/data.table/issues/4681 . Если после st_union я добавлю строку с кодом типа aa=aa[1:nrow (aa),] построение и запись работают так, как ожидалось.
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")
Комментарии:
1. Привет! Я думаю, вам может повезти больше, если вы создадите новую проблему в репозитории sf и задокументируете проблему (возможно, с меньшим примером).
Ответ №1:
Это потому, что ваш идентификатор 1 не имеет перекрытия и состоит из нескольких полигонов, поэтому он приводит его к мультиполигону вместо полигона. У остальных есть только два перекрывающихся многоугольника, поэтому они сводятся к одному полигону. Я протестировал это с помощью метода dplyr и получил аналогичную ошибку. Вы можете увидеть это в data.table ниже.
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,...
Это может быть ошибкой. Альтернативой может быть использование st_combine вместо st_union , если это подходит для ваших нужд. Или вы могли бы попытаться вернуть неправильные геометрии к правильной геометрии.
Редактировать: если посмотреть на это более внимательно, вам, вероятно, нужен не только st_combine метод. Кажется, что st_combine это объединит, но не объединит геометрию. Одним из способов обхода было использовать st_combine then st_union . Это работает, потому что объединяет все объекты по группам, а затем объединяет их по объектам.
# 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)
Комментарии:
1. Привет, спасибо за ответ. Я нашел часть базы решений с помощью github.com/Rdatatable/data.table/issues/4681 . Если после st_union я добавлю строку с кодом типа aa=aa[1:nrow (aa),] построение графика и запись работают так, как ожидалось.