Making a ‘heatmap’ of volume (or another inventory variable) is a straightforward process in R. Here’s how to go about turning your plot data into a mapped visualization of stand stocking.
Begin by making sure all of the packages you might need are installed and loaded into your R environment.
library(dplyr) library(rgdal) library(raster) library(gstat) library(ggplot2) library(ggmap)
You’ll want to make sure you have a file with your plot data in it- the file you start with should look something like this:
## PlotNumber sawVolume__doyle ## 1 1 1400.238 ## 2 2 5892.328 ## 3 3 9778.634 ## 4 4 8308.458 ## 5 5 5534.393 ## 6 6 5786.646
Next, read in a shapefile or kml file of your plot data. We’ll do that using the
rgdal package. There are numerous ways of working with GIS data in R- this is just one approach.
plotLocs <- readOGR(dsn = '/tmp/plot_locations.shp', layer = 'plot_locations') plotLocs@data <- dplyr::left_join(plotLocs@data, sawByPlot)
This results in a SpatialPointsDataFrame object, that includes coordinates for each plot as well as the saw volume for each plot from our column
sawVolume__doyle. The object also has a specific projection with it, which we can see by calling
## CRS arguments: ## +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
We’ll also make an
dataframe object with the plot locations and saw volume, for making a map later.
plotCoordDF <- data.frame(plotLocs@coords) %>% rename(x = coords.x1, y = coords.x2) plotCoordDF$PlotNumber <- plotLocs@data$PlotNumber plotCoordDF$sawVolume__doyle <- plotLocs@data$sawVolume__doyle head(plotCoordDF)
## x y PlotNumber sawVolume__doyle ## 1 -83.08880 37.47427 54 4420.695 ## 2 -83.09897 37.47427 42 5148.990 ## 3 -83.09388 37.47427 48 4431.422 ## 4 -83.09812 37.47427 43 1055.523 ## 5 -83.09304 37.47427 49 3354.248 ## 6 -83.09049 37.47427 52 3981.321
We also need to read in a shapefile or other GIS file with the boundary of the stand in which these plots are located. The projection of the boundary file must match the projection of the plots.
standBounds <- readOGR(dsn = '/tmp/stand_bounds.shp', layer = 'stand_bounds')
We will then ‘rasterize’ the stand bounds, to create the framework in which we’ll build the heatmap. You can vary the
ncol depending on the size of your stand and plots.
First, we create a raster that circumscribes the stand, with every pixel carrying a value of 1, and then define the raster to cover only the stand bounds by setting the extent. We also create a
SpatialPixelsDataFrame object from the raster to make predictions easier, and a
data.frame object so we can build a map inside R.
grid <- raster(ncol = 100, nrow = 100) extent(grid) <- extent(standBounds) rp <- rasterize(standBounds, grid, field = 1) rastSPDF <- as(rp, "SpatialPixelsDataFrame") rastDF <- as.data.frame(rastSPDF)
Now we’re ready to use the inverse distance weighting function to create a wall-to-wall surface of volume from the plot data. In the
idw function, the
nmax value can be changed to adjust the number of points used to calculate the inverse distance weighted mean- the smaller the
nmax, the more coarse the final map appears.
idwModel <- idw(sawVolume__doyle ~ 1, plotLocs, newdata = rastSPDF, nmax = 10)
## [inverse distance weighted interpolation]
rastSPDF@data$layer <- idwModel$var1.pred rastDF$totSawIDW <- rastSPDF$layer
We can use plotting functions to see how the idw predictions look. Here’s an example using the
maxSaw <- max(plotCoordDF$sawVolume__doyle) ggplot(data = rastDF, aes(x = x, y = y)) + geom_tile(aes(fill = totSawIDW)) + coord_equal() + geom_path(data = standBounds, aes(x = long, y = lat)) + geom_point(data = plotCoordDF, aes(color = sawVolume__doyle)) + scale_color_gradient(limits = c(0, maxSaw), "Saw Volume (doyle bd ft) per acre") + scale_fill_gradient(limits = c(0, maxSaw), "Saw Volume (doyle bd ft) per acre") + theme_nothing(legend = TRUE) + ggtitle('Saw volume heatmap predictions from IDW Model')