Chapter 14 Thermal Conversion

This method converts raw imagery from the Zenmuse H20T to single band temperature tiffs. We use this method to bypass the DJI thermal analysis tool app since we had thousands of images to process and the app would often crash when a few images were uploaded. This method uses the DJI thermal analysis tool’s software developer kit (SDK) instead to loop through a folder of images and return temperature rasters.

This method consists of:

1. Gathering all “*_T.JPG” imagery into a single folder to loop through

2. Calculating weighted averaged values for humidity and ambident temperature, which will be needed for the thermal analysis SDK to run

3. Looping through imagery and converting to temperature rasters, then attaching original exif data to the new single band temperature rasters

14.1 Folder Set-Up

This step goes over how we organized our imagery for the thermal conversion.

Here we:

• create a folder named:

  • “Merged_T” within the original “H20T” folder, this is the location the raw thermal images will be copied. This allow us to (a) keep a backup of the original data and (b) easily iterate over the images to convert them in later steps

  • “Temperature_rasters” within the Merged_T folder that will hold converted imagery

  • “Temperature_rasters_EXIF” within the Merged_T folder that will hold converted imagery that has original imagery exif data written to it

• locate all folders within the defined directory that contain raw thermal images, in our case these were folders that had DJI and H20T in the name

• Iterate over the folders found in the above step and:

  • find files that end in _T.JPG (the raw thermal images)
  • copy over the raw thermal images to the merged folder

Note: Our folder structure was:

• paste0(“I:\PARSER_Ext\”,site,“\Flights\”, data_date, “\1_Data\H20T”), where “site” is a string of the site name that we defined prior to the for loop and data_date is a string representing the flight date that is defined in the for loop. This was done so we could easily iterate over sites and dates. Feel free to change this to match your folder structure.

Click to show the code

library(exiftoolr)
library(dplyr)
library(lubridate)
library(OpenImageR)
library(raster)
library(stringr)

# Setting the site and flight date to convert thermal imagery
flight_dates <-  c("2023_05_10") # Here only one date is shown however this can be a list of however many dates you would like to iterate over
site <- "site name here" 
merged <- "Merged_T" #name of the 'merged' folder that will contain all the h20T images from the multiple folders ouput by the H20T

for(i in 1:length(flight_dates)){
  data_date <- flight_dates[i]
  print(data_date)
  dir <- paste0("I:\\PARSER_Ext\\",site,"\\Flights\\", data_date, "\\1_Data\\H20T")
  #create folders
  if (!dir.exists(paste0(dir, "\\",merged))) {
    dir.create(paste0(dir, "\\",merged),recursive = TRUE)
  }
  if (!dir.exists(paste0(dir, "\\",merged,"\\Temperature_rasters"))) {
    dir.create(paste0(dir, "\\",merged,"\\Temperature_rasters"),recursive = TRUE)
  }
  if (!dir.exists(paste0(dir, "\\",merged,"\\Temperature_rasters_EXIF"))) {
    dir.create(paste0(dir, "\\",merged,"\\Temperature_rasters_EXIF"),recursive = TRUE)
  }
  
  list_dir <- list.files(dir, full.names = TRUE,pattern = c("DJI.+-H20T")) #selecting folders in the H20T folder that have DJI and end in -H20T, this is how we named our folders, change this pattern to match your folders. You should be selecting all folders with H20T imagery ending in _T for that flight
  print(list_dir) #to check only correct directories are being read
  
  for (d in 1:length(list_dir)){
    folder_dir <- list_dir[d] #calling each directly separately
    T_files <- list.files(folder_dir, pattern = "_T\\.JPG$") #selecting all JPEGs ending in _T 
    file.copy(from = paste0(folder_dir,"\\", T_files),
              to = paste0(dir,"\\",merged,"\\", T_files), overwrite = FALSE) #copying h20T images ending in _T (aka unprocessed thermal images) into the thermal only folder for ease in DJI SDK step
  }
}

14.2 Weather Data

DJI thermal’s SDK requires inputs of humidity, reflection (or ambient temperature), distance, and emissivity (see section on parameters for more detail on each of the inputs). In this section we will take hourly weather data from an on-site HOBO climate logger and calculate weighted averages for humidity and ambient temperature.

First, we use the first and last image in the Merged_T folder to define the start and end time of the flight, along with the date of the flight as a date object.

Click to show the code

