!!-Closed-!! PhD GRA Opportunity at Northern Arizona University in Flagstaff, AZ in the School of Forestry’s Quantitative Ecology Lab

 Update 2/12/2018 – The GRA has been awarded and I am no longer soliciting application materials.

Graduate Research Assistantship (PhD) Opportunity at Northern Arizona University in Flagstaff, AZ

Data Fusion for Forest Planning and Implementation: Ecological Restoration, Remote Sensing, and Data Analytics

Are you interested in a PhD program that will provide you an opportunity to work in the frequent fire forests of the American Southwest and influence ecological restoration practices? These forests are in dire need of restoration, mainly due to a century of fire exclusion and subsequent, undesirable changes in forest structure and function. For example, the largest collaborative forest restoration project in the US, the Four Forest Restoration Initiative (4FRI), has a goal of implementing restoration treatments on approximately 1M ha of U.S. Forest Service lands in northern Arizona. Fundamental to these efforts are precise data on the amount and distribution of available resources, knowledge of how resources may change over time, and hazard assessments (e.g., wildfire potential); all of which require costly and resource intensive, spatially explicit data. As a result, managers are using more remote sensing data products (e.g., LiDAR), coupled with advanced forest inventory and data analysis techniques, to quantify existing conditions and support broad-scale analysis of forest ecosystems.

A PhD graduate research assistantship is available in the School of Forestry at Northern Arizona University, Flagstaff, AZ, focused on the development and assessment of data fusion techniques that will allow managers to better capitalize on major advancements in remote sensing to utilize more accurate data and enhance precision of landscape-scale analysis (e.g., >100,000 acres) project areas. Working alongside the Ecological Restoration Institute, the USDA Forest Service, USDI Fish and Wildlife Service, The Nature Conservancy, and Campbell Global; the successful applicant will focus on developing and statistically validating an open source big data, remote sensing, and inventory data fusion platform. This platform will provide enhanced forest structural and compositional information in support of forest resource decision-making.

The selected student will:

  1. Assess and statistically validate algorithms for identifying individual trees and species from remote sensing data of Southwestern forests using new and/or existing stemmapped, area and tree based sample data.
  2. Using these algorithms and data, design and implement a platform that integrates multiple data sources (data fusion) that are typically too large to analyze using traditional methods (big data) to provide detailed forest resource information at the tree-,stand-, and landscape levels.
  3. Assess the accuracy, precision, and statistical properties of forest resource estimates such as bias, consistency, error, spatial uncertainty, and use these to provide improved information for land management decision making.
  4. Apply the platform to Southwestern landscape-scale case studies to; quantify existing conditions, assess low-value biomass product availability, facilitate watershed treatment implementation, and monitor forest restoration treatments.

The position includes a full stipend, tuition waiver, health benefits and field support for 4 years.

Applications from quantitatively minded individuals with a practical approach to solving complex problems are welcome. Experience processing large remote sensing and inventory datasets using C++, R, and/or Python is preferred.


  • Master’s degree in forestry, geography, ecology, computer science, or related fields.
  • Demonstrable research experience, collaboration abilities, and English (written and oral) communication skills.
  • Competitive GRE scores (top 40th percentile).

 Information about NAU’s graduate program, including eligibility requirements, is available at http://nau.edu/CEFNS/Forestry/Degrees/.

NAU’s formal application deadline is for Fall 2018 is Feb 15 2018 and preferred start date is Summer 2018. However, interested candidates are encouraged to contact with Dr. Sanchez Meador as soon as possible using the information provided below or submit your CV, written statement of interest, and copies of unofficial degree transcripts to initiate a dialog via e-mail.  Andrew.SanchezMeador@nau.edu.

Contact Information:

Dr. Andrew Sánchez Meador 
School of Forestry 
Northern Arizona University
Flagstaff, AZ 86011-5018, USA 

Link to announcement PDF


Using rLiDAR and FUSION to delineate individual trees through canopy height model segmentation

example3Let’s face it. Light detection and ranging (LiDAR) is awesome and it’s completely arrived. Currently, the best sources for free nationwide LiDAR datasets are the United States Interagency Elevation Inventory, USGS Center for LIDAR Information Coordination and Knowledge, and NSF’s OpenTopography. Just as there are quite a few sources for datasets, your choices are equally diverse when it comes to tools for processing LiDAR data and for the detection and delineation of individual trees from LiDAR (see Jakubowski et al 2013 for a nice review). Personally, I use a combination of FUSION, Global Mapper’s LiDAR Module, LiForest’s implementation of Li et al’s 2012 point cloud segmentation method, Swetnam and Falk’s 2014 variable area local maxima algorithm (implemented in MatLab), and the local maximum with a fixed window size algorithm implemented in rLiDAR by Carlos Alberto Silva*.

