R – Generating Grid Shapefile Using Overlay and Vector Grid

overlayrvector-grid

I have the following shapefile of a 10×10 degree latitude/longitude grid that I created in QGIS. I can read it into R using the rgdal package.

Grid<-readOGR(".","GridShapeFile")

It has the following attributes and structure. 

summary(Grid)

Object of class SpatialPolygonsDataFrame
Coordinates:
   min max
x -180 190
y -100  90
Is projected: FALSE 
proj4string :
[+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0]
Data attributes:
       ID             XMIN           XMAX           YMIN           YMAX    
 Min.   :  0.0   Min.   :-180   Min.   :-170   Min.   :-100   Min.   :-90  
 1st Qu.:175.5   1st Qu.: -90   1st Qu.: -80   1st Qu.: -60   1st Qu.:-50  
 Median :351.0   Median :   0   Median :  10   Median : -10   Median :  0  
 Mean   :351.0   Mean   :   0   Mean   :  10   Mean   : -10   Mean   :  0  
 3rd Qu.:526.5   3rd Qu.:  90   3rd Qu.: 100   3rd Qu.:  40   3rd Qu.: 50  
 Max.   :702.0   Max.   : 180   Max.   : 190   Max.   :  80   Max.   : 90 

# An example row of the data
Grid[50,]
class       : SpatialPolygonsDataFrame 
features    : 1 
extent      : -60, -50, 70, 80  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84+towgs84=0,0,0 
variables   : 5
names       : ID, XMIN, XMAX, YMIN, YMAX 
min values  : 49,  -60,  -50,   70,   80 
max values  : 49,  -60,  -50,   70,   80

I would like to learn how to create this same grid/shapefile entirely within R, rather than in QGIS, as well as grids of other sizes (1×1,5×5,etc.). The ultimate goal is to use over() in the sp package to overlay this grid onto another file of polygons and count the number of grids intersected by each polygon. I already know how to do this next step (I think), I just want to know how to generate the grid layer in R.

Best Answer

Take a look at the raster function in the raster package. It will let you create a raster with a specified extent, number of rows/columns and resolution.

Here I will use characteristics of your data summary to create a 100x100 raster within the specified extent. I am passing an extent object to define the x and y limits. You can also use the specific arguments (xmn, xmx, ymn, ymx) within the raster function. 

library(raster) 
library(sp) 

r <- raster(extent(matrix( c(-180, -100, 190,  90), nrow=2)), nrow=100, ncol=100, 
            crs = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")            
  r[] <- 1:ncell(r)
  summary(r)
  print(r)
  plot(r)

It is simple to coerce raster objects to a gridded sp object using; 

sp.r <- as(r, "SpatialPixelsDataFrame")
  class(sp.r)   
  spplot(sp.r, "layer")
Related Question