Membalikkan kliping (menghapus) di R?

14

Klip terbalik hanya menyimpan bagian dari objek spasial Anda yang berada di luar batas objek lain, tidak seperti klip biasa yang menyimpan bagian-bagian yang ada di dalam objek lain.

Melakukan klip terbalik di ArcMap? menunjukkan cara melakukannya di ArcMap.

Bagaimana saya melakukan ini di R?

Contoh yang dapat direproduksi (pada mesin 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")

Yang ingin saya lakukan di sini adalah menyelamatkan seluruh Inggris kecuali London. Secara visual, saya ingin bentuk hitam pada gambar yang dihasilkan menjadi sebuah lubang.

masukkan deskripsi gambar di sini

RobinLovelace
sumber

Jawaban:

4

Jawaban untuk Fitur Sederhana:

Paket sf mengacu pada Open Source Engine Geometry, dan dengan demikian dapat mengakses daftar perintah seperti st_within dll.

Salah satu perintah tersebut, st_difference, akan melakukan pekerjaan:

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))

lihat juga di bagian bawah artikel ini

juga bisa dilakukan dengan memaksa Sp ke sf dengan st_as_sf. Mengindahkan peringatan karena atribut bisa rumit untuk dikelola!

Sam
sumber
12

Tampaknya aplikasi sederhana gDifferencedari rgeospaket:

> 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)

Peringatan proyeksi adalah karena LondonBoroughsshapefile tidak memiliki .prjfile.

Hanya untuk memastikan bahwa itu adalah lubang dan bukan garis besar atau poligon padat lainnya:

> gArea(lnd) + gArea(ukhole) - gArea(uk)
[1] 0
Spacedman
sumber
Sooo sederhana, terima kasih atas tanggapannya yang cepat. Akan tertarik melihat kode sumber fungsi-fungsi ini untuk melihat apa yang terjadi di bawah tenda.
RobinLovelace
Jauh di lubuk hati mereka hanya memanggil GEOS yang merupakan pustaka kode-C fungsi geometri trac.osgeo.org/geos
Spacedman
Menarik - dan membantu menjelaskan mengapa itu cukup cepat, saya kira. Berdasarkan halaman ini, tampaknya tidak sedang dikembangkan secara aktif, adakah yang bisa mengkonfirmasi / membantahnya? svn.osgeo.org/geos/branches/3.4/ChangeLog
RobinLovelace
1
Pasti itu dikembangkan. Lihat timeline trac.osgeo.org/geos/timeline atau arsip daftar mailing lists.osgeo.org/pipermail/geos-devel
user30184
5

Agak terlambat ke pesta, tetapi ada cara sederhana untuk melakukan ini dengan topeng menggunakan argumen 'terbalik';

ukhole <- mask(uk, lnd, inverse = TRUE)
Coding4Biology
sumber
Dari paket raster. Dan dari sf ada ide?
RobinLovelace