The following is a worked example using a forested LiDAR dataset, FUSION to build a CHM (called form R) and processing the resulting output using rLiDAR*. You need to install FUSION from here, and I also assume it’s install in “C:/FUSION/”. I hope you find the script useful…

# First things first, we need to set up some directories to keep the raw data 
# separate from that produced by the analysis. I basically have an input (.las) 
# and output (.las, .dtm, etc) directory in a dropbox folder.     
mainDir <- "C:/Dropbox/LiDAR_Prj"
inDir <- "inputDirectory"
outDir<- "outputDirectory"
dir.create(file.path(mainDir, outDir), showWarnings = FALSE)

# First, we process the data in FUSION

# Read in the .las file and use FUSION to produce a .las file of points that 
# approximate the ground's surface (bare-earth points). 
# http://forsys.cfr.washington.edu/fusion/FUSION_manual.pdf#page=94&zoom=auto,70,720
system(paste(file.path("C:","Fusion", "groundfilter.exe"),
       "/gparam:0 /wparam:1 /tolerance:1 /iterations:10",
       file.path(mainDir, outDir, "Example_GroundPts.las"),
       file.path(mainDir, inDir, "Example.las"),
       sep=" "))

# Next we use GridSurfaceCreate to compute the elevation of each grid cell using the 
# average elevation of all points within the cell. Check the manual for arguments and uasge 
# http://forsys.cfr.washington.edu/fusion/FUSION_manual.pdf#page=88&zoom=auto,70,720 
system(paste(file.path("C:","Fusion", "gridsurfacecreate.exe"),
       file.path(mainDir, outDir, "Example_DEM.dtm"),
       "1 M M 1 12 2 2",
       file.path(mainDir, outDir, "Example_GroundPts.las"),
       sep=" "))
# Next we use CanopyModel to create a canopy surface model using a LIDAR point cloud. 
# By default, the algorithm used by CanopyModel assigns the elevation of the highest return within 
# each grid cell to the grid cell center.	   
system(paste(file.path("C:","Fusion", "canopymodel.exe"),
       paste("/ground:",file.path(mainDir, outDir, "Example_DEM.dtm"), sep=""),
       file.path(mainDir, outDir, "Example_CHM.dtm"),
       "1 M M 1 12 2 2",
       file.path(mainDir, inDir, "Example.las"),
       sep=" "))

# Lastly, we use DTM2ASCII to convert the data stored in the PLANS DTM format into ASCII raster
# an file. Such files can be imported into GIS software such as ArcGIS or QGIS.
# http://forsys.cfr.washington.edu/fusion/FUSION_manual.pdf#page=88&zoom=auto,70,720
system(paste(file.path("C:","Fusion", "dtm2ascii.exe"),
       file.path(mainDir, outDir, "Example_CHM.dtm"),
       file.path(mainDir, outDir, "Example_CHM.asc"),
       sep=" "))

# Second, we process the resulting CHM in rLiDAR

#install.packages("rLiDAR", type="source")
#install.packages("raster", dependencies = TRUE)

# Import the LiDAR-derived CHM file that we just made in the above section and plot it
chm<-raster(file.path(mainDir, outDir, "Example_CHM.asc"))

Here’s the resulting canopy height model

# If we want, rLiDAR can smooth the CHM using a Gaussian filter
# Set the window size (ws)
ws<-3 # dimension 3x3
# Set the filter type
# filter<-"mean"
# Set the sigma value for the Gaussian filter
sCHM<-CHMsmoothing(chm, filter, ws, sigma) 

