Memperbaiki lubang yatim di R

Saya mencoba melakukan penyatuan pada bidang yang sama setelah menggabungkan dua bentuk yang berdekatan. Shapefile berakhir dengan setidaknya satu irisan tipis ruang di antara mereka. Ketika saya mencoba serikat pekerja saya mendapatkan kesalahan lubang yatim berikut:

Kesalahan pada createPolygonsComment (p): rgeos_PolyCreateComment: lubang yatim, tidak dapat menemukan mengandung poligon untuk lubang di indeks 17

Saya telah mengunggah contoh yang dapat direproduksi ke Dropbox di tautan ini .

Berikut adalah kode untuk membuat ulang masalah:

#loading required packages
require(sp)    
require(rgdal)
require(maptools)
require(rgeos)

#load example data, set "dsn=" to your working directory or specify the path
example <- readOGR(dsn=".",layer="ReproducibleExample")

#Attempting a UnionSpatialPolygons based on the COUNTY field
example.df <- as(example, "data.frame")
countycol <- example.df$COUNTY
example.diss <- unionSpatialPolygons(example, countycol)

Pengembalian:

Kesalahan pada createPolygonsComment (p): rgeos_PolyCreateComment: lubang yatim, tidak dapat menemukan mengandung poligon untuk lubang di indeks 17

Mencoba perbaikan yang diusulkan di sini dan di sini :

slot(example, "polygons") <- lapply(slot(example, "polygons"), checkPolygonsHoles)

Ini mengembalikan kesalahan yang sama yang berasal dari upaya serikat tetapi dengan nomor indeks berbeda:

rgeos_PolyCreateComment: lubang yatim, tidak dapat menemukan mengandung poligon untuk lubang pada indeks 30

Mencoba perbaikan yang diusulkan dalam tutorial bermanfaat Roger Bivand

fix <- slot(example, "polygons")
fixa <- lapply(fix, checkPolygonsHoles)

Mengembalikan kesalahan yang sama pada indeks 30 seperti di atas.

Yang lain telah mengangkat masalah ini di sini dan di sini , dan sementara solusi yang diajukan di atas tampaknya bekerja untuk beberapa kasus, kasus lain tidak terselesaikan. Satu pengguna menggunakan QGIS untuk mengatasi masalah, dan yang lain memiliki 2 dari 3 item yang diperbaiki, tetapi tidak ada resolusi untuk yang terakhir.

Tampaknya orang terus mengalami masalah meskipun kode ini berfungsi dari waktu ke waktu. Adakah yang menemukan solusi dalam R?

Saya telah melakukan "perbaikan geometri" alat di ArcGIS, dan itu memperbaiki masalah, tetapi sepertinya harus ada perbaikan di R.

Saya telah menganalisis masalah geometri pada data terlampir, dan sepertinya tidak HANYA dimiliki orphaned holestetapi juga geometry validity issues. Memang benar bahwa orphaned holeentah bagaimana masalah validitas geometri, tetapi rgeos tidak menanganinya dengan cara yang sama, seperti untuk lubang yatim, kesalahan muncul, bukannya peringatan sederhana. Seperti yang Anda tunjukkan, itu adalah petunjuk untuk memeriksa lubang poligon, tetapi itu tidak selalu berhasil ketika diterapkan untuk memperbaiki lubang yatim.

Jadi ayo:

1. bersihkan data Anda (yang diperlukan jika Anda ingin melakukan geoproses seperti gabungan)

2. gunakan data yang dibersihkan dengan proses penyatuan Anda

