MultiSpec_PS.R

#Process multispectral orthomosaics

rm(list=ls())

#Required
library(sf)
library(raster)
#https://github.com/isciences/exactextractr/blob/master/README.md
library(exactextractr)
library(dplyr)


#### Set Paths ####
#Directory containing multispectral orthomosaics
Input_Directory <- 'C:/Users/lleyt/Desktop/Shapefile'

#Shapefile
Shapefile_Path <- 'C:/Users/lleyt/Desktop/Shapefile/MC_NB_GDA94.shp'


##crs EPSG code##
#https://epsg.io/
#EPSG = 28356


#Indices - True or False, if not required put F
NDVI = T
OSAVI  = T
NDRE = T
Thermal = T
GroundCover = T


####Start####
Mosaic_list <- list.files(path = Input_Directory, pattern = c("Altum",".tif"))

sf <- st_read(Shapefile_Path)

df <- NULL

for(i in 1:length(Mosaic_list)){
  start = Sys.time()
  out_df <- NULL
  out_df <- as_tibble(sf$Plot_ID)
  names(out_df) <- "Plot_ID"
  
  Flight_split <- unlist(strsplit(Mosaic_list[i], '_'))
  out_df <- out_df %>%
    mutate("Trial" = Flight_split[1],
           "Date" = as.Date(Flight_split[4]),
           "Time" = Flight_split[5],
           "Sensor" = Flight_split[6],
           "Height" = Flight_split[8])
  
  band_order <- Flight_split[7]
  
  #Read mosaic bands
  if(band_order == "bgreni"){
  blue <- raster(paste0(Input_Directory,'/',Mosaic_list[i]), band = 1)
  green <- raster(paste0(Input_Directory,'/',Mosaic_list[i]), band = 2)
  red <- raster(paste0(Input_Directory,'/',Mosaic_list[i]), band = 3)
  rededge <- raster(paste0(Input_Directory,'/',Mosaic_list[i]), band = 4)
  nir <- raster(paste0(Input_Directory,'/',Mosaic_list[i]), band = 5)
  thermal <- raster(paste0(Input_Directory,'/',Mosaic_list[i]), band = 6)
  #xb <- raster(paste0(Input_Directory,'/',Mosaic_list[i]), band = 7)
  }
  
  #Get crs EPSG code
  crs <- unlist(sf::st_crs(sf))
  split <- unlist(strsplit(crs[2], '"'))
  length_spl <- length(split)
  EPSG <- split[length_spl]
  EPSG <- sub(']', '', EPSG)
  EPSG <- sub(']', '', EPSG)
  EPSG <- sub(',', '', EPSG)
  
  #NDVI
  if(NDVI == TRUE){
    NDVI_raster = ((nir - red)/(nir + red))
    crs(NDVI_raster) <- paste0("EPSG:",EPSG)
    NDVI_raster <- as(NDVI_raster, "SpatRaster")
    
    out_df <- out_df %>%
      mutate(NDVI = exact_extract(NDVI_raster, sf, 'mean'),
             NDVIsd = exact_extract(NDVI_raster, sf, 'stdev'))
    print("NDVI complete")
  }
  
  #OSAVI
  if(OSAVI == TRUE){
    OSAVI_raster = ((nir - red)/(nir + red + 0.16))
    crs(OSAVI_raster) <- paste0("EPSG:",EPSG)
    OSAVI_raster <- as(OSAVI_raster, "SpatRaster")
    
    out_df <- out_df %>%
      mutate(OSAVI = exact_extract(OSAVI_raster, sf, 'mean'),
             OSAVIsd = exact_extract(OSAVI_raster, sf, 'stdev'))
    rm(OSAVI_raster)
    print("OSAVI complete")
  }
  
  #NDRE
  if(NDRE == TRUE){
    NDRE_raster = ((nir - rededge)/(nir + rededge))
    crs(NDRE_raster) <- paste0("EPSG:",EPSG)
    NDRE_raster <- as(NDRE_raster, "SpatRaster")
    
    out_df <- out_df %>%
      mutate(NDRE = exact_extract(NDRE_raster, sf, 'mean'),
             NDREsd = exact_extract(NDRE_raster, sf, 'stdev'))
    rm(NDRE_raster)
    print("NDRE complete")
  }
  
  #Thermal band
  if(Thermal == TRUE){
    #Convert to degrees Celcius
    Thermal_raster <- (thermal/100)-273.15
    crs(Thermal_raster) <- paste0("EPSG:",EPSG)
    Thermal_raster <- as(Thermal_raster, "SpatRaster")
    
    out_df <- out_df %>%
      mutate(Temperature = exact_extract(Thermal_raster, sf, 'mean'),
             Temperaturesd = exact_extract(Thermal_raster, sf, 'stdev'))
    
    
    #Temperature relative to plant pixels using NDVI threshold
    Threshold = 0.25
    Plant_raster <- NDVI_raster
    
    #Create mask
    Plant_raster[Plant_raster<=Threshold] <- NA
    Plant_raster[Plant_raster>Threshold] <- 1
    
    Thermal_masked <- mask(Thermal_raster, Plant_raster)
  
    out_df <- out_df %>%
      mutate(Temperature2 = exact_extract(Thermal_masked, sf, 'mean'),
             Temperature2sd = exact_extract(Thermal_masked, sf, 'stdev'))
    print("Thermal complete")
  }
  
  #Ground Coverage based on threshold
  if(GroundCover == T){
  plant <- exact_extract(Plant_raster, sf, fun = 'count')
  total <- exact_extract(blue, sf, fun = 'count')
  out_df <- out_df %>%
    mutate(Ground_Cover = plant/total)
  print("GroundCover complete")
  }
  
  df <- rbind(df, out_df)
  
  print(paste("Completed", i, "of", length(Mosaic_list), "Orthomosaics"))
  end = Sys.time()
  print(end-start)
}

#Create output directory
output_dir = paste0(Input_Directory,'/Output')
if(dir.exists(output_dir) == F){
  dir.create(output_dir)
}

#Write results to csv
write.csv(df, paste0(output_dir,'/Output.csv'), row.names = F)