Создание пространственного многоугольника без использования шейп-файла в R

22

Итак, обычный способ, которым мы читаем шейп-файл в R, - это пакет maptools, например:

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

Тем не менее, у меня есть сценарий использования, при котором у меня нет shapefile.shp, но вместо этого у меня есть ряд многоугольных координат

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

и его соответствующая проекция

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

Как мне "создать" sfdata (который будет "объектом многоугольника") непосредственно из этих данных? (без обходного пути создания шейп-файла с этими данными и последующего чтения из вновь созданного шейп-файла)

Кельвин Ченг
источник

Ответы:

35

Сначала получите координаты в матрицу из 2 столбцов:

> 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

Затем создайте Polygon, оберните его в объект Polygons, затем оберните его в объект SpatialPolygons:

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

Причина такого уровня сложности заключается в том, что Polygon - это простое кольцо, объект Polygons может представлять собой несколько колец с идентификатором (в данном случае равным 1) (так же, как один объект в ГИС), а SpatialPolygons может иметь CRS. , Ох, я, наверное, должен установить это:

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

Если вы хотите превратить его в SpatialPolygonsDataFrame (что и происходит с readShapeSpatial, когда шейп-файл является полигоном), выполните:

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

давая это:

> 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
источник
+1 Очень хорошая, понятная экспозиция. Замечательно видеть код разбитым по объяснениям, а не предложенным в виде монолитного блока!
whuber
Отлично ... приятно видеть, как эти объекты соединяются! Нужно видеть больше страниц справки R, написанных ясно, как это.
Симбамангу
Это то, чему я должен заново учить себя каждый раз, когда я хочу это делать, поэтому я пользуюсь любой возможностью, чтобы учить других людей!
Spacedman
1
отлично ... как бы я добавил несколько уникальных многоугольников id (f) во фрейм данных?
Мга,
2
Чтобы этот ответ имел более общую достоверность, не могли бы вы показать, как это сделать в случае нескольких полигонов? Это немного сложно.
Томас
2

Чтобы завершить отличный ответ Spacedman для случая, когда ваши данные будут содержать несколько полигонов, вот код, использующий 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)

Просто для удовольствия, вы можете сравнить с графиком, полученным с ggplot2использованием исходного фрейма данных:

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

Обратите внимание, что в приведенном выше коде предполагается, что у вас есть только один polyogn для каждого идентификатора. Если у некоторых идентификаторов были непересекающиеся многоугольники, я думаю, что нужно добавить еще один столбец в набор данных, сначала group_byсубидид, а затем использоватьgroup_by(upper-id) вместоrowwise

Тот же код, используя purrr::map функции:

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
источник