Instanciando polígono espacial sem usar um shapefile em R

22

Portanto, a maneira usual de ler um shapefile no R é através do pacote maptools, assim:

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

No entanto, tenho um caso de uso em que não tenho um shapefile.shp, mas tenho uma série de coordenadas poligonais

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

e sua projeção correspondente

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

Como instanciar sfdata (que será um "objeto de polígono") diretamente desses dados? (sem seguir uma maneira indireta de criar um shapefile com esses dados e depois ler o shapefile recém-criado)

Calvin Cheng
fonte

Respostas:

35

Primeiro, obtenha as coordenadas em uma matriz de 2 colunas:

> 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

Em seguida, crie um polígono, envolva-o em um objeto Polígonos e, em seguida, em um objeto SpatialPolygons:

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

A razão para esse nível de complexidade é que um polígono é um anel simples, um objeto de polígonos pode ter vários anéis com um ID (aqui definido como 1) (é como um único recurso em um GIS) e um SpatialPolygons pode ter um CRS . Ooh, eu provavelmente deveria configurá-lo:

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

Se você deseja transformá-lo em um SpatialPolygonsDataFrame (que é o que vem de readShapeSpatial quando o shapefile é polígono), faça:

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

dando isso:

> 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
fonte
+1 Exposição muito agradável e clara. É ótimo ver o código dividido por explicações, em vez de ser oferecido como um bloco monolítico!
whuber
Excelente ... ótimo ver como esses objetos são colocados juntos! Precisa ver mais das páginas de ajuda do R escritas claramente dessa forma.
Simbamangu
É algo que eu tenho que me re-ensinar toda vez que quero fazer isso, então aproveito a oportunidade para ensinar outras pessoas!
Spacedman
1
excelente ... como eu adicionaria vários polígonos de identificação (f) exclusivos ao quadro de dados?
mga
2
Para que esta resposta tenha validade mais geral, você poderia mostrar como fazer isso no caso de vários polígonos? Isso é um pouco complicado.
Tomas
2

Para concluir a excelente resposta de Spacedman para o caso em que seus dados conteriam vários polígonos, aqui está um código usando 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)

Apenas por diversão, você pode comparar com o gráfico obtido com o ggplot2uso do quadro de dados inicial:

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

Observe que o código acima pressupõe que você tenha apenas um polognio por ID. Se alguns IDs tiverem polígonos separados, acho que se deve adicionar outra coluna no conjunto de dados, primeiro group_byo sub-ID e, em seguida, use em group_by(upper-id)vez derowwise

O mesmo código usando a purrr::mapfunção:

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
fonte