#getting image capture date from exif data of the first image in the H20T folder

data_date <- flight_dates[1]

time_file_list <- list.files(paste0(dir,"\\",merged), pattern = "_T\\.JPG$")

#Retrieving the time the FIRST H20T thermal image was taken
(img_1 <- time_file_list[1]) #find the first thermal image and print the name
start_exif_flight_date <- data.frame(exif_read(paste0(dir,"\\",merged,"\\",img_1),
                                               tags = "CreateDate",
                                               quiet = TRUE)) #find the time the image was taken/written from the CreateDate tag in the exif data
(start_flight_date_time <- start_exif_flight_date$CreateDate) #Start of the flight

#Retrieving the time the LAST H20T thermal image was taken
(img_last <- tail(time_file_list, n=1))#End flight time
end_flight_date_time <- data.frame(exif_read(paste0(dir,"\\",merged,"\\",img_last),
                                             tags = "CreateDate",
                                             quiet = TRUE))
(end_flight_date_time <- end_flight_date_time$CreateDate)#End of the flight


(flight_date <- as.Date(start_exif_flight_date$CreateDate, '%Y:%m:%d %H:%M:%S')) #the flight date as a date object
(flight_date_filt <- as.Date(flight_date, format = '%Y-%m-%d')) #flight date in the format of 2024-03-10, YYYY-mm-dd, this ill be used in the next step to filter weather data

Next, since our weather station provided hourly measurements, we calculate a weighted average of temperature and humidity values based on time intervals. This was done to account for the fact that variables will change throughout the flight. The hourly variables are weighted in proportion to how much of the hour fell within the flight window which gives us a better estimate of the overall weather variables than if a plain average was taken.

To begin, we read in the weather data and set a start time and end time of the flight. Start time is rounded down to the nearest hour (i.e. 10:15am becomes 10:00am) and end time is rounded up to the nearest hour (i.e. 11:30am becomes 12:00pm). These values will be used to filter the weather data.

Click to show the code

cfs_weather <- "path to weather data" #loading in HOBO weather station data
cfs_weather$Date <- as.Date(cfs_weather$Date, format = "%m/%d/%y") #converting date in the format m/d/Y to a date object written as YYY-mm-dd

start_time_avg <- "10:00:00" #Based off of start_flight_date_time of 10:15am
stop_time_avg <- "12:00:00" #Based off of end_flight_date_time of 11:30am

For simplicity we kept weights as 0.25 (15min), 0.5 (30min), 0.75 (45min) or 1 (1 hr) and rounded the start time to the earliest 15 min marker and the end time to the latest 15 min marker. For example, for a flight that started at 10:03am and ended at 11:09am, we rounded 10:03am to 10:00am and 11:09am to 11:15am.

Weights for start, middle, and end hours:

• The start weight corresponds to how much of the first hour is covered by the flight time. For example, if the flight starts at 10:15, 45 minutes (or 0.75 of the hour) are included in that hour, so the start weight is 0.75.

• The middle weight (if the flight spans more than one hour) represents the full hour covered by the flight. In the example of a flight from 10:15 to 11:30, the middle weight would be 1 which takes into account that the drone was flying at 11am when the climate logger logged data.

• The end weight applies similarly for the last hour of the flight. If the flight ends at 11:30, 30 minutes (or 0.5 of the hour) are included in that hour, so the end weight is 0.5.

Click to show the code

#Ie: Flight starts at 10:15am and ends at 11:30am
start_weight <- 0.75 # 75% of 10am
mid_weight <- 1 # 100% of  11am 
end_weight <- 0.5 # 50% of 12pm

sum_weight <- sum(start_weight,
                  mid_weight,#comment out the mid_weight if no mid weight
                  end_weight) 

#creating df with weighted average air temperature and humidity values
(cfs_weather_time_avg <- cfs_weather%>%
  filter(Date == flight_date_filt ) %>%
  filter(Time <= stop_time_avg & Time >= start_time_avg)%>% #filtering for values taken between the defined start (10am) and end (12pm) times
  mutate(Average_hum = (start_weight*RH[1]
                         + mid_weight*RH[2] #comment out this line if no mid_weight and just below index of [3] to [2]
                         + end_weight*RH[3])/sum_weight,
         Average_air_temp = (start_weight*Temperature[1] 
                             + mid_weight*Temperature[2]#comment out this line if no mid_weight and just below index of [3] to [2]
                             + end_weight*Temperature[3])/sum_weight))

