Superposición de polígono espacial con cuadrícula y comprobación de qué coordenadas específicas del elemento de cuadrícula se encuentran utilizando R [cerrado]

¿Cómo se puede usar R para

1. dividir un archivo de forma en cuadrados de 200 metros / subpolígonos,
2. trace esta cuadrícula (incluidos los números de identificación para cada cuadrado) sobre el mapa original a continuación, y
3. evaluar en qué cuadrado se ubican las coordenadas geográficas específicas .

Soy un principiante en SIG y esta es quizás una pregunta básica, pero no he encontrado un tutorial sobre cómo hacer esto en R.

Lo que he hecho hasta ahora es cargar un archivo de forma de Nueva York y trazar algunas coordenadas geográficas ejemplares.

Estoy buscando un ejemplo (código R) de cómo hacerlo con los datos a continuación.

# 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")

Aquí hay un ejemplo usando un SpatialGridobjeto:

### 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
# ...

Ahora puede usar el overmétodo implementado para obtener las ID de celda:

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

Para trazar el archivo de forma y la cuadrícula con las ID de celda:

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

o sin color / clave de color:

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

El conjunto de datos de Nueva York proporcionado en la pregunta ya no está disponible para descargar. Utilizo el conjunto de datos nc del paquete sf para demostrar una solución usando el paquete sf:

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)

Si no ha mirado el paquete ráster R, tiene herramientas para convertir a / desde objetos SIG vectoriales, por lo que debería ser capaz de a) crear un ráster (cuadrícula) con celdas de 200x200m y b) convertirlo en un conjunto de polígonos con Una identificación lógica de algún tipo. A partir de ahí, miraría el paquete sp para ayudar con la intersección de los puntos y la cuadrícula de polígonos. Esta http://cran.r-project.org/web/packages/sp/vignettes/over.pdf podría ser un buen comienzo. Paseando por los documentos del paquete sp, puede comenzar con la clase SpatialGrid y simplemente omitir la parte ráster por completo.

El "universo SIG" es complejo y tiene muchos estándares que sus datos deben cumplir. Todas las "herramientas SIG" interoperan con los estándares SIG . Todos los "datos SIG serios" de hoy (2014) se almacenan en una base de datos .

La mejor manera de "usar R" en un contexto SIG, con otras herramientas FOSS , está integrada en SQL. Las mejores herramientas son PostgreSQL 9.X (ver PL / R ) y PostGIS .

Contestas:

• Para importar / exportar archivos de formas: use shp2pgsqlypgsql2shp .
• A "dividir un archivo de formas en 200 cuadrados metros / sub-polígonos": véase ST_SnapToGrid(), ST_AsRaster()etc. Necesitamos comprender mejor sus necesidades de expresar en una "receta".
• dices que necesitas "coordenadas geográficas ubicadas" .. quizás ST_Centroid()de los cuadrados (?) ... Puedes expresar "más matemáticamente" así que entiendo.

... Tal vez no necesite ninguna conversión de trama, solo una matriz de puntos regurlar muestreados.

Una forma primitiva es usar R sin PL / R , en su compilador externo habitual: solo convierta sus polígonos y exporte como forma o como WKT (ver ST_AsText), luego convierta los datos con awk u otro filtro al formato R.