Vjylku

How can one use R to

1. 
   split a shapefile in 200 meter squares/sub-polygons,


2. 
   plot this grid (incl. ID numbers for each square) over the original map below, and


3. evaluate in which square specific geographic coordinates are located.

I am a beginner in GIS and this is perhaps a basic question, but I haven't found a tutorial on how to do this in R.

What I have done so far is loading a shapefile of NYC and plotting some exemplary geographic coordinates.

I appreciate any general advice, but I am looking for an example (R code) how to this with the data below.

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

Here is an example using a SpatialGrid object:

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

# ...

Now you can use the implemented over-method to obtain the cell IDs:

over(poi, sp_grd)

#   id

# 1 57

# 2 51

# 3 38

# 4 39

# 5 14

# 6 28

To plot the shapefile and the grid with the cell IDs:

library("lattice")

spplot(sp_grd, "id",

       panel = function(...) {

         panel.gridplot(..., border="black")

         sp.polygons(shp)

         sp.points(poi, cex=1.5)

         panel.text(...)

       })

or without colour/colour key:

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

       })

The New York dataset provided in the question is no longer available for download. I use the nc dataset from sf package to demonstrate a solution using sf package:

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)

Have you looked at the R raster package? It has tools to convert to/from vector GIS objects so you should be able to a) create a raster (grid) with 200x200m cells and b) convert it to a set of polygons with a logical id of some kind. From there I would look at the sp package to help with intersecting the points and the polygon grid. This http://cran.r-project.org/web/packages/sp/vignettes/over.pdf page might be a good start. Wandering through the sp package docs you might be able to start with the SpatialGrid-class and just skip the raster part entirely.

The "GIS universe" is complex and have many standards that your data must be compliant. All "GIS tools" interoperates by GIS-standards. All "serious GIS data" today (2014) are stored in a database.

The best way to "use R" in a GIS context, with other FOSS tools, is embedded into SQL. The best tools are PostgreSQL 9.X (see PL/R) and PostGIS.

You answer:

• To import/export shape files: use shp2pgsql and pgsql2shp.


• To "split a shape file in 200 meter squares/sub-polygons": see ST_SnapToGrid(), ST_AsRaster(), etc. We need understand better your needs to express into a "recipe".


• you say that need "geographic coordinates are located" .. perhaps ST_Centroid() of the squares (?)... You can express "more mathematically" so I understand.

... Perhaps you not need any raster convertion, only a matrix of regurlar-sampled points.

A primitive way is use R without PL/R, in a your usual external compiler: only convert your polygons and export as shape or as WKT (see ST_AsText), then convert data with awk or another filter to the R format.

Is there way to dropping empty grids? or just griding the boundary of interest?