# Edited - Have a look at the comments section in this post to see the discussion that triggered this fix... basically, before we can call FindTreesCHM(), it needs to be
# replaced with an edited, custom version. Why? - well as of 04/06/2016 the FindTreesCHM funciton has a couple of bugs - namely one form a call to SpatialPoints() with specifying the projection
# and one in the call to colnames() (really, it's cbind(), but who's counting) 
FindTreesCHM<-function(chm, fws = 5, minht = 1.37) 
    if (class(chm)[1] != "RasterLayer") {
        chm <- raster(chm)
    if (class(fws) != "numeric") {
        stop("The fws parameter is invalid. It is not a numeric input")
    if (class(minht) != "numeric") {
        stop("The minht parameter is invalid. It is not a numeric input")
    w <- matrix(c(rep(1, fws * fws)), nrow = fws, ncol = fws)
    chm[chm < minht] <- NA
    f <- function(chm) max(chm)
    rlocalmax <- focal(chm, fun = f, w = w, pad = TRUE, padValue = NA)
    setNull <- chm == rlocalmax
    XYmax <- SpatialPoints(xyFromCell(setNull, Which(setNull == 
        1, cells = TRUE)), proj4string = crs(chm))                # Edited
    htExtract <- over(XYmax, as(chm, "SpatialGridDataFrame"))
    treeList <- cbind(coordinates(XYmax), htExtract)              # Edited
    colnames(treeList) <- c("x", "y", "height")

# Setting the fixed window size (fws)
fws<-3 # dimention 3x3
# Set the specified height above ground for the detection break
# Create the individual tree detection list and summarize it
loc<-FindTreesCHM(sCHM, fws, minht)

# Set the maxcrown parameter - maximum individual tree crown radius expected
# Set the exclusion parameter - A single value from 0 to 1 that represents the % of pixel
# exclusion. E.g. a value of 0.5 will exclude all of the pixels for a single tree that has 
#a height value of less than 50% of the maximum height from the same tree. Default value is 0.3.
# Compute canopy areas for the individual tree detections
canopy<-ForestCAS(sCHM, loc, maxcrown, exclusion)

# Retrieve the boundary for individual tree detection and canopy area calculation

# Retrieve the list of individual trees detected for canopy area calculation
canopyList<-canopy[[2]] # list of ground-projected areas of individual tree canopies
summary(canopyList)     # summary
canopyList$crad<-sqrt(canopyList$ca/pi) # Compute the corresponding crown radii

# Write the output to a CSV file. This will make bringing it into ArcGIS or QGIS by others easier
write.csv(canopyList, file.path(mainDir, outDir, "Example_Out.csv"))

# Plot the results in ggplot

# Let's convert the tree results into spatial points to be used in ggplot
XY<-SpatialPoints(canopyList[,1:2])    # Spatial points
XY<-data.frame(XY) # Converted to a dataframe
# We're not dealing with lat and long, so let's rename the columns

# Rasters are a little problematic, so convert the values into rows in a dataframe with 
# the corresponding information
CHMdf <- rasterToPoints(sCHM); CHMdf <- data.frame(CHMdf)
colnames(CHMdf) <- c("X","Y","Ht")
# Build the breaks for plotting
b.chm <- seq(0,50,10)

# Plotting the individual tree canopy boundary over the CHM
ggplot(CHMdf) +
  geom_raster(data=CHMdf,aes(X,Y,fill=Ht)) +
  scale_fill_gradientn(name="Canopy Height",colours = terrain.colors(length(b.chm))[length(b.chm):1],breaks=b.chm) +
  geom_polygon(data = fortify(boundaryTrees), aes(x=long, y=lat,
                                                  group = group),colour='black', fill='transparent')+
  geom_point(data=XY, aes(x=x, y=y), color="black", shape=3, size=0.5)+
  scale_shape_discrete(name = "Tree Locations", labels=c("Tree Centroid","")) + 
  scale_alpha(range = c(0, 0.5)) +
  ggtitle("Location map showing individual trees and associated crown areas \n segmented from a LiDAR-derived canopy height model") +
  coord_equal() + theme_bw()

You can download the script from here and here’s the resulting figure.

Rplot01* – As per Carlos suggestion in the comments, here’s how one might cite rLiDAR – “Silva, C.A., Crookston, N.L., Hudak, A.T., and Vierling, L.A. 2015. rLiDAR: An R package for reading, processing
and visualizing LiDAR (Light Detection and Ranging) data, version 0.1, accessed Oct. 15 2015.”



More Fun with LiDAR

LiDARPlotView Just a little animated GIF showing LiDAR from the North Kaibab. Notice the penetration achieved and the resolution it provides for picking out large trees, small bushes, and topography.

Fun with LiDAR

So I finally got my hands on the North Kaibab LiDAR acquisition and it’s huge (~2TB). I was able to find a few minutes to play with it (in ENVI) to produce this simple figure of the landscape surrounding Matt Tuten’s research sites. Cool beans…