Odwrotne obcinanie (kasowanie) w R?

14

Odwrotny klip zapisuje tylko część obiektu przestrzennego, która znajduje się poza granicami innego obiektu, w przeciwieństwie do zwykłego klipu, który zapisuje części znajdujące się wewnątrz innego obiektu.

Wykonujesz odwrotny klip w ArcMap? pokazuje, jak to zrobić w ArcMap.

Jak to zrobić w R?

Powtarzalny przykład (na komputerach z systemem Linux):

system("wget 'https://github.com/Robinlovelace/Creating-maps-in-R/archive/master.zip' -P /tmp/")
unzip("/tmp/master.zip", exdir = "/tmp/master")
uk <- readOGR("/tmp/master/Creating-maps-in-R-master/data/", "ukbord")
lnd <- readOGR("/tmp/master/Creating-maps-in-R-master/data/", "LondonBoroughs")
plot(uk)
plot(lnd, add = T, col = "black")

Chcę tutaj uratować całą Wielką Brytanię oprócz Londynu. Wizualnie chcę, aby czarny kształt na wynikowym obrazie był dziurą.

wprowadź opis zdjęcia tutaj

RobinLovelace
źródło

Odpowiedzi:

4

Odpowiedź na proste funkcje:

Pakiet sf korzysta z Geometry Engine Open Source i jako taki może uzyskać dostęp do listy poleceń, takich jak st_within itp.

Jedno takie polecenie, st_difference, wykona zadanie:

require(sf)

# make a square simple feature
s <- rbind(c(1,1), c(1,5), c(5,5), c(5,1), c(1,1))
s.sf <-st_sfc(st_polygon(list(s)))
s.pol = st_sf(ID = "sq", s.sf)

# make a smaller triangle simple feature
s2 <- rbind(c(2,2), c(3,4), c(4,2), c(2,2))
s2.sf <-st_sfc(st_polygon(list(s2)))
s2.pol = st_sf(ID = "tr", s2.sf)

# find the 'difference', i.e. reverse of st_intersection
t <- st_difference(s.pol,s2.pol)

plot(t)

# have a look at the new geometry, a well known text format with exterior followed by hole
st_geometry(t)[[1]]
POLYGON((1 1, 1 5, 5 5, 5 1, 1 1), (2 2, 4 2, 3 4, 2 2))

zobacz także na dole tego artykułu

można to również zrobić, zmuszając Sp do sf za pomocą st_as_sf. Zwróć uwagę na ostrzeżenia, ponieważ zarządzanie atrybutami może być trudne!

Sam
źródło
12

Wydaje się, że jest to prosta aplikacja gDifferencez rgeospakietu:

> require(rgeos)
> ukhole = gDifference(uk, lnd)
Warning message:
In RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, "rgeos_difference") :
  spgeom1 and spgeom2 have different proj4 strings
> plot(ukhole)

Ostrzeżenie o rzutowaniu jest spowodowane tym LondonBoroughs, że .prjplik kształtu nie ma pliku.

Aby upewnić się, że jest to dziura, a nie kontur lub inny solidny wielokąt:

> gArea(lnd) + gArea(ukhole) - gArea(uk)
[1] 0
Spacedman
źródło
Bardzo proste, dzięki za szybką odpowiedź. Byłbym zainteresowany spojrzeniem na kod źródłowy tych funkcji, aby zobaczyć, co się dzieje pod maską.
RobinLovelace
Głęboko nazywają po prostu GEOS, który jest biblioteką kodów C funkcji geometrii trac.osgeo.org/geos
Spacedman,
Ciekawe - i pomaga wyjaśnić, dlaczego jest dość szybki. Na podstawie tej strony wydaje się, że nie jest ona aktywnie rozwijana. Czy ktoś może to potwierdzić / obalić? svn.osgeo.org/geos/branches/3.4/ChangeLog
RobinLovelace
1
Na pewno jest rozwinięty. Zobacz oś czasu trac.osgeo.org/geos/timeline lub archiwa list mailingowych lists.osgeo.org/pipermail/geos-devel
user30184
5

Trochę późno na imprezę, ale istnieje prosty sposób, aby to zrobić za pomocą maski, używając argumentu „odwrotnego”;

ukhole <- mask(uk, lnd, inverse = TRUE)
Coding4Biology
źródło
Z pakietu rastrowego. A jeśli masz jakieś pomysły?
RobinLovelace,