1.8 Creating a Hexagonal Polygon Grid Over a Study Area

Numerous research objectives require the need for creating a grid system of equal size over a study site such as studies on resource selection and disease epidemiology. Grid systems overlayed on a study site typically are shapefiles that can either be created and imported from GIS software or created in R. Considering we have already learned how to import shapefiles,
we will explore how to create grids in R for this section. Grids can be of any size and shape but should be based on something biologically meaningful to the animal or system you are studying. For example, disease epidemiology studies often base the size of the grid cell on the dailiy movement distance or home range of the study animal if that data is known (Farnsworth et al. 2006, Rees et al. 2011).

1. Exercise 1.8 - Download and extract zip folder into your preferred location
2. Set working directory to the extracted folder in R under File - Change dir...
3. First we need to load the packages needed for the exercise
   
   library(sp)
   library(lattice)
   library(rgdal)
   library(rgeos)
   library(raster)
4. Now open the script "GridScripts.R" and run code directly from the script
5. Also need to import several shapefiles for mule deer from Section 1.3
   
   study.counties<-readOGR(dsn=".",layer="MDcounties")
   str(study.counties) #Identifies 5 slots for the shapefile (data, polygons, order, bbox,
   and proj4string)
   class(study.counties)#Shows class and package used
   proj4string(study.counties) #Shows projection information
   plot(study.counties)#plots study sites on map
   study.counties@data$StateCO #Displays labels for counties in plot
   #Labels each county with @plotOrder of each polygon (i.e., county)
   text(coordinates(study.counties), labels=sapply(slot(study.counties, "polygons"),
   function(i) slot(i, "ID")), cex=0.8)
   
   #NOTE: This can be any column or label within your shapefile
   
   muleys <-read.csv("muleysexample.csv", header=T)
   str(muleys)
   
   #Remove outlier locations
   

   newmuleys <-subset(muleys, muleys$Long > -110.50 & muleys$Lat > 37.8
   & muleys$Long < -107)
   muleys <- newmuleys

6. Identify the columns with coordinates then make a spatial data frame of locations after removing outliers
   
   coords<-data.frame(x = muleys$Long, y = muleys$Lat)
   crs<-"+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
   coords
   deer.spdf <- SpatialPointsDataFrame(coords= coords, data = muleys, proj4string =
   CRS(crs))
   deer.spdf[1:5,]
   class(deer.spdf)
   proj4string(deer.spdf)
   points(deer.spdf,col="red")
7. Rename labels by county name otherwise plot order would be used because duplicate counties within each state (i.e., CO, UT) occured in original shapefile from ArcMap
   
   row.names(study.counties)<-as.character(study.counties$StateCO)
   str(study.counties@polygons[3], max.level=3)
   
   #Now add labels of State and County to Map
   text(coordinates(study.counties), labels=sapply(slot(study.counties, "polygons"),
   function(i) slot(i, "ID")), cex=0.3)
8. Now lets extract counties within the extent of the mule deer locations
   
   int <- gIntersection(study.counties,deer.spdf)#requires rgeos library
   clipped <- study.counties[int,]
   MDclip <- as(clipped, "SpatialPolygons")
   plot(MDclip,pch=16)
   #Now add labels of State and County to Map
   text(coordinates(MDclip), labels=sapply(slot(MDclip, "polygons"), function(i) slot(i, "ID")),
   cex=0.8)
9. We also can create a hexagonal grid across the study site
   
   HexPts <-spsample(MDclip, , n=1000, offset=c(0,0))
   HexPols <- HexPoints2SpatialPolygons(HexPts)
   proj4string(HexPols) <- CRS(crs)
   plot(HexPols, add=T)
10. Let's create this hexagonal grid across our study site by zooming into deer locations from Section 1.3.
    
    #Import the study site zoomed in shapefile
    study.zoom<-readOGR(dsn=".",layer="MDzoom")
    plot(study.zoom,pch=16)
    points(deer.spdf,col="red")
    
    #Create new hexagonal grid
    HexPts2 <-spsample(study.zoom, , n=500, offset=c(0,0))
    HexPols2 <- HexPoints2SpatialPolygons(HexPts2)
    proj4string(HexPols2) <- CRS(crs)
    plot(HexPols2, add=T)
    #Now add labels to each hexagon for unique ID
    text(coordinates(HexPols2), labels=sapply(slot(HexPols2, "polygons"), function(i) slot(i,
    "ID")), cex=0.3)
