Creación de un polígono espacial sin usar un archivo de formas en R

Question

Entonces, la forma habitual de leer un shapefile en R es a través del paquete maptools, de esta manera:

sfdata <- readShapeSpatial("/ruta/a/mi/shapefile.shp", proj4string=CRS("+proj=longlat"))

Sin embargo, tengo un caso de uso en el que no tengo un shapefile.shp sino que tengo una serie de coordenadas de polígono

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

y su proyección correspondiente

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

¿Cómo "instancio" sfdata (que será un "objeto de polígono") directamente a partir de estos datos? (sin tener que crear un shapefile con estos datos y luego leerlo desde el shapefile recién creado)

Answer 1

Primero obtén las coordenadas en una matriz de 2 columnas:

> 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

Luego crea un Polígono, envuélvelo en un objeto Polygons, luego envuélvelo en un objeto SpatialPolygons:

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

La razón de este nivel de complejidad es que un Polígono es un anillo simple, un objeto Polygons puede ser varios anillos con un ID (aquí establecido en 1) (así que es como una sola característica en un SIG) y un objeto SpatialPolygons puede tener un SRC. Oh, probablemente debería configurarlo:

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

Si deseas convertirlo en un SpatialPolygonsDataFrame (que es lo que sale de readShapeSpatial cuando el archivo de forma son polígonos) entonces haz:

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

dando esto:

> 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

Answer 2

Para completar la excelente respuesta de Spacedman para el caso en que sus datos contienen múltiples polígonos, aquí hay algo de 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)

Solo por diversión, puedes comparar con el gráfico obtenido con ggplot2 usando el marco de datos inicial:

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

Tenga en cuenta que el código anterior asume que solo tiene un polígono por id. Si algunos ids tuvieran polígonos disjuntos, supongo que se debería agregar otra columna en el conjunto de datos, primero group_by el sub-id, luego usar group_by(upper-id) en lugar de rowwise

Mismo código usando la función 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)}