#check values 
(avg_hum <- cfs_weather_time_avg$Average_hum)
(avg_temp <- cfs_weather_time_avg$Average_air_temp)
(date <- as.Date(cfs_weather_time_avg$Date, '%d-%m-%Y'))

#creating a data frame  with weighted averaged values
Canoe_flight_table <- data.frame(date, avg_hum, avg_temp)
Canoe_flight_table <- Canoe_flight_table %>% distinct()#removing duplicate rows

#resetting index in table to go from 1-N
rownames(Canoe_flight_table) <- NULL

#setting variables for weighted average of humidity and ambient temperature, these variables will be used as inputs into DJI's thermal anlysis tool SDK to convert rasters to temperature
flight_humidity <- Canoe_flight_table[Canoe_flight_table$date == flight_date, "avg_hum"]
flight_humidity

flight_amb_temp <- Canoe_flight_table[Canoe_flight_table$date == flight_date, "avg_temp"]
flight_amb_temp

14.3 DJI Thermal SDK: temperature conversion

14.3.1 Parameters

The parameters required as inputs for DJI thermal’s SDK are:

• Reflected Temperature : The reflected temperature of the target takes into account the impact of radiation reflected from nearby objects that can significantly impact the temperature reading of the object. This parameter is important if there are any objects nearby with either an extremely high or low temperature. Otherwise the reflected temperature can be set to the ambient air temperature. Since none of our sites had an object with extreme high or low temperature nearby we used ambient air temperature (flight_amb_temp) from on-site HOBO climate loggers.

• Relative Humidity: The relative humidity during the flight. The default value is set to 70% and the allowed range is 20-100%. We used a weighted average humidity value (flight_humidity) calculated using data from the on-site HOBO climate loggers.

• Distance: The distance from the sensor to the target (i.e. tree crowns). The H20T is an infrared thermal sensor that measures the infrared radiation received from objects. Hence, the further away the object, the more the radiation attenuates and the less accurate the temperature measurement. The maximum distance that the thermal analysis tools allows is 25m. In our case, we are flying higher (~40m away from the crowns, depending on the site) and thus use the 25m maximum and interpret the temperature values as relative to each other while understanding that the absolute temperature likely contains error given the larger distance between the sensor and the object.

• Emissivity: The emissivity of a material is a measure of its ability to emit energy as thermal radiation. The default value is 1. We used a value of 0.98 as an estimate for all conifers (Rubio et al., 1997).

Click to show the code

#setting parameters:
(flight_amb_temp) # already set above, () to print out the value in the console to ensure it is the correct one
(flight_humidity) # already set above
distance <- 25 # max distance that can be set
emissivity <- 0.98 # approximate value for trees

14.3.2 Running the SDK

The SDK is run using a batch script that calls the SDK in a for loop to iterate over the imagery. To set up, first we copy the original .bat file into the Merged_T folder. We do this so that the original script does not get edited and instead each flight will have their own edited .bat file. This allows you to go back and see what parameters were used if needed.

The original .bat file can be found on GitHub and looks like this:

Click to show the code

REM to change dir to the location of this file
REM start time used to average humidity and temperature values: start_time_avg
REM end time used to average humidity and temperature values: stop_time_avg
REM weights (in 15min intervals ie 0.25, 0.5, 0.75 or 1, round down for start and up for end time): start w: start_weight, mid w: mid_weight (might not have one), end w: end_weight

cd /D "%~dp0"

mkdir DJI_SDK_raw

for %%i in (*.JPG) do (
echo Working on %%i...
C:\dji_thermal_sdk_v1.3_20220517\utility\bin\windows\release_x64\dji_irp.exe -s %%i -a measure -o DJI_SDK_raw/%%~ni.raw --humidity hum_change --distance dist_change --emissivity emis_change --reflection reflec_change
)
pause

This script makes a folder called DJI_SDK_raw in the directory the .bat is saved in (which will become the Merged_T folder once we copy it over in the below code) and converts the imagery in the Merged_T folder ending in .JPG to temperature using DJI’s thermal anlysis tool’s SDK.

Before continuing to the next steps:

1. save the .bat file from GitHub and update the omp_dir path in the below R code to the folder where this file was saved to

2. download the dji thermal sdk here

