This tutorial shows how to divide a spatial polygon into a grid. There are various ways to achieve this, and we will introduce several of them.

We use the London Borough of Camden as an example and the map projection is the British National Grid. Our method can also be applied to the unprojected spatial objects, for example, under the coordination system of WGS84.

First, the shapefile is read into R:

require(rgdal)
require(raster)
camden <- readOGR("camdenBoundary_nationalGrid.shp", stringsAsFactors = F)
## OGR data source with driver: ESRI Shapefile 
## Source: "camdenBoundary_nationalGrid.shp", layer: "camdenBoundary_nationalGrid"
## with 1 features
## It has 9 fields
plot(camden)

Dividing the grid

We would like the cells that have a length of 500 m on east-west direction and a length of 200 m on north-west direction. The codes are:

# note the order of resolution
grid <- raster(extent(camden), resolution = c(500,200), crs = proj4string(camden))

# class of grid
class(grid)
## [1] "RasterLayer"
## attr(,"package")
## [1] "raster"
# convert to SpatialPolygonsDataFrame
gridPolygon <- rasterToPolygons(grid)
class(gridPolygon)
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
plot(camden)
plot(gridPolygon, add = T)

The result is not perfect. The grid fails to cover some parts of the region, such as tips on the north and east. An improvement is to extend the grid on each direction by one cell:

grid <- raster(extent(camden), resolution = c(500,200), crs = proj4string(camden))
grid <- raster::extend(grid, c(1,1))

# convert to SpatialPolygonsDataFrame
gridPolygon <- rasterToPolygons(grid)
class(gridPolygon)
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
plot(camden)
plot(gridPolygon, add = T)

It is nearly done! Now you can save gridPolygon as a shapefile, via the writeOGR function.

Coordinates of the cell centers

If you want to get the coordinates of the cell centers and then create a Spatial object, it can be done by the rasterToPoints function.

# a SpatialPoints object. The same projection as grid is used.
rasterCent <- rasterToPoints(grid, spatial = T)
plot(gridPolygon)
plot(rasterCent, add = T, col = "red", pch = 18)

# a matrix
matRasterCent <- rasterToPoints(grid, spatial = F)
head(matRasterCent)
##             x        y
## [1,] 523701.7 187703.7
## [2,] 524201.7 187703.7
## [3,] 524701.7 187703.7
## [4,] 525201.7 187703.7
## [5,] 525701.7 187703.7
## [6,] 526201.7 187703.7

Clipping the grid to fit the polygon boundary

The grid has been generated. If you want to keep only the cells that only intersect the polygon, or to clip the cells according to the polygon boundary, it is possible.

Clipping the intersected cells can be done as follows:

intersectGridClipped <- raster::intersect(gridPolygon, camden)
## Loading required namespace: rgeos
plot(intersectGridClipped)

Keeping the interested cells without clipping them can be done, via adding the id attribute to the cells:

gridPolygon$id <- 1:nrow(gridPolygon)
intersectGridClipped <- raster::intersect(gridPolygon, camden)
intersectGrid <- gridPolygon[gridPolygon$id %in% intersectGridClipped$id, ]
plot(intersectGrid)