Superponga un polígono espacial con una cuadrícula y verifique en qué elemento de la cuadrícula se ubican las coordenadas específicas

30

¿Cómo puede uno usar R para

  1. dividir un shapefile en cuadrados / sub-polígonos de 200 metros,
  2. traza esta cuadrícula (incluidos los números de identificación para cada cuadrado) sobre el mapa original a continuación, y
  3. evalúa en qué cuadrado se ubican las coordenadas geográ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 shapefile de NYC y trazar algunas coordenadas geográficas ejemplares.

Agradezco cualquier consejo general, pero Estoy buscando un ejemplo (código R) sobre 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")

    
pregunta majom 06.03.2014 - 23:02

4 respuestas

34

Aquí hay un ejemplo que usa un objeto SpatialGrid :

### 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 método implementado over 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 shapefile y la cuadrícula con los identificadores 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(...)
       })

osincolor/tecladecolor:

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

    
respondido por el rcs 10.03.2014 - 22:53
5

El conjunto de datos de Nueva York que se proporciona en la pregunta ya no está disponible para descargar. 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)

    
respondido por el sebdalgarno 05.03.2018 - 02:47
2

¿Has mirado el paquete R raster? Tiene herramientas para convertir a / desde objetos vectoriales GIS, por lo que debería a) crear un ráster (cuadrícula) con 200x200m de celdas yb) convertirlo en un conjunto de polígonos con una identificación lógica de algún tipo. Desde allí miraría el paquete sp para ayudar a intersectar los puntos y la cuadrícula de polígonos. Esta página enlace puede ser un buen comienzo. Deambulando por los documentos del paquete sp puede comenzar con la clase SpatialGrid y simplemente omitir la parte ráster por completo.

    
respondido por el cokrzys 10.03.2014 - 05:25
0

El "universo GIS" es complejo y tiene muchos estándares que sus datos deben cumplir. Todas las "herramientas de GIS" interactúan mediante normas de GIS . Todos los "datos GIS serios" de hoy (2014) están almacenados en una base de datos .

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

Usted responde:

  • Para importar / exportar archivos de formas: use shp2pgsql y pgsql2shp .
  • Para "dividir un archivo de forma en cuadrados / sub-polígonos de 200 metros": vea ST_SnapToGrid() , ST_AsRaster() , etc. Necesitamos entender mejor sus necesidades para expresar en una "receta".
  • dices que necesitas "las coordenadas geográficas están ubicadas" ... quizás ST_Centroid() de los cuadrados (?) ... Puedes expresar "más matemáticamente" para que yo entienda.

... Tal vez no necesite ninguna conversión ráster, solo una matriz de puntos de muestra regurlar.

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

    
respondido por el Peter Krauss 09.03.2014 - 22:20

Lea otras preguntas en las etiquetas