Обрезать полигон и сохранить данные?

13

У меня есть эти два полигона:

library(sp); library(rgeos); library(maptools)

coords1 <- matrix(c(-1.841960, -1.823464, -1.838623, -1.841960, 55.663696,
                    55.659178, 55.650841, 55.663696), ncol=2)
coords2 <- matrix(c(-1.822606, -1.816790, -1.832712, -1.822606, 55.657887,
                    55.646806, 55.650679, 55.657887), ncol=2)
p1 <- Polygon(coords1)
p2 <- Polygon(coords2)
p1 <- Polygons(list(p1), ID = "p1")
p2 <- Polygons(list(p2), ID = "p2")
myPolys <- SpatialPolygons(list(p1, p2))
spdf1 = SpatialPolygonsDataFrame(myPolys, data.frame(variable1 = c(232,
                                                                   242), variable2 = c(235, 464), row.names = c("p1", "p2")))
proj4string(spdf1) <- CRS("+proj=longlat +datum=WGS84 +ellps=WGS84
+towgs84=0,0,0")
plot(spdf1, col="red")

coords1a <- matrix(c(-1.830219, -1.833753, -1.821154, -1.830219, 55.647353,
                     55.656629, 55.652122, 55.647353), ncol=2)
p1a <- Polygon(coords1a)
p1a <- Polygons(list(p1a), ID = "p1a")
myPolys1 <- SpatialPolygons(list(p1a))
spdf2 = SpatialPolygonsDataFrame(myPolys1, data.frame(variable1 = c(2),
                                                      variable2 = c(3), row.names = c("p1a")))
proj4string(spdf2) <- CRS("+proj=longlat +datum=WGS84 +ellps=WGS84
+towgs84=0,0,0")
plot(spdf2, col="yellow", add=T)

введите описание изображения здесь

Я хочу вырезать части spdf1, которые пересекаются spdf2. Тем не менее, я хочу spdf1остаться как SpatialPolygonsDataFrame и сохранить любую информацию, содержащуюся в spdf1@data.

Я попробовал gDifference следующим образом: он вырезает части spdf1, которые пересекаются spdf2, но затем конвертируется spdf1в SpatialPolygons (т.е. отбрасывает информацию, содержащуюся в spdf1@data).

gDifference(spdf1, spdf2, byid=T)

Как я могу нарезают spdf1с , spdf2но сохраняют содержащиеся данные в spdf1@data? Я проверил и попробовал этот похожий вопрос без того, как наложить многоугольник на SpatialPointsDataFrame и сохранить данные SPDF?

Лучиано
источник

Ответы:

9

Самый простой подход будет

  library(raster)
  x <- spdf1 - spdf2

  # or, more formally
  y <- erase(spdf1,  spdf2)

Смотрите? 'Raster-package' (раздел XIV) для большего количества функций, которые имеют дело с наложением полигонов. Эти функции используют базовые функции rgeos под капотом, в функциях «пользовательского уровня» (в отличие от «уровня разработчика»).

Роберт Хейманс
источник
Что вы подразумеваете под «в функциях уровня пользователя (в отличие от функций уровня разработчика)»?
Лучано
rgeosобеспечивает геометрические операции, но не имеет дело с атрибутами данных. Таким образом, использование этих функций требует много работы, чтобы все было вместе. Растровые функции упрощают это и ведут себя как аналогичные функции в программном обеспечении ГИС,
Роберт Хейманс
Да, но это может привести к ошибке в SpatialPolygonsDataFrame (part2, x @ data [match (row.names (part2),: row.names данных и идентификаторы полигонов не совпадают)
jebyrnes
Это было бы ошибкой. Можете ли вы привести пример этого?
Роберт Хейманс
4

Обходной путь может состоять в том, чтобы повторно добавить атрибуты после выполнения клипа при преобразовании из SpatialPolygonsв SpatialPolygonsDataFrame.

sp3 <- gDifference(spdf1, spdf2, byid = TRUE)
row.names(sp3) <- row.names(spdf1)

spdf3 <- SpatialPolygonsDataFrame(sp3, data = spdf1@data)

spdf3@data

   variable1 variable2
p1       232       235
p2       242       464

plot(spdf3, col="red")

введите описание изображения здесь

Андре Силва
источник
Это похоже на проблему, которую я имею, только выходной клип в моем конкретном экземпляре создает имена строк из spdf1, которые не существуют (как простой gsub, чтобы избавиться от второй цифры в строке. Имена должны заботиться о совпадении имен строк, нет? )
Jebyrnes
Ошибка в SpatialPolygonsDataFrame (clip, data = as.data.frame (all_spdfs_together @ data)): Несоответствие длины объекта: клип содержит 1718 объектов Polygons, но as.data.frame (all_spdfs_together @ data) имеет 86 строк
jebyrnes
Конечно, у меня есть куча полигонов в океане. Некоторые из них неправильно размещены на земле или перекрываются с землей. Я хочу вырезать это, чтобы у меня были только области, которые находятся в океане. У меня есть шейп-файл береговой линии для сравнения. Вот код - gist.github.com/jebyrnes/c2e8d2b6c82849dad3a813d952ab8bb0
jebyrnes
1
неважно - решение raster :: erase работает (раньше по какой-то странной причине его не было)
jebyrnes