# Overlay a spatial polygon with a grid and check in which grid element specific coordinates are located

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)

# OLD URL (no longer working)

tmp    <- tempfile(fileext=".zip")
files <- unzip(tmp, exdir=getwd())

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

• – Peter Krauss Mar 22 '14 at 7:51
• Is there way to dropping empty grids? or just griding the boundary of interest? – iHermes Apr 6 at 22:53
• If you have a new question, please ask it by clicking the Ask Question button. Include a link to this question if it helps provide context. - From Review – BBG_GIS Apr 6 at 23:16

Here is an example using a `SpatialGrid` object:

``````### read shapefile
library("rgdal")

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

• This looks like an answer to me, but in case you are looking for something different. Try the r tag in stackoverflow stackoverflow.com/search?q=R+tag – Brad Nesom Mar 14 '14 at 2:26
• @rcs this code looks just like what I am trying to do but my shapefile is in a different projection: `proj4string (DK_reg1) [1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"` does anyone have any suggestions on how to break this shapefiles of this projection to 1000 equal sized grid cells? and then randomly select 100 of them and highlight them? – I Del Toro Sep 13 '14 at 11:34

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

• Thanks. I updated the link in my question to relfect the changes on their webpage. Now it should work again. – majom Mar 5 '18 at 15:18
• I really need to start using the `sf` package. This is awesome! – philiporlando Mar 23 '18 at 20:41
• Is there an easy way to only plot the grid cells that intersect with the state polygon? – philiporlando Mar 23 '18 at 22:09
• st_intersection(grid_50, nc) should do it – sebdalgarno Mar 24 '18 at 23:18
• Is there a way to replicate the same, but the points in the centre of each grid, so a grid is being drawn with the lat/long as the centre of the grid @sebdalgarno – Vijay Ramesh Aug 15 at 17:58

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.

• Thanks for your comment. I already looked at those packages, but failed to implement this. This was why I was asking for an example. – majom Mar 10 '14 at 9:05

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.

• 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.
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.
• Thanks for your help. However, I would strongly prefer a solution which relies completely on R and existing packages. When I am able to split the shape file in 200m*200m subpolygons I can check with `point.in.polygon` which coordinates are in which polygons. My problem is to split the original shapefile in those sub-polygons. – majom Mar 9 '14 at 21:38