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

Question

¿Cómo se puede utilizar R a

1. dividir un archivo shapefile en cuadrados/subpolígonos de 200 metros,
2. trazar esta cuadrícula (incluidos los números de identificación de cada casilla) sobre el mapa original que figura a continuación, y
3. evaluar en qué plaza específica las coordenadas geográficas se encuentran .

Soy principiante en SIG y quizás sea 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 shapefile de NYC y trazar algunas coordenadas geográficas ejemplares.

Estoy buscando un ejemplo (código R) cómo esto 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")

Answer 1

He aquí un ejemplo en el que se utiliza un SpatialGrid objeto:

### 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 utilizar el over -para obtener los ID de las celdas:

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

Para trazar el shapefile y la cuadrícula con los 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")
       })

Answer 2

El conjunto de datos de Nueva York proporcionado en la pregunta ya no está disponible para su descarga. Utilizo el conjunto de datos nc del paquete sf para demostrar una solución utilizando 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)

Answer 3

Si no has mirado el paquete R raster, tiene herramientas para convertir a/desde objetos SIG vectoriales, así que deberías ser capaz de a) crear un raster (cuadrícula) con celdas de 200x200m y b) convertirlo en un conjunto de polígonos con un id lógico de algún tipo. A partir de ahí, yo miraría el paquete sp para ayudar con la intersección de los puntos y la cuadrícula de polígonos. Este http://cran.r-project.org/web/packages/sp/vignettes/over.pdf puede ser un buen comienzo. Deambulando por la documentación del paquete sp podrías empezar con la clase SpatialGrid y saltarte la parte raster por completo.

Answer 4

El "universo SIG" es complejo y tiene muchas normas que sus datos deben cumplir. Todas las "herramientas SIG" interoperan mediante Normas SIG . Todos los "datos SIG serios" actuales (2014) son almacenados en una base de datos .

La mejor manera de "utilizar R" en un contexto de SIG, con otros FOSS herramientas, está integrado en SQL. Las mejores herramientas son PostgreSQL 9.X (véase PL/R ) y PostGIS .

---

Contesta tú:

• Para importar/exportar archivos shape: utilice shp2pgsql y pgsql2shp .
• Para "dividir un archivo shape en cuadrados/subpolígonos de 200 metros": véase ST_SnapToGrid() , ST_AsRaster() etc. Necesitamos entender mejor sus necesidades para expresarlas en una "receta".
• usted dice que necesita "coordenadas geográficas se encuentran" tal vez ST_Centroid() de los cuadrados (?)... Usted puede expresar "más matemáticamente" para que yo entienda.

... Tal vez usted no necesita ninguna conversión raster, sólo una matriz de puntos regurlar-muestreados.

---

Una forma primitiva es utilizar R sin PL/R en su compilador externo habitual: sólo convierta sus polígonos y expórtelos como shape o como WKT (véase ST_AsText ) y, a continuación, convertir los datos con awk u otro filtro al formato R.