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

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

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)

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

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