Crea una instancia de polígono espacial sin usar un shapefile en R

19

Por lo tanto, la forma habitual en que leemos un shapefile en R es a través del paquete maptools, como este:

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

Sin embargo, tengo un caso de uso por el que no tengo un shapefile.shp pero en cambio tengo una serie de coordenadas de polígonos

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 puedo "instanciar" sfdata (que será un "objeto polígono") directamente desde estos datos? (sin ir de una manera indirecta de crear un shapefile con estos datos y luego leer el nuevo shapefile creado)

    
pregunta Calvin Cheng 30.12.2011 - 05:02

2 respuestas

30

Primero obtenga 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 Polígonos, luego envuélvelo en un objeto SpatialPolygons:

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

El motivo de este nivel de complejidad es que un Polígono es un anillo simple, un objeto Polígonos puede ser varios anillos con una ID (aquí establecido en 1) (por lo tanto, es como una única característica en un SIG) y un Polígono espacial puede tener un CRS. Ooh, probablemente debería configurarlo:

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

Si desea convertirlo en un SpatialPolygonsDataFrame (que es el resultado de readShapeSpatial cuando el shapefile es polígonos), haga lo siguiente:

> 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 
    
respondido por el Spacedman 30.12.2011 - 11:16
2

Para completar la excelente respuesta de Spacedman en el caso de que sus datos contengan múltiples polígonos, aquí hay un código que usa 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, puede 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 polyogn por id. Si algunos identificadores tenían polígonos separados, creo que uno debería agregar otra columna en el conjunto de datos, primero group_by el sub-id, luego use group_by(upper-id) en lugar de rowwise

El mismo código que usa 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)}
    
respondido por el Matifou 08.10.2016 - 22:04

Lea otras preguntas en las etiquetas