3. edit the location of the dji_irp.exe in your saved .bat file to match where the program is saved on your computer

   • i.e. replace the “C:\dji_thermal_sdk_v1.3_20220517\utility\bin\windows\release_x64\dji_irp.exe” in the above script to the directory you saved the dji_irp.exe file from step 2

Click to show the code

omp_dir <- "C:/Users/owaite/Documents/Scripts/DJI_SDK_TemperatureConversion_Script" #change to the directory where you saved the original .bat script from the above GitHub link to

file.copy(from = paste0(omp_dir,"\\", "DJI_SDK_loop.bat"), #copy the original omp.bat file into the Merged_T folder
          to = paste0(dir,"\\",merged,"\\","DJI_SDK_loop.bat"), overwrite = FALSE)

bat_old <- file(paste0(dir,"\\",merged,"\\","DJI_SDK_loop.bat")) #opening .bat file to edit
bat_new <- readLines(bat_old) #reading in lines of .bat file to edit and saving them out
close(bat_old)

#replacing values
comment out mid_weight line if necessary
bat_new <- gsub("hum_change",flight_humidity, bat_new)
bat_new <- gsub("dist_change",distance, bat_new)
bat_new <- gsub("emis_change",emissivity,bat_new)
bat_new <- gsub("reflec_change",flight_amb_temp,bat_new)
bat_new <- gsub("start_time_avg",start_time_avg, bat_new)
bat_new <- gsub("stop_time_avg",stop_time_avg,bat_new)
bat_new <- gsub("start_weight", start_weight, bat_new)
bat_new <- gsub("mid_weight", mid_weight, bat_new) #if no mid_weight, comment this line out
bat_new <- gsub("end_weight", end_weight, bat_new)

#Write the new file
fileConn <- file(paste0(dir,"\\",merged,"\\DJI_SDK_loop_updated.bat"))
writeLines(bat_new, fileConn)
close(fileConn)

Next, we run the updated .bat file in windows terminal using the R code below. This will automatically output a folder of new raw images with temperature in binary saved to the Merged_T/DJI_SDK_raw folder.

Click to show the code

shell.exec(paste0(dir,"\\",merged,"\\DJI_SDK_loop_updated.bat"))

14.4 Convert to Raster

The following R script:

• reads in the binary data outputted by the thermal SDK

• converts the data to rasters

• updates the new temperature raster with exif data from original _T.JPG H20T imagery

Click to show the code

dir_raw <- paste0(dir,"\\",merged,"\\DJI_SDK_raw\\")
raw_list <- list.files(dir_raw) #check this is a list of the images saved out from the SDK step above
start.time <- Sys.time()#to log how long it takes for code to run 

