Recently researchers have been creating grids for analyses of various shapes. We already explored how to create a hexagonal grid but now we will learn how to create a square grid within the extent of a pre-defined study area. This method requires a few more steps but square polygon grids and the resulting adjacency matrix are common in disease epidemiology and will be used in future exercises.
- Exercise 1.9 - Download and extract zip folder into your preferred location
- Set working directory to the extracted folder in R under File - Change dir...
- First we need to load the packages needed for the exercise
library(sp)
library(rgdal)
library(raster)
library(adehabitatMA) - Now open the script "GridSystem2Script.R" and run code directly from the script
- We need to have all layers in same projection so import, create, and remove outliers for mule deer locations then project all to the Albers projection as we did previously.
muleys <-read.csv("muleysexample.csv", header=T)
#Make a spatial data frame of locations after removing outliers
summary(muleys$id)
str(muleys)
#Remove outlier locations
newmuleys <-subset(muleys, muleys$Long > -110.50 & muleys$Lat > 37.8
& muleys$Long < -107)
muleys <- newmuleys
coords<-data.frame(x = muleys$Long, y = muleys$Lat)
crs<-"+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
head(coords)
deer.spdf <- SpatialPointsDataFrame(coords= coords, data = muleys, proj4string =
CRS(crs))
head(deer.spdf)
proj4string(deer.spdf)
#Project the deer.spdf to Albers as in previous exercise
Albers.crs <-CRS("+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-
96+x_0=0+y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs")
deer.albers <-spTransform(deer.spdf, CRS=Albers.crs)
proj4string(deer.albers)
bbox(deer.albers)
# min max
#x -1145027 -1106865
#y 1695607 1729463 - Create points for x and y from the bounding box of all mule deer locations with 1500 m spacing between each point.
plot(deer.albers)
## create vectors of the x and y points
x <- seq(from = -1145027, to = -1106865, by = 1500)
y <- seq(from = 1695607, to = 1729463, by = 1500) - Create a grid of all pairs of coordinates (as a data.frame) using the "expand grid" function and then make it a gridded object.
xy <- expand.grid(x = x, y = y)
class(xy)
str(xy)
#Identifiy projection before creating Spatial Points Data Frame
crs2<-"+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0
+y_0=0+ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
grid.pts<-SpatialPointsDataFrame(coords= xy, data=xy, proj4string = CRS(crs2))
plot(grid.pts)
gridded(grid.pts)
class(grid.pts)
#Make
points a gridded object (i.e., TRUE or FALSE)
gridded(grid.pts) <- TRUE
gridded(grid.pts)
str(grid.pts)
plot(grid.pts) - Make the grid of points into a Spatial Polygon then convert the spatial polygons to a SpatialPolygonsDataFrame.
grid <- as(grid.pts, "SpatialPolygons")
plot(grid)
str(grid)
class(grid)
summary(grid)
gridspdf <- SpatialPolygonsDataFrame(grid, data=data.frame(id=row.names(grid),
row.names=row.names(grid)))
names.grd<-sapply(gridspdf@polygons, function(x) slot(x,"ID"))
text(coordinates(gridspdf), labels=sapply(slot(gridspdf, "polygons"), function(i) slot(i,
"ID")), cex=0.3)
points(deer.albers, col="red")
str(gridspdf@polygons) - Similar to the hexagonal grid, identify the cell ID that contains each mule deer location.
o = over(deer.albers,gridspdf)
head(o)
new = cbind(deer.albers@data, o)
head(new) - We get some NA errors because our grid does not encompass all mule deer locations so expand the grid then re-run the code over from xy through new2 again.
x <- seq(from = -1155027, to = -1106865, by = 1500)
y <- seq(from = 1695607, to = 1739463, by = 1500)
##BE SURE TO RUN CODE FROM XY CREATION THROUGH NEW2 AGAIN THEN BEFORE CODE BELOW!!
o2 = over(deer.albers,gridspdf)
head(o2)
new2 = cbind(deer.albers@data, o2)#No more NAs causing errors!
new2[1:15,] - Now we can load a vegetation raster layer textfile clipped in ArcMap to summarize vegetation categories within each polygon grid cell.
veg <-raster("ExtentNLCD2.txt")
plot(veg)
class(veg) - Clip the raster within the extent of the newly created grid
bbclip <- crop(veg, gridspdf)
plot(bbclip)
points(deer.albers, col="red")
plot(gridspdf, add=T)
#Cell size of raster layer
xres(bbclip)#shows raster cell size
#Create histogram of vegetation categories in bbclip
hist(bbclip)
#Calculate the size of each cell in your square polygon grid
ii <- calcperimeter(gridspdf)#requires adehabitatMA package
as.data.frame(ii[1:5,])#Identifies size of only the first 5 grid cells - We can extract the vegetation characteristics within each polygon of the grid.
table = extract(bbclip,gridspdf) - We can then tabulate area of each vegetation category within each polygon by extracting vegetation within each polygon by ID then appending the results back to the extracted table by running it twice but with different names. Summarizing the vegetation characteristics in each cell will be used in future resource selection analysis or disease epidemiology.
table = extract(bbclip,gridspdf)
str(table)
area = extract(bbclip,gridspdf)
combine=lapply(area,table)
combine
combine[[1]]#Shows vegetation categories and numbers of cells in grid #1
combine[[27]]21 22 31 42 52 82 90 38 7 23 392 1883 11 146
21 42 52 81 82 101 69 279 5 2046
