Überlagern eines räumlichen Polygons mit einem Gitter und Prüfen, in welchen Gitterelementen sich spezifische Koordinaten befinden, mit R [closed]

32

Wie kann man R verwenden, um

  1. Teilen eines Shapefiles in 200-Meter-Quadrate / Unterpolygone,
  2. Zeichnen Sie dieses Raster (inkl. ID-Nummern für jedes Quadrat) über die Originalkarte unten und
  3. Auswerten, in welchem ​​Quadrat sich bestimmte geografische Koordinaten befinden .

Ich bin ein Anfänger in GIS und dies ist vielleicht eine grundlegende Frage, aber ich habe kein Tutorial dazu in R gefunden.

Bisher habe ich ein Shapefile von NYC geladen und einige beispielhafte geografische Koordinaten aufgezeichnet.

Ich suche ein Beispiel (R-Code), wie das mit den Daten unten geht.

# Load packages 
library(maptools)

# Download shapefile for NYC
# OLD URL (no longer working)
# shpurl <- "http://www.nyc.gov/html/dcp/download/bytes/nybb_13a.zip"
shpurl <- "https://www1.nyc.gov/assets/planning/download/zip/data-maps/open-data/nybb_13a.zip"

tmp    <- tempfile(fileext=".zip")
download.file(shpurl, destfile=tmp)
files <- unzip(tmp, exdir=getwd())

# Load & plot shapefile
shp <- readShapePoly(files[grep(".shp$", files)])
plot(shp)

# Define coordinates 
points_of_interest <- data.frame(y=c(919500, 959500, 1019500, 1049500, 1029500, 989500), 
                 x =c(130600, 150600, 180600, 198000, 248000, 218000),
                 id  =c("A"), stringsAsFactors=F)

# Plot coordinates
points(points_of_interest$y, points_of_interest$x, pch=19, col="red")

Bildbeschreibung hier eingeben

Majom
quelle
Siehe auch stackoverflow.com/q/17801398/287948
Peter Krauss

Antworten:

36

Hier ist ein Beispiel mit einem SpatialGridObjekt:

### read shapefile
library("rgdal")
shp <- readOGR("nybb_13a", "nybb")

proj4string(shp)  # units us-ft
# [1] "+proj=lcc +lat_1=40.66666666666666 +lat_2=41.03333333333333 
# +lat_0=40.16666666666666 +lon_0=-74 +x_0=300000 +y_0=0 +datum=NAD83
# +units=us-ft +no_defs +ellps=GRS80 +towgs84=0,0,0"

### define coordinates and convert to SpatialPointsDataFrame
poi <- data.frame(x=c(919500, 959500, 1019500, 1049500, 1029500, 989500),
                  y=c(130600, 150600, 180600, 198000, 248000, 218000),
                  id="A", stringsAsFactors=F)
coordinates(poi) <- ~ x + y
proj4string(poi) <- proj4string(shp)

### define SpatialGrid object
bb <- bbox(shp)
cs <- c(3.28084, 3.28084)*6000  # cell size 6km x 6km (for illustration)
                                # 1 ft = 3.28084 m
cc <- bb[, 1] + (cs/2)  # cell offset
cd <- ceiling(diff(t(bb))/cs)  # number of cells per direction
grd <- GridTopology(cellcentre.offset=cc, cellsize=cs, cells.dim=cd)
grd
# cellcentre.offset 923018 129964
# cellsize           19685  19685
# cells.dim              8      8

sp_grd <- SpatialGridDataFrame(grd,
                               data=data.frame(id=1:prod(cd)),
                               proj4string=CRS(proj4string(shp)))
summary(sp_grd)
# Object of class SpatialGridDataFrame
# Coordinates:
#      min     max
# x 913175 1070655
# y 120122  277602
# Is projected: TRUE
# ...

Jetzt können Sie die implementierte overMethode -method verwenden, um die Zellen-IDs abzurufen:

over(poi, sp_grd)
#   id
# 1 57
# 2 51
# 3 38
# 4 39
# 5 14
# 6 28

So zeichnen Sie das Shapefile und das Raster mit den Zellen-IDs:

library("lattice")
spplot(sp_grd, "id",
       panel = function(...) {
         panel.gridplot(..., border="black")
         sp.polygons(shp)
         sp.points(poi, cex=1.5)
         panel.text(...)
       })

spplot1

oder ohne Farbe / Farbschlüssel:

library("lattice")
spplot(sp_grd, "id", colorkey=FALSE,
       panel = function(...) {
         panel.gridplot(..., border="black", col.regions="white")
         sp.polygons(shp)
         sp.points(poi, cex=1.5)
         panel.text(..., col="red")
       })

spplot2

