Instantiating poligon spasial tanpa menggunakan shapefile di R

22

Jadi, cara biasa kita membaca shapefile di R adalah melalui paket maptools, seperti ini:

sfdata <- readShapeSpatial("/path/to/my/shapefile.shp", proj4string=CRS("+proj=longlat"))

Namun, saya memiliki kasus penggunaan di mana saya tidak memiliki shapefile.shp tetapi saya memiliki serangkaian koordinat poligon

16.484375 59.736328125,17.4951171875 55.1220703125,24.74609375 55.0341796875,22.5927734375 61.142578125,16.484375 59.736328125

dan proyeksi yang sesuai

coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 

Bagaimana cara "instantiate" sfdata (yang akan menjadi "objek poligon") langsung dari data ini? (tanpa berputar-putar cara membuat shapefile dengan data ini dan kemudian membaca dari shapefile yang baru dibuat)

Calvin Cheng
sumber

Jawaban:

35

Pertama-tama dapatkan koordinat ke dalam matriks 2-kolom:

> xym
         [,1]     [,2]
[1,] 16.48438 59.73633
[2,] 17.49512 55.12207
[3,] 24.74609 55.03418
[4,] 22.59277 61.14258
[5,] 16.48438 59.73633

Kemudian buat Polygon, bungkus itu menjadi objek Polygons, kemudian bungkus itu menjadi objek SpatialPolygons:

> library(sp)
> p = Polygon(xym)
> ps = Polygons(list(p),1)
> sps = SpatialPolygons(list(ps))

Alasan tingkat kerumitan ini adalah bahwa Poligon adalah cincin sederhana, objek Poligon dapat berupa beberapa cincin dengan ID (di sini diatur ke 1) (jadi seperti fitur tunggal dalam GIS) dan SpatialPoligon dapat memiliki CRS . Ooh, saya mungkin harus mengaturnya:

> proj4string(sps) = CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")

Jika Anda ingin mengubahnya menjadi SpatialPolygonsDataFrame (yang merupakan milik readShapeSpatial kami ketika shapefile adalah poligon) maka lakukan:

> data = data.frame(f=99.9)
> spdf = SpatialPolygonsDataFrame(sps,data)
> spdf

memberikan ini:

> summary(spdf)
Object of class SpatialPolygonsDataFrame
Coordinates:
       min      max
x 16.48438 24.74609
y 55.03418 61.14258
Is projected: FALSE 
proj4string :
[+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0]
Data attributes:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   99.9    99.9    99.9    99.9    99.9    99.9 
Spacedman
sumber
+1 Eksposisi yang sangat bagus dan jelas. Sangat menyenangkan melihat kode dipecah oleh penjelasan daripada ditawarkan sebagai blok monolitik!
whuber
Luar biasa ... bagus untuk melihat bagaimana benda-benda ini disatukan! Perlu melihat lebih banyak halaman bantuan R yang ditulis dengan jelas seperti ini.
Simbamangu
Ini adalah sesuatu yang saya harus ajarkan kembali pada diri saya sendiri setiap kali saya ingin melakukannya, jadi saya mengambil kesempatan untuk mengajar orang lain!
Spacedman
1
unggul ... bagaimana saya bisa menambahkan beberapa id unik (f) poligon ke bingkai data?
mga
2
Agar jawaban ini memiliki validitas yang lebih umum, dapatkah Anda menunjukkan bagaimana melakukannya dalam beberapa poligon? Ini agak sulit.
Tomas
2

Untuk melengkapi jawaban Spacedman yang luar biasa untuk kasus di mana data Anda akan mengandung banyak poligon, berikut adalah beberapa kode menggunakan dplyr:

library(dplyr)
library(ggplot2)
library(sp)
## use data from ggplot2:::geom_polygon example:
positions <- data.frame(id = rep(factor(c("1.1", "2.1", "1.2", "2.2", "1.3", "2.3")), each = 4),
                    x = c(2, 1, 1.1, 2.2, 1, 0, 0.3, 1.1, 2.2, 1.1, 1.2, 2.5, 1.1, 0.3,
                          0.5, 1.2, 2.5, 1.2, 1.3, 2.7, 1.2, 0.5, 0.6, 1.3),
                    y = c(-0.5, 0, 1, 0.5, 0, 0.5, 1.5, 1, 0.5, 1, 2.1, 1.7, 1, 1.5,
                          2.2, 2.1, 1.7, 2.1, 3.2, 2.8, 2.1, 2.2, 3.3, 3.2)) %>% as.tbl


df_to_spp <- positions %>%
  group_by(id) %>%
  do(poly=select(., x, y) %>%Polygon()) %>%
  rowwise() %>%
  do(polys=Polygons(list(.$poly),.$id)) %>%
  {SpatialPolygons(.$polys)}

## plot it
plot(df_to_spp)

Hanya untuk bersenang-senang, Anda dapat membandingkan dengan plot yang diperoleh dengan ggplot2menggunakan frame data awal:

ggplot(positions) + 
  geom_polygon(aes(x=x, y=y, group=id), colour="black", fill=NA)

Perhatikan bahwa kode di atas mengasumsikan bahwa Anda hanya memiliki satu polyogn per id. Jika beberapa id telah memisahkan poligon, saya kira kita harus menambahkan kolom lain dalam dataset, pertama group_by- tama sub-id, kemudian gunakan group_by(upper-id)alih-alihrowwise

Kode yang sama menggunakan purrr::mapfungsi:

df_to_spp <- positions %>%
  nest(-id) %>%
  mutate(Poly=purrr::map(data, ~select(., x, y)  %>% Polygon()),
         polys=map2(Poly, id, ~Polygons(list(.x),.y))) %>%
  {SpatialPolygons(.$polys)}
Matifou
sumber