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:
head(sawByPlot)
## 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
plotLocs@proj4string
## 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 nrow and 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 ggplot2 package.
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')
Finally, we can export the raster heatmap as a .tif file in order to read it into an external GIS system.
outRasterIDW <- rastSPDF
outRasterIDW <- raster(outRasterIDW)
writeRaster(outRasterIDW, "/tmp/heatmap_from_IDW.tif",
format = "GTiff", overwrite = TRUE)
That’s it! We have created a visual representation of the volume in the stand, using the plot data. The same basic steps can be followed for other variables, such as basal area, or for other products or individual species.