rcs
quelle
Das sieht nach einer Antwort für mich aus, aber für den Fall, dass Sie etwas anderes suchen. Probieren Sie das r-Tag in stackoverflow stackoverflow.com/search?q=R+tag
Brad Nesom
@rcs Dieser Code sieht genauso aus wie das, was ich versuche, aber mein Shapefile befindet sich in einer anderen Projektion: proj4string (DK_reg1) [1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" Hat jemand Vorschläge, wie man diese Shapefiles dieser Projektion in 1000 gleich große Gitterzellen aufteilt? und dann zufällig 100 davon auswählen und hervorheben?
I Del Toro
9

Der in der Frage angegebene New York-Datensatz steht nicht mehr zum Download zur Verfügung. Ich verwende den nc-Datensatz aus dem sf-Paket, um eine Lösung mit dem sf-Paket zu demonstrieren:

library(sf)
library(ggplot2)

# read nc polygon data and transform to UTM 
nc <- st_read(system.file('shape/nc.shp', package = 'sf')) %>%
  st_transform(32617)

# random sample of 5 points
pts <- st_sample(nc, size = 5) %>% st_sf

# create 50km grid - here you can substitute 200 for 50000
grid_50 <- st_make_grid(nc, cellsize = c(50000, 50000)) %>% 
  st_sf(grid_id = 1:length(.))

# create labels for each grid_id
grid_lab <- st_centroid(grid_50) %>% cbind(st_coordinates(.))

# view the sampled points, polygons and grid
ggplot() +
  geom_sf(data = nc, fill = 'white', lwd = 0.05) +
  geom_sf(data = pts, color = 'red', size = 1.7) + 
  geom_sf(data = grid_50, fill = 'transparent', lwd = 0.3) +
  geom_text(data = grid_lab, aes(x = X, y = Y, label = grid_id), size = 2) +
  coord_sf(datum = NA)  +
  labs(x = "") +
  labs(y = "")

# which grid square is each point in?
pts %>% st_join(grid_50, join = st_intersects) %>% as.data.frame

#>   grid_id                 geometry
#> 1      55 POINT (359040.7 3925435)
#> 2      96   POINT (717024 4007464)
#> 3      91 POINT (478906.6 4037308)
#> 4      40 POINT (449671.6 3901418)
#> 5      30 POINT (808971.4 3830231)

Bildbeschreibung hier eingeben

sebdalgarno
quelle
Vielen Dank. Ich habe den Link in meiner Frage aktualisiert, um die Änderungen auf der Webseite zu korrigieren. Jetzt sollte es wieder funktionieren.
Majom
Ich muss wirklich anfangen, das sfPaket zu benutzen . Das ist fantastisch!
philiporlando
Gibt es eine einfache Möglichkeit, nur die Gitterzellen zu zeichnen, die sich mit dem Zustandspolygon schneiden?
philiporlando
st_intersection (grid_50, nc) sollte es tun
sebdalgarno
Gibt es eine Möglichkeit, dasselbe zu replizieren, aber die Punkte in der Mitte jedes Gitters, so dass ein Gitter mit dem Lat / Long als Mittelpunkt des Gitters gezeichnet wird? @Sebdalgarno
Vijay Ramesh
2

Wenn Sie sich das R-Rasterpaket nicht angesehen haben, verfügt es über Tools zum Konvertieren in / aus Vektor-GIS-Objekten, sodass Sie a) ein Raster (Raster) mit 200 x 200 m Zellen erstellen und b) es in eine Gruppe von Polygonen mit konvertieren können eine logische ID irgendeiner Art. Von dort aus würde ich mir das sp-Paket ansehen, um beim Schneiden der Punkte und des Polygonrasters zu helfen. Diese http://cran.r-project.org/web/packages/sp/vignettes/over.pdf Seite könnte ein guter Anfang sein. Wenn Sie durch die sp-Paketdokumente wandern, können Sie möglicherweise mit der SpatialGrid-Klasse beginnen und den Raster-Teil einfach vollständig überspringen.

cokrzys
quelle
-1

Das "GIS-Universum" ist komplex und hat viele Standards, denen Ihre Daten entsprechen müssen. Alle "GIS-Tools" arbeiten nach GIS-Standards zusammen . Alle "seriösen GIS-Daten" von heute (2014) werden in einer Datenbank gespeichert .

Die beste Möglichkeit, R in einem GIS-Kontext mit anderen FOSS- Tools zu verwenden, ist die Einbettung in SQL. Die besten Tools sind PostgreSQL 9.X (siehe PL / R ) und PostGIS .


Du antwortest:

  • So importieren / exportieren Sie Formdateien: Verwenden Sie shp2pgsqlundpgsql2shp .
  • Um „aufzuteilen eine Formdatei in 200 Metern Quadrate / Subpolygone“: siehe ST_SnapToGrid(), ST_AsRaster()etc. Wir müssen besser verstehen , Ihre Bedürfnisse in ein „Rezept“ zum Ausdruck bringen.
  • Sie sagen, dass "geografische Koordinaten gefunden werden müssen" .. vielleicht ST_Centroid()von den Quadraten (?) ... Sie können "mehr mathematisch" ausdrücken, damit ich verstehe.

... Möglicherweise benötigen Sie keine Rasterkonvertierung, sondern nur eine Matrix von Punkten mit Regurlar-Abtastung.


Ein primitiver Weg ist die Verwendung von R ohne PL / R in einem üblichen externen Compiler: Konvertieren Sie nur Ihre Polygone und exportieren Sie sie als Form oder als WKT (siehe ST_AsText). Konvertieren Sie dann Daten mit awk oder einem anderen Filter in das R-Format.

Peter Krauss
quelle
1
Danke für Ihre Hilfe. Ich würde jedoch eine Lösung vorziehen, die sich vollständig auf R und vorhandene Pakete stützt. Wenn ich die Formdatei in 200 mx 200 m Teilpolygone aufteilen kann, kann ich überprüfen, mit point.in.polygonwelchen Koordinaten sich welche Polygone befinden. Mein Problem ist, das ursprüngliche Shapefile in diese Unterpolygone aufzuteilen.
9.