This R-Markdown document is a demonstration of loading, correcting, extracting and analysing thermal images collected at Holden Arboretum

To edit and extend these analyses, open ‘process_thermal_images.Rmd’.

I’ve worked up an example of how to work with thermal images. My code is rough and will certainly contain many mistakes. Please use this as your starting point and improve it! If you find mistakes please tell me :)

I’ve annotated my code with comments to make it easier to follow, but refer to the Thermimage website for more details and examples.

To run the code below you’ll need to install and load the ‘Thermimage’ and ‘fields’ packages.

# read in micromet data for correction
met <- read.csv('../data/micromet/Tower_top_071818_interim.csv') # this is from Chris Chambers on the 18th, newer data should be available from Meter

# calculate RH from Vapor Pressure
met$es <- 0.6108 * exp(17.27 * met$Ta_C/(met$Ta_C + 237.3)) # saturation vapor pressure
met$RH <- (met$VapP_kPa / met$es) * 100

# format 'DateTime'
met$DateTime <- as.POSIXct(met$DateTime, format = '%m/%d/%Y %H:%M')

# interpolate Met data to 1 minute intervals
new_times <- seq(min(met$DateTime), max(met$DateTime), 60)
new_data <- data.frame(new_times, NA, NA)
names(new_data) <- c('DateTime', 'Ta', 'RH')

new_data$Ta <- spline(met$DateTime, met$Ta_C, xout = new_data$DateTime)$y
new_data$RH <- spline(met$DateTime, met$RH, xout = new_data$DateTime)$y


# read in TIR images
files <- dir('../data/images_edit/', pattern = '\\.seq$', full.names = T)
img <- list()

# to save time with this demonstration, I'm only opening every 6th image (every 30 minutes)
saver <- seq(1, 486, 6)

# load thermal images
for (i in 1:81){
    img[[i]] <- readflirJPG(files[saver[i]], exiftoolpath="installed")
}; print('Images imported')
## [1] "Images imported"
# extract thermal parameters from thermal image metadata. The 'flirsettings' function requires exiftool to be installed. Refer to the thermimage website for more info. Place 'exiftool.exe' in your /Rmd/ directory.
cams <- list()
dateOriginal <- NULL; PlanckR1 <- NULL; PlanckB <- NULL;
PlanckF<- NULL; PlanckO <- NULL; PlanckR2 <- NULL
ReflT <- NULL; AtmosT <- NULL

for (i in 1:81){
  tmp <- flirsettings(files[saver[i]], exiftoolpath="installed", camvals="")
  dateOriginal[i] <- tmp$Dates$DateTimeOriginal   # Original date/time extracted from file
  PlanckR1[i] <-    tmp$Info$PlanckR1         # Planck R1 constant for camera  
  PlanckB[i] <-     tmp$Info$PlanckB          # Planck B constant for camera  
  PlanckF[i] <-     tmp$Info$PlanckF          # Planck F constant for camera
  PlanckO[i] <-     tmp$Info$PlanckO          # Planck O constant for camera
  PlanckR2[i] <-    tmp$Info$PlanckR2         # Planck R2 constant for camera
  ReflT[i] <-       tmp$Info$ReflectedApparentTemperature  # Reflected apparent temperature
  AtmosT[i] <-      tmp$Info$AtmosphericTemperature        # Atmospheric temperature
}; print('Parameters extracted')
## [1] "Parameters extracted"
# get dates and times from images and format
img_times <- strptime(dateOriginal, format = '%Y-%m-%d %H:%M') # wrong timezone, but not a problem

img_times <- img_times + 3*60*60 # add 3 hours to image times to correct timezone

# subset met data to match image capture times
img_met <- subset(new_data, new_data$DateTime %in% img_times)
# correct thermal data using micromet data

# we're using the 'raw2temp' function from Thermimage. Go look at the webpage for more information: https://github.com/gtatters/Thermimage

# unfortunately, we don't have downwelling longwave radiation as we didn't have a pyrgeometer on the micromet station. it's possible to estimate it with a complicated model, but for simplicity here I'm going to use air temperature.

temperature_corr <- list()
  
for (i in 1:length(img)){
  temperature_corr[[i]] <- raw2temp(img[[i]], 
            E = 0.97, OD = 1.2, 
            RTemp = new_data$Ta[i],
            20, new_data$Ta[1], 
            1, 
            RH = new_data$RH[i], PlanckR1[i], PlanckB[i], PlanckF[i], 
            PlanckO[i], PlanckR2[i])
}

# for some reason the images are sideways - rotate them using the code below
rotate <- function(x) t(apply(x, 2, rev))
temperature_corr <- lapply(temperature_corr, rotate)

# temperature<-raw2temp(img, ObjectEmissivity, OD, ReflT, AtmosT, IRWinT, IRWinTran, RH,
#                       PlanckR1, PlanckB, PlanckF, PlanckO, PlanckR2)
# visualize an image for fun
# I'm going to use a small ROI to demonstrate how to extract and analyse data. Note that the camera moved during the time it was on the tower, so any ROI's you use will need to be modified whenever the camera moves. I've looked through all the images and identified when I think the camera moved - see 'holden_thermal_info.txt'


plotTherm(temperature_corr[[1]], h=240, w=320, minrangeset=min(temperature_corr[[1]]), maxrangeset=max(temperature_corr[[1]]))

test <- temperature_corr[[1]][278:307,191:214] # extract ROI using [columns, rows]

plotTherm(test, h = dim(test)[2], w = dim(test)[1], minrangeset=min(temperature_corr[[1]]), maxrangeset=max(temperature_corr[[1]]))

# You can see that this ROI includes some branches (very hot). To improve data quality so that it is only 'leaf temperature' you might need to figure out a way to filter out branches.
# extract ROIs in each image. A more elegant way to do this might be to convert the list into a rasterstack.

ext_temp1 <- lapply(temperature_corr, function(x){y <- x[278:307,191:214]})
ext_temp2 <- lapply(temperature_corr, function(x){y <- x[38:77,44:68]})

leaf_temp1 <- unlist(lapply(ext_temp1, mean))
leaf_temp2 <- unlist(lapply(ext_temp2, mean))
img_met$Tl1 <- leaf_temp1
img_met$Tl2 <- leaf_temp2
# plot some data!

{
  plot(img_met$Tl1 ~ img_met$DateTime, type = 'l', xlab = 'Time', ylab = expression(Temperature~(degree~C)),
       main = 'Canopy and Air Temperatures')
  lines(img_met$Tl2 ~ img_met$DateTime, lty = 2)
  lines(img_met$Ta ~ img_met$DateTime, lty = 2, col = 'red')
  legend("topright", c('ROI1', 'ROI2', 'Air'), lty = c(1,2,2), col = c('black', 'black', 'red'))
}

plot(img_met$RH ~ img_met$DateTime, type = 'l', xlab = 'Time', ylab = 'RH (%)', main = 'Relative Humidity')