for (i in 1:length(raw_list)){#loop that only writes out a folder with the Temperature rasters
  wh <- c(640, 512)
  print("new")
  str_m <- raw_list[i]
  #getting name of .raw file but without the .raw extension so it can be saved as a tiff with the same name as the .raw (which has the same name as the original H20T images) so they can be swapped after alignment in metashape
  name <- substr(str_m,1,nchar(str_m)-4)
  print(name)
  #If else below allows you to rerun this code from where it left off if R aborts the session
  test_exists <- paste0(dir,"\\",merged,"\\Temperature_rasters_EXIF\\", name,".tiff")
  if (file.exists(test_exists)){
    print("exists")
  } else {
    #divide by 2 because 16 means 2 bit integers (?)
    #print(paste0(dir_raw, str_m))
    
    file.info(paste0(dir_raw, str_m))$size/2 
    prod(wh) #getting product of width and height and checking it is correct
    v <- readBin(paste0(dir_raw,str_m), what = "integer",  #getting numerical values
                 n = prod(wh), size = 2,  
                 signed = TRUE, endian = "little")
    
    v_temp <- v/10 #temp in deg Celsius, DJI requires we divide by 10 here for Celsius
    range(v_temp) #check the values
    matrix <- matrix(v_temp, wh[1], wh[2])[wh[1]:1,] # creating a matrix of the proper size
    mat_rotated = rotateFixed(matrix, 90)#rotate the matrix by 90 to get correct orientation of image, requires library OpenImageR
    #to view image
    #image(matrix(v_temp, wh[1], wh[2])[wh[1]:1,], useRaster = TRUE, col = grey.colors(256))
    rast <- raster(mat_rotated) #creating a raster 
   
    #matching original image to measure.raw file
    dir_h20t <- paste0(dir,"\\",merged,"\\") #original h20T images:
    str_h20t <- paste0(substr(str_m, 1, nchar(str_m)-3),"JPG") #since names are the same, here we call the name of the raw file -.raw and add .JPG
    r_h20t <- raster(paste0(dir_h20t,str_h20t)) #reading in jpg h20T raster
    r_h20t_path <- paste0(dir_h20t,str_h20t) #for later exiftool
    print(r_h20t_path) #checking the path
    
    #need to add these before exif because exif command used won't overwrite already written ones
    extent(rast) <- extent(r_h20t)# the extent is bound by the resolution that you have assigned to it - so must change extent before resolution
    res(rast) <- res(r_h20t)
    #crs(rast) <- crs(r_h20t)
    print(paste0(dir,"\\",merged,"\\Temperature_rasters\\", name, ".tiff"))
    writeRaster(rast, paste0(dir,"\\",merged,"\\Temperature_rasters\\", name,".tiff"), overwrite = TRUE) #writing out temperature raster, that does not have exif data attached
    r_temperature_path <- paste0(dir,"\\",merged,"\\Temperature_rasters\\", name,".tiff")
    
    #Use EXIFTOOL in terminal to add all remaining/missing exif data from the original H20T images to the new single band temperature TIFFS
    out_dir <- paste0(dir,"\\",merged,"\\Temperature_rasters_EXIF\\", name,".tiff")
    shell(paste0("exiftool -wm cg -tagsfromfile ",r_h20t_path," -all:all ",r_temperature_path, " -o ",out_dir)) #shell lets you run windows terminal cmds from R  
    }
  percent <- i/length(raw_list)*100
  print(paste0("percent complete: ", round(percent, 3))) #gives you the % of photos done to give an idea of code speed/progress
}
end.time <- Sys.time()
(time.taken <- round(end.time - start.time,2)) #time taken for the script to run

The output are single band thermal tiffs that can be loaded into metashape to create a thermal orthomosaic.

14.5 Examples

Below are two examples of data acquired using the H20T to show both the radiometric and spatial resolution of the H20T within the crowns as well as some important factors to keep in mind when interpreting crown temperatures.

Figure 14.1 shows the temperature and NIR values of two mature (~25 year old) Douglas-fir tree crowns, where the crowns delineate the upper 25th percentile of the tree, from a flight that occurred on July 27th, 2022 under heatdome conditions. These flights occurred in full sun conditions. Looking at the NIR plots, we can see shadowed areas (see chapter 10 for shadow masking with NIR) towards the north east sides of both crowns. These shadowed areas roughly align with the areas of cooler temperature seen in the crown temperature plots. As you can see in the daily temperature plot on the left, the MicaSense flight occurred from 1:30-2:30pm and was followed by the thermal flight, flown with the H20T, from roughly 2:40-3:40pm. We see a greater area of cooler crown temperatures in the thermal imagery compared to the NIR proxy for shadow likely due to the increase in shadow from the MicaSense flight to the thermal flight as the solar zenith angle decreases into the afternoon.

Figure 14.1: Left: plot of hourly ambient tempatures from an on-site HOBO climate logger at Canoe for July 27th, 2022. Middle column: near-infrared (NIR) values for two Douglas-fir crowns aquired from a MicaSense flight flown directly before the thermal flight. Right column: temperature values within two mature Douglas-fir crowns flown under heatdome conditions. Pixels plotted for both the NIR and temperature rasters are those within a delineated crown representing the top 25th percentile of the tree.

Figure 14.2 shows the temperature values for two young (~5 year old) Douglas-fir trees from an overcast day. Whereas we saw higher temperatures in areas with more direct light in Figure 14.1, we see the opposite pattern below. This is likely due to the heat of the ground radiating (i.e. warming) the lower branches of the young trees which are in view given the relatively small crown sizes and lack of crown closure. This is an important factor to keep in mind when analyzing crown temperatures.

Figure 14.2: Left: section of a P1 orthomosaic of W45 showing two neighbouring tree crowns. Right: temperature values for the two neighbouring young Douglas-fir.

References:

Rubio, E., Caselles, V., & Badenas, C. (1997). Emissivity measurements of several soils and vegetation types in the 8–14, μm Wave band: Analysis of two field methods. Remote Sensing of Environment, 59(3), 490–521. https://doi.org/10.1016/S0034-4257(96)00123-X