################################
### Get tree positions and heights from two surveys: 
###    LiDAR in 2010,  producing a DTM and a DSM, and
###    a drone (UAV) survey in 2016 producing a DSM.
### Calculate volume increase and where trees were harvested.
### 
### Modified from tutorial https://cran.r-project.org/web/packages/ForestTools/vignettes/treetopAnalysis.html
###           by Francesco Pirotti University of Padova - TESAF francesco.pirotti@unipd.it
################################

library(ForestTools)
library(rgdal)
library(raster)
library(rstudioapi)    

## Step 1: read the folder path to this script-file and set working directory to this folder
##         as my DTM/DSM files are in the same folder of the script
dn<-dirname(rstudioapi::getActiveDocumentContext()$path) 
setwd(dn)

## Step 2: load all DSMs and  DTMs from drone and from lidar - 
##         NB: if your raster files are in a different folder, or have a different name,
##             change the text inside the double quotes
dsm.lidar<-raster("lidar_dsm_32632.tif")
dtm.lidar<-raster("lidar_dtm_32632.tif")
dsm.drone<-raster("drone_dsm_32632.tif")
plot(dsm.lidar)


##################  Step 3:  CREATE THE CANOPY HEIGHT MODEL CHM - ALSO KNOWN AS THE NORMALIZED DIGITAL SURFACE MODEL (nDSM)
## apply the math: DSM - DTM to create the canopy height model - CHM, derived from LiDAR and from drone (UAV)
##         NB: we use the DTM from lidar because we assume that the terrain morphology does not change significantly
##             between 2010 and 2016.

###  LIDAR CHM = LIDAR DSM - LIDAR DTM
chm.lidar<-dsm.lidar - dtm.lidar  

###  DRONE CHM = DRONE DSM - LIDAR DTM
## force resolution of raster to 30 cm (DRONE DSM does not have square cells, which is required in the next functions)
res(dsm.drone)<- c(0.3, 0.3)
## NOT REQUIRED, only if dsm.lidar has a different CRS from dsm.drone:
##  - reproject the coordinate reference system of the dsm.drone to match the one from dsm.lidar
#  dsm.drone <- projectRaster(dsm.drone, crs = crs(dsm.lidar))

## apply math: DSM - DTM to create the canopy height model CHM
chm.drone<-dsm.drone - dtm.lidar  
### YOU GET ERROR!!! Because extents are different: how to fix it?

## resample and crop (2 function together) of the DSM drone to make sure it is identical to the resolution of the DTM lidar
## and allow math between rasters (subtraction)
dsm.drone.2m<- crop( resample(dsm.drone, dtm.lidar), dsm.drone)
## PS we could have done the opposite, match the extents and resolution of the dtm.lidar to the dsm.drone

## apply math: DSM - DTM to create the canopy height model CHM
chm.drone<-dsm.drone.2m - dtm.lidar 


## function for the moving window algorithm - x= CHM pixel value (height of tree canopy) -- see :
#    Popescu S.C., Wynne R.H., (2004) – Seeing the trees in the forest:Using lidar and multispectral data fusion with local filtering and variable window size for estimating tree height. Photogrammetric Engineering & Remote Sensing, 70: 589-604.
#    Pirotti, F. et al 2012. Laser Scanner Applications in Forest and Environmental Sciences. Italian Journal of Remote Sensing, 44, 109–123. doi:10.5721/ItJRS20124419)
lin <- function(x){x * 0.05 + 0.6}
## apply the function to extract tree positions from the CHM LiDAR only above 5 m!
ttops.lidar <- vwf(CHM = chm.lidar, winFun = lin, minHeight = 5)
## write a shapefile with tree positions so we can open in QGIS
shapefile(ttops.lidar,"treePosFromLidar", overwrite=TRUE) 
 
## get treetops again from the DRONE CHM
ttops.drone <- vwf(CHM = chm.drone, winFun = lin, minHeight = 5)
## write a shapefile with tree positions so we can open in QGIS
shapefile(ttops.drone,"treePosFromDrone", overwrite=T) 

## uncomment line below if you want to draw CHM
# plot(chm.drone)
## uncomment line below if you want to draw tree points
# plot(ttops, col = "blue", pch = 20, cex = 0.5, add = TRUE)

## uncomment line below if you want to get crown polygons and plot them 
# crowns <- mcws(treetops = ttops.lidar, CHM = chm.drone, minHeight = 2.5, verbose = FALSE)
# plot(crowns, col = sample(rainbow(50), length(unique(crowns[])), replace = TRUE) )

## get crown polygons and save them as shapefile polygons! 
crownsPoly<- mcws(treetops = ttops.lidar, CHM = chm.drone, format = "polygons", minHeight = 5, verbose = FALSE)
## write a shapefile with tree crown areas so we can open in QGIS
shapefile(crownsPoly,"treeCrownsFromDrone", overwrite=T) 

## uncomment line below if you want to get some statistics on your data and plot a couple of variables
#sp_summarise(crownsPoly, variables = c("crownArea", "height"))
#plot(crownsPoly$height, crownsPoly$crownArea)