1. Membersihkan geometri Memperbaiki geometri dalam R terkadang sulit, jadi saya telah mencoba membuat paket R eksperimental (lihat https://github.com/eblondel/cleangeo ) yang bermaksud memfasilitasi pembersihan spobjek (saat ini terbatas pada bentuk poligonal). Anda dapat menginstal paket dengan:

require(devtools)
install_github("eblondel/cleangeo")
require(cleangeo)

Untuk memulai, ada baiknya Anda melihat apa saja masalah geometri dengan data sumber Anda. Untuk ini, Anda dapat menjalankan yang berikut ini (data Anda besar sehingga dapat memakan waktu):

#get a report of geometry validity & issues for a sp spatial object
report <- clgeo_CollectionReport(sp)
summary <- clgeo_SummaryReport(report)
issues <- report[report$valid == FALSE,]

Dengan ini, Anda akan melihat bahwa data Anda memiliki 2 jenis masalah: orphaned holesdan geometry validity issues. Baik (dan tidak hanya lubang yatim) cenderung membuat unionproses gagal, sehingga data harus dibersihkan sebelumnya, dengan cara otomatis jika memungkinkan. Untuk reproduksi cepat, kode sampel pertama di bawah ini hanya mengambil subset fitur yang ditandai sebagai mencurigakan (kecuali yang terbaru, dengan indeks = 9002 dalam data asli - lihat catatan saya di bawah ini tentang ini)

#get suspicious features (indexes)
nv <- clgeo_SuspiciousFeatures(report)
mysp <- sp[nv[-14],]

#try to clean data
mysp.clean <- clgeo_Clean(mysp, print.log = TRUE)

#check if they are still errors
report.clean <- clgeo_CollectionReport(mysp.clean)
summary.clean <- clgeo_SummaryReport(report.clean)

Jika clgeo_Cleanberhasil dengan baik, Anda harus mendapatkan semua geometri yang valid sekarang. Anda dapat menerapkan ini pada dataset lengkap (kecuali indeks fitur = 9002)

#try to clean data
mysp <- sp[-9002,]
mysp.clean <- clgeo_Clean(mysp, print.log = TRUE)

#check if they are still errors
report.clean <- clgeo_CollectionReport(mysp.clean)
summary.clean <- clgeo_SummaryReport(report.clean)

2. Proses penyatuan Sekarang, mari kita lihat apakah unionbekerja pada dataset ini:

#Attempting a UnionSpatialPolygons based on the COUNTY field
mysp.df <- as(mysp, "data.frame")
countycol <- mysp.df$COUNTY
mysp.diss <- unionSpatialPolygons(mysp.clean, countycol)

Catatan: seperti yang dikatakan sebelumnya, saya telah menghapus satu fitur (indeks = 9002). Dengan memplotnya:, plot(sp[9002,])Anda akan melihat bahwa fitur ini sangat (sangat) kompleks. Saya telah mengecualikannya dari sampel hanya karena memeriksa lubang terlalu banyak memakan waktu. Mari kita lihat sekarang jika masalah yang sama terjadi menggunakan readShapePoly(dari maptools) untuk membaca data ...

3. Beralih ke readShapePoly vs. readOGR untuk membaca data (PEMBARUAN)

readOGRbukan satu-satunya fungsi yang tersedia untuk membaca shapefile. Anda juga dapat menggunakan readShapePolydari maptoolspaket, umumnya lebih berkinerja daripada yang pertama:

require(maptools)
mysp <- readShapePoly("ReproducibleExample.shp")

Selain berjalan lebih cepat:

• jika Anda menggunakan kode di atas berdasarkan clgeo_CollectionReport, tidak ada masalah lubang yatim, tetapi masih masalah geometri.

• Membersihkan geometri dengan clgeo_Cleanjuga berjalan dengan baik, dan sekarang tidak terjebak dengan indeks fitur 9002

• Dan ... proses persatuan bekerja.

Lihat di bawah hasil plot:

#plot the result
plot(mysp, border= "lightgray")
plot(mysp.diss, border="red", add = TRUE)

Kesimpulan : lebih suka maptools untuk membaca data shapefile Anda, dan pertimbangkan untuk menggunakan cleangeo untuk membersihkan data Anda sebelum melakukan geoproses.

Solusi mudah yang terus bekerja untuk saya di R adalah menerapkan buffer nol lebar :

#loading required packages
require(sp)    
require(rgdal)
require(maptools)
require(rgeos)

#load example data, set "dsn=" to your working directory or specify the path
example <- readOGR(dsn=".",layer="ReproducibleExample")

#project your data (I'm using California Albers here) and apply a zero-width buffer
example <- spTransform(example, CRS("+init=epsg:3310"))
example <- gBuffer(example, byid = T, width = 0)

#Attempting a UnionSpatialPolygons based on the COUNTY field
example.df <- as(example, "data.frame")
countycol <- example.df$COUNTY
example.diss <- unionSpatialPolygons(example, countycol)

unionSpatialPolygons membutuhkan beberapa saat dengan kumpulan data ini, tetapi tampaknya berfungsi dengan baik.