11. We can intersect the mule deer locations with the polygon shapefile (i.e., county) they occured in if needed
    
    o = over(deer.spdf,study.counties) #By county locations occurs in
    new = cbind(deer.spdf@data, o)
    head(o)
    head(deer.spdf)
    head(new)
    
    #Used to rename labels by hexagonal grid ID only for visualization only!
    row.names(HexPols2)<-as.character(HexPols2@plotOrder)
    
    
12. As an aside, let's explore how to assign the area a location occurs in by intersecting points within the polygon shapefile.
    
    o2 = over(deer.spdf,HexPols2)
    o2
    new2 = cbind(deer.spdf@data,o2)
    head(new2)
    new2
    deer.spdf@data[1:10,]
    HexPols2
    
    #Now plot with new grid IDs
    plot(study.zoom,pch=16)
    points(deer.spdf,col="red")
    plot(HexPols2, add=T)
    
    #Now add labels of State and County to Map
    text(coordinates(HexPols2), labels=sapply(slot(HexPols2, "polygons"), function(i) slot(i,
    "ID")), cex=0.3)
13. As an alternative to importing a polygon that we created in ArcMap, we can create a polygon in R using the coordinates of the boundary box of the area of interest. In our case here, the bounding box will be the mule deer locations.
    
    #First we need to create the polygon within the extent of our mule deer locations
    proj4string(deer.spdf)
    bbox(deer.spdf@coords)
    bb <- cbind(x=c(-108.69984,-108.69984,-109.14286,-109.14286, -108.69984),
    y=c(37.63163, 37.92621,37.92621,37.63163,37.63163))
    SP <- SpatialPolygons(list(Polygons(list(Polygon(bb)),"1")),
    proj4string=CRS(proj4string(MDclip)))
    plot(SP)
    proj4string(SP)
    points(deer.spdf,col="red")
    
    
14. Now let's make practical use of the new bounding box we created by clipping a larger raster dataset. A smaller raster dataset runs analyses faster, provides a zoomed in view of mule deer locations and vegetation, and is just easier to work with.
    
    #Load vegetation raster layer textfile from ArcMap
    veg <-raster("extentnlcd2.txt")
    plot(veg)
    class(veg)
    
    #Clip using the raster imported with "raster" package
    bbclip <- crop(veg, SP)
    veg
    #WON'T WORK because projections are not the same, WHY?
    
    #Let's check projections of layers we are working with now.
    proj4string(MDclip)
    proj4string(deer.spdf)
    proj4string(SP)
    proj4string(veg)
15. We need to have all layers in same projection so let's project the deer.spdf to Albers and then clip vegetation layer with new polygon we created in the Albers projection.
    
    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)
    points(deer.albers, col="red")
    class(deer.albers)
    proj4string(deer.albers)
    head(deer.spdf)
    head(deer.albers)
    
    #Now determine the new coordinates and create a new polygon to clip the raster.
    bbox(deer.albers)
    bb1 <- cbind(x=c(-1106865,-1106865,-1145027,-1145027, -1106865),y=c(1695607,
    1729463,1729463,1695607,1695607))
    AlbersSP <- SpatialPolygons(list(Polygons(list(Polygon(bb1)),"1")),
    proj4string=CRS(proj4string(deer.albers)))
    
    #Check to see all our layers are now in Albers projection
    plot(AlbersSP)
    proj4string(veg)
    proj4string(deer.albers)
    proj4string(AlbersSP)
    
    plot(veg)
    points(deer.albers, col="red")
    
    #Clip using the raster imported with "raster" package
    bbclip <- crop(veg, AlbersSP)
    plot(bbclip)
    points(deer.albers, col="red")
    plot(AlbersSP, lwd=5, add=T)
    #text(coordinates(AlbersSP), labels="Colorado Mule Deer")