Share

1.9 Creating a Square Polygon Grid Over a Study Area

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.

  1. Exercise 1.9 - 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(rgdal)
    library(raster)
    library(adehabitatMA)
  4. Now open the script "GridSystem2Script.R" and run code directly from the script
  5. 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)
    summary(muleys$id)
    str(muleys)
    #Remove outlier locations
    newmuleys <-subset(muleys, muleys$Long > -110.50 & muleys$Lat > 37.8
    & muleys$Long < -107)
    muleys <- newmuleys

    #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"
    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
  6. 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)
  7. 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)
  8. 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)
  9. 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)
  10. 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,]
  11. 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)
  12. 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
  13. We can extract the vegetation characteristics within each polygon of the grid.

    table = extract(bbclip,gridspdf)
  14. 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
    21 22 31 42 52 82 90
    38 7 23 392 1883 11 146
    combine[[27]]
    21 42 52 81 82
    101 69 279 5 2046