Czy przyciąć wielokąt i zachować dane?

13

Mam te dwa wielokąty:

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)

wprowadź opis zdjęcia tutaj

Chcę wyciąć części, spdf1które są przecięte spdf2. Chcę spdf1jednak pozostać jako SpatialPolygonsDataFrame i zachować wszelkie zawarte w nim informacje spdf1@data.

Próbowałem gDifference w następujący sposób, który wycina części, spdf1które są przecinane spdf2, ale następnie przekształca spdf1się w SpatialPolygons (tj. Odrzuca zawarte w nim informacje spdf1@data).

gDifference(spdf1, spdf2, byid=T)

Jak mogę pokroić spdf1w spdf2ale zachowują dane zawarte w spdf1@data? Sprawdziłem i wypróbowałem to podobne pytanie bez nakładania wielokąta na SpatialPointsDataFrame i zachowania danych SPDF?

luciano
źródło

Odpowiedzi:

9

Najprostsze podejście byłoby

  library(raster)
  x <- spdf1 - spdf2

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

Zobacz? 'Pakiet raster' (sekcja XIV), aby uzyskać więcej funkcji dotyczących nakładki wielokątów. Funkcje te wykorzystują podstawowe funkcje rgeo pod maską, w funkcjach „na poziomie użytkownika” (w przeciwieństwie do „poziomu programisty”).

Robert Hijmans
źródło
Co rozumiesz przez „na poziomie użytkownika” (w przeciwieństwie do funkcji „na poziomie programisty”)?
luciano,
rgeoszapewnia operacje geometryczne, ale nie zajmuje się atrybutami danych. Dlatego korzystanie z tych funkcji wymaga dużo pracy, aby utrzymać wszystko razem. Funkcje rastrowe upraszczają to i zachowują się jak podobne funkcje w oprogramowaniu GIS,
Robert Hijmans,
Tak, ale może to prowadzić do błędu w SpatialPolygonsDataFrame (part2, x @ data [match (row.names (part2),: row.names danych i ID wielokątów nie pasują
jebyrnes)
To byłby błąd. Czy możesz podać przykład?
Robert Hijmans
4

Obejściem tego problemu byłoby ponowne dodanie atrybutów po wykonaniu klipu podczas konwersji z SpatialPolygonsna 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")

wprowadź opis zdjęcia tutaj

Andre Silva
źródło
Wygląda na to, że mam problem, tylko klip wyjściowy w mojej konkretnej instancji tworzy nazwy rownów ze spdf1, które nie istnieją (jako prosty gsub, aby pozbyć się drugiej cyfry z rzędu. Nazwy powinny zadbać o dopasowanie nazw rownów, nie? )
jebyrnes,
Błąd SpatialPolygonsDataFrame (zacisk dane = as.data.frame (all_spdfs_together @ dane)): długość niedopasowanie obiekt ma zacisk 1718 wielokątów przedmioty, ale as.data.frame (all_spdfs_together @ danych) ma 86 rzędów
jebyrnes
Jasne - mam mnóstwo wieloboków w oceanie. Niektóre są nieprawidłowo umieszczone na lądzie lub pokrywają się z ziemią. Chcę to wyciąć, więc mam tylko obszary znajdujące się w oceanie. Mam plik kształtu linii brzegowej do porównania. Oto kod - gist.github.com/jebyrnes/c2e8d2b6c82849dad3a813d952ab8bb0
jebyrnes
1
Nevermind - prace rozwiązanie raster :: erase (nie miał wcześniej do jakiegoś dziwnego powodu)
jebyrnes