Instanziieren eines räumlichen Polygons ohne Verwendung eines Shapefiles in R

Die übliche Art, wie wir ein Shapefile in R lesen, ist über das maptools-Paket:

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

Ich habe jedoch einen Anwendungsfall, bei dem ich kein shapefile.shp habe, sondern stattdessen eine Reihe von Polygonkoordinaten

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

und seine entsprechende Projektion

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

Wie kann ich sfdata (die ein "Polygonobjekt" sein werden) direkt aus diesen Daten "instanziieren"? (ohne auf Umwegen ein Shapefile mit diesen Daten zu erstellen und dann aus dem neu erstellten Shapefile zu lesen)

Holen Sie sich zuerst die Koordinaten in eine 2-Spalten-Matrix:

> 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

Erstellen Sie dann ein Polygon, wickeln Sie es in ein Polygons-Objekt und anschließend in ein SpatialPolygons-Objekt um:

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

Der Grund für diese Komplexität ist, dass ein Polygon ein einfacher Ring ist, ein Polygon-Objekt mehrere Ringe mit einer ID (hier auf 1 gesetzt) ​​sein kann (wie ein einzelnes Feature in einem GIS) und ein SpatialPolygon ein CRS haben kann . Oh, ich sollte es wahrscheinlich einstellen:

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

Wenn Sie es in einen SpatialPolygonsDataFrame verwandeln möchten (was bei readShapeSpatial der Fall ist, wenn das Shapefile aus Polygonen besteht), gehen Sie wie folgt vor:

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

dies geben:

> 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 

Um die ausgezeichnete Antwort von Spacedman für den Fall zu vervollständigen, dass Ihre Daten mehrere Polygone enthalten würden, verwenden Sie folgenden Code 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)

Nur zum Spaß können Sie mit dem Diagramm vergleichen, das mit ggplot2dem anfänglichen Datenrahmen erhalten wurde:

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

Beachten Sie, dass der obige Code davon ausgeht, dass Sie nur ein Polyogn pro ID haben. Wenn einige IDs disjunkte Polygone hatten, sollte man dem Datensatz vermutlich eine weitere Spalte hinzufügen, zuerst group_bydie Unter-ID, dann group_by(upper-id)statt verwendenrowwise

Gleicher Code mit der purrr::mapFunktion:

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)}