Early Detection of Myrtle Rust on Pōhutukawa Using Indices Derived from Hyperspectral and Thermal Imagery

Abstract

Myrtle rust is a very damaging disease, caused by the fungus Austropuccinia psidii, which has recently arrived in New Zealand and threatens the iconic tree species pōhutukawa (Metrosideros excelsa). Canopy-level hyperspectral and thermal images were taken repeatedly within a controlled environment, from 49 inoculated (MR treatment) and 26 uninoculated (control treatment) pōhutukawa plants. Measurements were taken prior to inoculation and six times post-inoculation over a 14-day period. Using indices extracted from these data, the objectives were to (i) identify the key thermal and narrow-band hyperspectral indices (NBHIs) associated with the pre-visual and early expression of myrtle rust and (ii) develop a classification model to detect the disease. The number of symptomatic plants increased rapidly from three plants at 3 days after inoculation (DAI) to all 49 MR plants at 8 DAI. NBHIs were most effective for pre-visual and early disease detection from 3 to 6 DAI, while thermal indices were more effective for detection of disease following symptom expression from 7 to 14 DAI. Using results compiled from an independent test dataset, model performance using the best thermal indices and NBHIs was excellent from 3 DAI to 6 DAI (F1 score 0.81–0.85; accuracy 73–80%) and outstanding from 7 to 14 DAI (F1 score 0.92–0.93; accuracy 89–91%).

1. Introduction

The extent and impact of invasive fungal pathogens on the world’s forest ecosystems is increasing, in part driven by growing connectivity and human-assisted spread into new areas [1]. Commonly known as myrtle rust, the disease caused by the fungus Austropuc–cinia psidii [2], basionym Puccinia psidii, affects plants in the Myrtaceae family. Native to South and Central America [3,4], invasive strains of the fungus have invaded new areas over the past two decades, having been reported in North America, parts of Asia, Australia, New Caledonia, South Africa and New Zealand [5]. Combined with the pathogen’s wide host range within the Myrtaceae [6,7], the spread of A. psidii has been facilitated by the long-distance dispersal of urediniospores either via wind or anthropogenic movement of spores and infected plant material. This pathogen is a major threat to Myrtaceae in native forest ecosystems as well as myrtle-related industries around the world [8,9,10,11,12].

Due to its impact, and the threat of invasion into new areas by strains currently restricted in range [13], biosecurity restrictions on the movement of Myrtaceae and infected plants are recommended. To minimise the risk of spreading the pathogen, anyone working with or managing Myrtaceae needs to be able to reliably identify infected plants. Field diagnosis of A. psidii requires identification of the host and symptoms. The disease develops on young actively growing tissues while older tissues develop ontogenic resistance to infection [14]. The first visible symptoms appear as distortion, blistering and small purple/red coloured lesions or spots. Eventually, uredinia or telia erupt containing characteristic powdery spores. Identification after visible symptoms have progressed to the stage of uredinia, containing bright yellow spores, is straightforward. However, it is difficult to confirm a diagnosis during the very early stages before uredinia erupt, and not possible during the latent period before symptoms develop. The length of the latent period varies depending on the host and environmental conditions. Under ideal conditions, early symptoms typically occur within 4–7 days and uredinia erupt 5–10 days after infection [14].

The diagnosis of infected Myrtaceae while asymptomatic or with early-stage symptoms is challenging for plant producers. This is particularly the case with large volumes of stock where close visual inspection is unlikely to be practical or cost-effective and may be biased or inaccurate. Diagnosis can be confirmed by examining spores under the microscope or using DNA-based diagnostic assays that can detect the fungus, even in symptomless tissue [15,16,17,18,19]. However, while these methods are easy to apply when sampling is guided by visual detection of symptoms, a random sampling from asymptomatic plants would likely be ineffective and expensive. Hyperspectral and thermal imagery provides a promising method of pre-visual and early disease detection that could be spatially scaled for rapid screening. Physiological responses in plants subject to infection by pathogens often include changes in chlorophyll and other pigments, leaf water content, and key photosynthetic variables such as stomatal conductance, transpiration and photosynthesis that may be detected by shifts in leaf temperature and changes in indices or plant functional traits [20]. Pre-visual or early detection of a range of diseases, in agricultural and forestry crops, have been demonstrated using hyperspectral and thermal imagery at the plant level, often at broad scales [21,22,23,24,25,26].

A recent study [27] found that indices derived from hyperspectral data and thermal imagery were effective under controlled conditions for the pre-visual and early detection of myrtle rust on the model host rose apple (Syzygium jambos (L.) Alston). Excellent pre-visual classification (F1 score range = 0.89–0.94) was achieved from 3 to 1 Days Before Symptoms (DBSs) using narrow-band hyperspectral indices (NBHIs) derived from leaf-level data measured on older resistant leaves using a spectroradiometer. Models using indices derived from thermal imagery could perfectly (F1 score = 1.0) distinguish control plants from myrtle rust-inoculated plants one DBSs were visually detected and through the entire early symptom period. Indices that were most important in models for pre-visual and early detection characterized changes in chlorophyll, water stress and a reduction in leaf temperature, which were associated with measured increases in transpiration. However, the recently developed Lemon Myrtle–Myrtle Rust (LMMR) spectral index [10] for the classification of lemon myrtle (Backhousia citriodora F. Muell. 1859) leaves with advanced disease symptoms was not a significant predictor of pre-visual or early symptoms on rose apple. These findings reinforce the need to test the utility of hyperspectral and thermal indices for the detection of myrtle rust on a range of Myrtaceae and at different stages of disease development. Further research should use canopy-level, rather than leaf-level, hyperspectral data, as this matches the scale of imagery acquisitions when these models are applied for detection purposes within a nursery.

Pōhutukawa (Metrosideros excelsa Sol. ex Gaertn.) is an endemic tree growing up to 25 m, predominantly found in coastal vegetation in northern Aotearoa New Zealand [28]. Of great historical and cultural significance to Māori [29,30], pōhutukawa is also revered as the ‘New Zealand Christmas Tree’ due to the proliferation of red flowers produced in December-January [28]. Planted outside of its natural range in urban areas, pōhutukawa is a native myrtle commonly stocked and propagated in nurseries throughout New Zealand [31]. Having undergone a significant decline following human settlement, the remaining natural pōhutukawa stands are threatened by invasive species and land use change [28]. The most recent threat to pōhutukawa is the pandemic strain of A. psidii, first detected on the New Zealand mainland in 2017, on pōhutukawa at a nursery in the North Island, that has since spread rapidly through its predicted climatic range [32,33,34]. In New Zealand, the current recorded host range for A. psidii (including native, exotic and hybrid taxa) comprises 36 Myrtaceae [35], of which pōhutukawa is one of the most frequently reported hosts [33].

This study examines the accuracy of indices derived from thermal and visible near-infrared (VNIR) hyperspectral imagery to detect myrtle rust on pōhutukawa plants, growing in a controlled environment. The objectives of this research were to (i) identify the key thermal and narrow-band hyperspectral indices (NBHIs) associated with pre-visual and early expression of myrtle rust on pōhutukawa and (ii) develop a classification model to detect the disease.

2. Materials and Methods

2.1. Experimental Setup

The experiment was undertaken within a dark room in which environmental conditions were regulated from late November to mid December 2023 at Scion (Rotorua, New Zealand). A total of 80 healthy pōhutukawa, that did not show any visual signs of infection, were selected for the experiment. These plants, which were sourced from a Rotorua nursery, were grown in 2.5 litre pots and averaged 35 cm in height at the initiation of the experiment. The plants were well watered by hand over the course of the experiment.

The plants were randomly allocated to a control and a myrtle rust (MR) treatment group on the 29th November. These two groups were separated into different plastic enclosures that had misting capability. Over the measurement period, the following conditions were maintained within each enclosure: 16 h photoperiod under 30 W LED 3000 K grow lights, 22 °C temperature, and relative humidity of 70–80%. The plants remained under the enclosures for most of the time and were only briefly removed for measurements. Within the room, the environmental conditions were similar to those of the enclosure and relatively constant, with an average air temperature of 22 °C and relative humidity of 61%.

During the measurement period, four of the control plants developed minor symptoms and were immediately removed from the enclosure. One of the plants allocated to the MR treatment developed symptoms of another disease and was also removed from the enclosure. Following these removals there were a total of 75 plants used for analyses, of which 26 plants were allocated to the control and 49 plants were allocated to the MR treatment. No further MR symptoms were observed on remaining control plants for 7 days after measurements ended, so they were considered disease-free.

Following the method described in [27] the plants in the MR treatment were inoculated on the 4th December 2023. Inoculum used in the experiment was mass-produced on rose apple and harvested using a portable vacuum pump (Mini Cyclone Spore Collector, Tallgrass Solutions, Manhattan, KS, USA) into 00 gelatin capsules and stored for 24 h at low humidity over silica beads. Prior to inoculation, urediniospores were suspended in a solution comprising 0.05% Tween 20 in sterile distilled water, achieving a concentration of 1 × 10⁵ spores/mL, as quantified using a haemocytometer. The inoculum was uniformly dispersed across the surfaces of young leaves and stems utilizing a gravity-fed airbrush connected to a compressor, operating at a pressure of 2 bar. Subsequent to inoculation, plants were relocated to the enclosures, where they were kept in darkness at a temperature of 18 °C and a relative humidity exceeding 90% for a duration of 24 h to foster optimal infection conditions. Control plants received a mock inoculation employing a sterile 0.05% Tween 20 solution, with all environmental conditions mirrored to eliminate any potential physiological discrepancies attributable to environmental alterations between treatment groups. Following this process, conditions in the enclosures for plants in both the control and MR treatment were returned to a 16 h photoperiod, 22 °C and a relative humidity of 70–80% for the remaining duration of the experiment.

Pre-inoculation measurements were taken from both treatments on the 29th November. Following inoculation, daily measurements were taken from 3–8 days after inoculation (DAI), with a final measurement made 14 DAI on the 18th December. Of the 49 MR plants used within the analyses, the number of symptomatic plants increased rapidly from 3 plants (at 3 DAI) to 12 (4 DAI), 47 (5 DAI), 48 (6 DAI), and 48 (7 DAI), with all 49 plants showing symptoms by 8 DAI.

2.2. Visual Assessment of Symptoms

Following the method described in [27], we recorded the presence or absence of signs and symptoms for each new leaf and stem, that was distinguishable from mature growth in colour and texture. These symptoms included slight bumps or blistering on the leaf surface, red/purple spots, yellow discolouration, and developing pustules yet to break through the epidermis or fully emerged pustules with spores erupted on the surface.

2.3. Thermal Measurements

Canopy-level thermal data were captured using an infrared thermal camera (FLIR A700SC, Teledyne FLIR LLC, Wilsonville, OR, USA). This instrument is characterized by a spectral sensitivity spanning 7.5 µm to 14 µm and can detect thermal gradients as subtle as 30 mK. The imaging resolution of the camera is 640 by 480 pixels. The thermal camera was mounted on a 2 m high tripod, directly above the plant, with a nadir orientation. FLIR Research Studio (Version 3.2.0) was used to control the camera (Figure 1).

Imagery was acquired in a 32-bit TIFF format. To enhance the precision of each data capture, air temperature, relative humidity, emissivity, and target distance were input into the Research Studio software, version 3.2.0. The target distance was set as 2 m and following [26], the emissivity was set to 0.98. The ambient air temperature and relative humidity for the surrounding environment was determined using an HMP155A probe (Vaisala, Vantaa, Finland) connected to a CR1000 datalogger (Campbell Scientific, Logan, UT, USA) that took continuous measurements of these variables at a 1 min interval.

Two calibration targets were positioned either side of the plant within the thermal camera field of view (Figure 1, shown immediately to the left of the bottom right plant). One of these targets consisted of a container with damp soil while black cotton cloth was stretched across the top of the other container. During the capture, the temperatures of these calibration targets were measured using the thermal camera and a handheld IR thermometer. Linear regression was used to develop a relationship between the temperature recorded by the camera (x-variable) and thermometer (y-variable) using data pooled by capture day. This regression was used during post-processing to calibrate the measurements of canopy temperature made by the camera.

The canopy of each individual plant was delineated using object-based image segmentation algorithms incorporated within the Fiji module of ImageJ software, version 1.53t [36]. For every image captured, the average temperature across the segmented crown area was determined. These canopy temperatures were then used as the independent variable (x variable) in the previously described calibration equation to predict the calibrated canopy temperature, denoted as T_c. These T_c values were then paired with the ambient air temperature, T_a, that was concurrently recorded with each measurement and used to determine, for each tree, the normalised canopy temperature, as T_c − T_a.

2.4. Hyperspectral Measurements

2.4.1. Data Acquisition

Plant-level hyperspectral data were acquired on all measurement dates under artificial lighting across the visible near-infrared (VNIR) range using a FX10 hyperspectral camera (Specim, Spectral Imaging Ltd., Oulu, Finland). This push-broom camera captures 448 bands across the 400–1000 nm range with a full width half maximum of 5.5 nm and spectral sampling interval of 2.7 nm. The FX10 has a high signal-to-noise ratio of 600:1 and the spectral sampling includes 1024 pixels captured within a 38° field of view.

Images were recorded in nadir view and the camera was mounted 2 m above the ground on a cross beam (Figure 1). A conveyor belt was used to move the plants through the field of view. The speed of the conveyor belt was adjusted to fit the frame rate of the camera, which in turn was dependent on the exposure time. During the measurements, the conveyor belt speed and frame rate were kept constant and the exposure time was adjusted to avoid over or undersaturation. A diffuse white reference standard (Spectralon, North Hutton, NH, USA) was placed so that it was visible in every frame, allowing calibration of the imagery as a function of the changing illumination conditions.

The conveyor was located inside the dark room and six ASD illuminator reflectance lamps (Malvern Panalytical technologies, Worcestershire, UK) mounted on the frame were used to light the scene (Figure 1). Each of these lamps is a quartz–tungsten–halogen light source with an integrated reflector and in combination they provided stable illumination across the 350–2500 nm range.

2.4.2. Processing of Data and Extraction of Narrowband Hyperspectral Indices

All pre-processing of the hyperspectral data was carried out using Matlab (The MathWorks, Inc., Natick, MA, USA) following the methods described in [37]. The white reference bar was placed close to the edge of the field of view in the hyperspectral images. The brightest pixels in each image line belonged to the white reference bar. These were used to convert digital numbers to reflectance factors by dividing all values by the mean of the five brightest pixels and multiplying the result by the white reference reflectance. Due to wavelength shifts across the field of view, the reflectance signal was distorted. An additional hyperspectral image of a white reference plate filling the complete field of view using the same geometric set-up was recorded. Dividing this image by its own values at the same column numbers where the white reference bar had been detected in the original image resulted in a correction spectrum that could be used to remove the distortions from the reflectance signals [37].

In the corrected reflectance images, pixels containing vegetation were identified using two simple rules: Pixels with NDVI ≥ 0.4 and reflectance at 780 nm ≥ 0.1 were selected as vegetation pixels. The resulting masks were checked visually and coincided well with the plant material. The mean spectrum for all vegetation pixels was extracted and used for further analyses. Following this masking, the plant level spectra were interpolated to a 1 nm resolution from 397–1004 nm. Narrow-band hyperspectral indices (NBHIs) that have been linked to the xanthophyll pigment, red/green/blue bands, plant disease, water content, curvature index and plant structure were derived from the VNIR spectral region for use in analyses. Equations and references for these 81 NBHIs are given in

Table A1. Thermal and narrow-band hyperspectral indices (NBHIs) that were used in analyses.

Indices	Equation	Ref.
Thermal indices
Normalised canopy temperature	Tc − Ta	[66]
Standard deviation normalised temp.	TSD = std dev (Tc − Ta)	[66]
Xanthophyll indices
Photochemical Refl. Index (570)	${PRI}_{570}={(R}_{570}-R_{531})/(R_{570}+R_{531})$	[67]
Photochemical Refl. Index (515)	${PRI}_{515}={(R}_{515}-R_{531})/(R_{515}+R_{531})$	[68]
Photochemical Refl. Index (528)	${PRI}_{528}={(R}_{528}-R_{567})/(R_{528}+R_{567})$	[67]
Photochemical Refl. Index (550)	${PRI}_{550}={(R}_{550}-R_{531})/(R_{550}+R_{531})$	[67]
Photochemical Refl. Index m1	${PRI}_{m1}={(R}_{512}-R_{531})/(R_{512}+R_{531})$	[68]
Photochemical Refl. Index m2	${PRI}_{m2}={(R}_{600}-R_{531})/(R_{600}+R_{531})$	[67]
Photochemical Refl. Index m3	${PRI}_{m3}={(R}_{670}-R_{531})/(R_{670}+R_{531})$	[67]
Photochemical Refl. Index m4	${PRI}_{m4}={(R}_{570}-R_{531}-R_{670})/{(R}_{570}+R_{531}+R_{670})$	[68]
Normalized Photoch. Refl. Index	${PRI}_n={PRI}_{570}/[RDVI·{(R}_{700}/R_{670})]$	[69]
Ratio of PRI to Simple Ratio	$DRI PRI={PRI}_{570}/(R_{800}/R_{670})$	[70]
Carotenoid/Chlorophyll Ratio Index	$PRI·CI={(R}_{570}-R_{530})/{(R}_{570}+R_{530})·\left(\right(R_{760}/R_{700})-1)$	[71]
R/G/B indices
Redness Index	$R=R_{700}/R_{670}$	[72]
Greenness Index	$G=R_{570}/R_{670}$	[47]
Greenness Index 2	$GI=R_{554}/R_{677}$	[73]
Blue Index	$B=R_{450}/R_{490}$	[47]
Blue/green indices	$BGI1=R_{400}/R_{550}$	[73]
	$BGI2=R_{450}/R_{550}$	[73]
Blue/red indices	$BRI1=R_{400}/R_{690}$	[74]
	$BRI2=R_{450}/R_{690}$	[74]
BF1	$BF1=R_{400}/R_{410}$	[25]
BF2	$BF2=R_{400}/R_{420}$	[25]
BF3	$BF3=R_{400}/R_{430}$	[25]
BF4	$BF4=R_{400}/R_{440}$	[25]
BF5	$BF5=R_{400}/R_{450}$	[25]
Red/green index	$RGI=R_{690}/R_{550}$	[73]
Ratio Analysis of Reflectance Spectra	$RARS=R_{746}/R_{513}$	[75]
Lichtenthaler indices	$LIC1=\left(R_{800}-R_{680}\right)/\left(R_{800}+R_{680}\right)$	[76]
	$LIC2=R_{440}/R_{690}$	[76]
	$LIC3=R_{440}/R_{740}$	[76]
Plant disease indices
Cercospora leaf spot index	$CLS=\frac{R_{698}-R_{570}}{R_{698}+R_{570}}-R_{734}$	[22]
Healthy-index	$HI=\frac{R_{534}-R_{698}}{R_{534}+R_{698}}-\frac{1}{2}·R_{704}$	[22]
Powdery mildew index	$PMI=\frac{R_{520}-R_{584}}{R_{520}+R_{584}}+R_{724}$	[22]
Sugar beet rust–index	$SBRI=\frac{R_{570}-R_{513}}{R_{570}+R_{513}}+\frac{1}{2}·R_{704}$	[22]
Water Indices
Floating position Water Band index	${fWBI=R}_{900}/{min}_{R930-R980}$	[77]
Water Band Index	${WBI=R}_{970}/R_{900}$	[78]
Water Index	${WI=R}_{900}/R_{970}$	[79]
Curvature index
Curvature index	$CUR=\left(R_{675} R_{690}\right)/{R_{683}}^2$	[80]
Structural indices
Normalized Difference Veg. Index	$NDVI={(R}_{800}-R_{670})/(R_{800}+R_{670})$	[81]
Renormalized Difference Veg. Index	$RDVI={(R}_{800}-R_{670})/\sqrt{(R_{800}+R_{670})}$	[82]
Optimized Soil-Adjusted Veg. Index	$OSAVI=(\left(1+0.16\right)·(R_{800}-R_{670})/(R_{800}+R_{670}+0.16\left)\right)$	[83]
Modified Soil-Adjusted Vegetation Index	$MSAVI=\frac{2·R_{800}+1-\sqrt{{(2·R_{800}+1)}^2-8(R_{800}-R_{670})}}{2}$	[84]
Triangular Vegetation Index	$TVI=0.5·[120·\left(R_{750}-R_{550}\right)-200·\left(R_{670}-R_{550}\right)]$	[85]
Modified Triangular Veg. Index 1	$MTVI1=1.2[1.2\left(R_{800}-R_{550}\right)-2.5\left(R_{670}-R_{550}\right)]$	[86]
Modified Triangular Veg. Index 2	$MTVI2=\frac{1.5\left[1.2\left(R_{800}-R_{550}\right)-2.5\left(R_{670}-R_{550}\right)\right]}{\sqrt{{(2R_{800}+1)}^2-\left(6R_{800}-5\sqrt{R_{670}}\right)-0.5}}$	[86]
Chlorophyll Abs. Reflectance Index	$CARI=\left(R_{700}-R_{670}\right)-0.2\left(R_{700}-R_{550}\right)$	[87]
Modified Chlorophyll Abs. Index	$MCARI=\left[\left(R_{700}-R_{670}\right)-0.2\left(R_{700}-R_{550}\right)\right]·(R_{700}/R_{670})$	[86]
Modified Chlorophyll Abs. Index 1	$MCARI1=1.2\left[2.5\left(R_{800}-R_{670}\right)-1.3\left(R_{800}-R_{550}\right)\right]$	[86]
Modified Chlorophyll Abs. Index 2	$MCARI2=\frac{1.5\left[2.5\left(R_{800}-R_{670}\right)-1.3\left(R_{800}-R_{550}\right)\right]}{\sqrt{{(2R_{800}+1)}^2-\left(6R_{800}-5\sqrt{R_{670}}\right)-0.5}}$	[86]
Modified Chlorophyll Abs. Index 3	$MCARI3=\left[\left(R_{750}-R_{705}\right)-0.2\left(R_{750}-R_{550}\right)\right]·(R_{750}/R_{705})$	[88]
Simple Ratio	$SR=R_{800}/R_{670}$	[89]
Modified Simple Ratio	$MSR=\frac{R_{800}/R_{670}-1}{{(R_{800}/R_{670})}^{0.5}+1}$	[90]
Enhanced Vegetation Index	$EVI=2.5·\left(R_{800}-R_{670}\right)/\left(R_{800}+6·R_{670}-7.5·R_{800}+1\right)$	[91]
Pigment indices
Vogelmann indices	$VOG1=R_{740}/R_{720}$	[92]
	$VOG2={(R}_{734}-R_{747})/{(R}_{715}+R_{726})$	[92]
	$VOG3={(R}_{734}-R_{747})/{(R}_{715}+R_{720})$	[92]
Gitelson & Merzlyak indices	$GM1=R_{750}/R_{550}$	[53]
	$GM2=R_{750}/R_{700}$	[53]
	$GM4=R_{750}/R_{555}$	[93]
Transformed Chlorophyll Absorption in Reflectance Index	$TCARI=3·[\left(R_{700}-R_{670}\right)-0.2·\left(R_{700}-R_{550}\right)\left(\frac{R_{700}}{R_{670}}\right)]$	[94]
TCARI/OSAVI	$\frac{TCARI}{OSAVI}$	[94]
Chlorophyll Index Red Edge	$CI=R_{750}/R_{710}$	[94]
Simple Ratio Pigment Index	$SRPI=R_{430}/R_{680}$	[51,95]
Normalized Phaeophytinization Index	$NPQI={(R}_{415}-R_{435})/{(R}_{415}+R_{435})$	[51,95]
Normalized Pigments Index	$NPCI={(R}_{680}-R_{430})/{(R}_{680}+R_{430})$	[95]
Carter indices	$CTR1=R_{695}/R_{420}$	[96]
	$CAR=R_{695}/R_{760}$	[97]
Reflectance band ratio indices	$DCabCxc=R_{672}/{(R}_{550}·R_{708})$	[98]
	$DNIRCabCxc=R_{860}/{(R}_{550}·R_{708})$	[98]
Structure-Insensitive Pigment Index	$SIPI=(R_{800}-R_{445})/(R_{800}-R_{680})$	[95]
Carotenoid Reflectance Indices	${CRI}_{550}=\left(1/R_{510}\right)-\left(1/R_{550}\right)$	[99,100]
	${CRI}_{700}=\left(1/R_{510}\right)-\left(1/R_{700}\right)$	[99]
	${CRI}_{550 515}=\left(1/R_{515}\right)-\left(1/R_{550}\right)$	[99]
	${CRI}_{700 515}=\left(1/R_{515}\right)-\left(1/R_{700}\right)$	[99]
	${RNIR·CRI}_{550}=(\left(1/R_{510}\right)-\left(1/R_{550}\right))·R_{770}$	[99,100]
	${RNIR·CRI}_{700}=(\left(1/R_{510}\right)-\left(1/R_{700}\right))·R_{770}$	[99,100]
Plant Senescing Reflectance Index	$PSRI={(R}_{680}-R_{500})/R_{750}$	[101]
Pigment Specific Simple Ratio Chlorophyll a	$PSSRa=R_{800}/R_{675}$	[102]
Pigment Spec. Simple Ratio Chl. b	$PSSRb=R_{800}/R_{650}$	[102]
Pigment Specific Simple Ratio Carotenoid	$PSSRc=R_{800}/R_{500}$	[102]
Pigment Specific Normalized Difference	$PSNDc=(R_{800}-R_{470})/(R_{800}+R_{470})$	[102]
Reciprocal reflectance	$RR=1/R_{700}$	[103]

.

2.5. Physiological Measurements

Physiological measurements were taken prior to inoculation on the 28th November and from symptom-free leaves, after inoculation on the 13th December at 9 DAI. A subset of the total plants was measured on both dates over the course of one day, with pre-inoculation and post-inoculation measurements comprising, respectively, 22 and 26 Control plants and 28 and 31 MR plants. Photosynthetic parameters, including net photosynthetic rate (A), stomatal conductance (g_s), and transpiration rate (E), were quantified using the GFS-3000 gas exchange analyser (M-Series, Heinz Walz GmbH, Effeltrich, Germany) following a pre-illumination period of 120 s, under an ambient carbon dioxide concentration of 400 ppm. The GFS-3000 gas exchange analyser was connected to an Imaging-PAM chlorophyll fluorometer (Heinz Walz GmbH, Effeltrich, Germany) and was equipped with a supplementary CO₂ cartridge to ensure stable CO₂ levels during the measurements.

2.6. Data Analysis

2.6.1. Treatment Differences in Measured Variables

All analyses described in this section were undertaken using R version 4.2.3. [38]. Reflectance was plotted against wavelength, by treatment, for all captures, and a t-test was used to identify the wavelengths with significant treatment differences for each capture. Differences between treatments for all 83 hyperspectral and thermal indices were examined using a t-test for all post-inoculation captures and p values for all these tests were documented. Boxplots showing treatment variation against DAI for the variables that most frequently differed between treatments, or were most important in the classification model (see Section 2.6.2), were constructed. A t-test was used to determine if there were significant treatment differences for each of the physiological variables during both measurement dates.

2.6.2. Classification Model

Random forests was used to classify the two treatments (Control, MR plants) from the measured thermal and hyperspectral indices using scikit learn, version 0.23.2. [39] which was implemented in python. Random Forests is an ensemble learning method that creates a collection of decision trees during training. Each tree in the forest is built from a sample drawn with replacement (bootstrap sample) from the training set, and split decisions are made at each node based on a random subset of the features. This approach enhances diversity among the trees, which reduces overfitting and improves generalisation to unseen data. The final classification decision is made based on the majority voting principle, where the class label predicted by the majority of the trees is chosen as the final output, thereby enhancing the accuracy and robustness of the model.

Separate random forests models were created using data for each of the seven post-inoculation measurement dates, with features comprising either (i) the two thermal indices, (ii) the 81 NBHIs or (iii) the thermal indices and NBHIs. Thus, there were 21 classification models created that included the factorial combination of the seven measurement dates and three sets of predictor variables.

Prior to model fitting, recursive feature elimination (RFE) was used with the random forest algorithm to subset the predictors to the most important variables for each of the 14 models that included either NBHIs or thermal indices/NBHIs. The RFE undertook a stratified cross validation and was constrained so that no more than 10 features, that had correlations of <0.9, were selected to avoid redundancy and overfitting in the final models. Both thermal indices were included in the seven thermal models as there was little collinearity between these variables for all DAI (R² ≤ 0.06). The final variables included in the models are shown in

Table 1. Variables that were included in the classification models by days after inoculation (DAI). Models were created using either thermal indices, narrow-band hyperspectral indices (NBHIs) or a combination of both types (NBHIs + Thermal). The indices are sorted in order of importance for each model and the indices that appear most frequently in the NBHI and NBHI + thermal models are bolded. Full names, equations and references for all indices are presented in Table A1.

Data	DAI	Variables in the Model
Thermal	3, 8	Tc − Ta
indices	All others	Tc − Ta, TSD
NBHIs	3	PRI528, CUR, PRI CI, RARS, PRIm1, VOG3
	4	HI, BF1, BGI1, B, fWBI
	5	B, PRI528, NPQI, R, DRI PRI
	6	CUR, PRI528, RGI, RR
	7	PRI528, NPQI, CUR, WBI, B, CTR1, PRI CI, BF4
	8	PRI528, HI, RDVI, CRI700 515, RR, WI, R, SIPI
	14	DRI PRI, PRIn, EVI, BF2, RR, BF1
NBHIs +	3	PRI CI, PRI528, CUR
Thermal	4	HI, BF1, BGI1, fWBI, PRI528, RGI, R
indices	5	NPQI, PRI528, STD, B, DRI PRI, RR, PRI CI, Tc − Ta
	6	PRI528, CUR, Tc − Ta
	7	Tc − Ta, PRI528, B, CUR, BF4
	8	Tc − Ta
	14	Tc − Ta, DRI PRI, PRIn, EVI, BF5

. The number of variables for models with NBHIs and NBHIs/thermal indices, ranged from, respectively, 4–8 (mean = 6) and 1–8 (mean = 5).

Following feature selection for each model, a stratified split (with respect to treatment) was used to divide observations into a training dataset comprising 80% of the observations (60 data points) with the remaining 20% (15 observations) retained as a test dataset. The random forest model was fitted to the training dataset using a stratified cross validation with predictions from this fitted model made on the independent test dataset. This process was repeated 49 times using a different train/test split during each iteration, and performance statistics were averaged over predictions on all 50 test datasets. This approach removed any potential bias that may have resulted from an unbalanced choice of test dataset using a single train/test split providing a more rigorous evaluation of true model performance.

Three hyperparameter grids were evaluated for model fitting that varied in complexity, with the simplest using default parameters and the most complex optimising model performance across a range in five parameters (number of estimators, max depth, max features, min samples split, min samples leaf). Decisions around which of these three hyperparameter grids to use for each of the 21 models were based on mean model performance extracted from the 50 iterations on the test dataset to ensure generality. For almost all of the 21 models the default hyperparameter grid resulted in the best model performance.

Using the mean predictions made on all 50 independent test datasets, a confusion matrix was constructed for each model. As is common practice in experiments where disease is predicted [40], the MR plants were designated positive and values in the confusion matrix quantified the percentage of true positives (TPs), true negatives (TNs), false positives (FPs) and false negatives (FNs).

As the dataset was unbalanced, model performance was predominantly assessed using precision, recall and the F1 score (for equations see [27]). Precision quantifies the fraction of correct positive predictions whereas recall identifies the fraction of true positives that were accurately detected. High values of precision mean there are few false positives while high values of recall indicate low numbers of false negatives. The F1 score represents the harmonic mean of precision and recall, serving as a balanced measure of the model discriminatory power. F1 scores span a continuum from 0 to 1, with categorizations from 0.5 to 0.7, 0.7 to 0.8, 0.8 to 0.9, and >0.9 signifying poor, acceptable, excellent, and outstanding levels of discrimination, respectively.

Accuracy is defined as the ratio of correctly predicted observations to the total observations examined. While accuracy is not as useful as the F1 score for evaluating model performance using unbalanced datasets, this metric was incorporated into the analysis to provide a comprehensive assessment.

2.6.3. Relationships between Physiological Variables and Indices

Using the post-inoculation data, linear models were constructed between the three physiological variables (g_s, E and A) and the 83 hyperspectral and thermal indices listed in Table A1 (i.e., 249 models comprising 3 dependent variables × 83 predictor variables). Model strength was assessed for all relationships that were significant using the coefficient of determination (R²).

2.7. Software Used

The software that was used within the analyses is shown in

Table 2. Software and modules used within analyses.

Software	Modules	Methods Sections
Matlab version 2022a	None	2.4.2.
R version 4.2.3	ggplot2, dplyr, tidyverse, broom, gridExtra	2.6.1., 2.6.3.
Python 3.8.5.	pandas, numpy, sklearn	2.6.2.

.

3. Results

3.1. Disease Symptoms

The first symptoms occurred 3 DAI and by 8 DAI all plants were symptomatic. Between 3 and 5 DAI only bumps were observed. These bumps typically became pronounced and sometimes developed yellow chlorotic colouration between 4 and 5 DAI. There were no red spots or lesions observed at any stage during the experiment. Pustules were developing between 4 and 5 DAI and open from 6 DAI onwards. The extent and severity of symptoms varied widely between plants, and the proportion of new leaves on each plant with symptoms typically increased over time, as did the area with symptoms on leaves and stems (Figure 2).

3.2. Hyperspectral Spectra

Very little treatment variation was noted in plant level reflectance across the wavelength range for pre-treatment or post-treatment measurements at 3 and 4 DAI. For measurements taken at 5, 7 and 8 DAI, there was slight treatment divergence in reflectance above the upper end of the red edge at ca. 740 nm with values in the MR treatment exceeding those in the control (Figure 3). These differences in mean reflectance diverged markedly over this spectral range at 14 DAI with values for plants in the MR treatment exceeding those of control plants most markedly at 912 nm. Apart from significant treatment differences at 4 DAI that occurred at 404 nm, significant differences between the treatments were not noted until 14 DAI (Figure 3). At 14 DAI, highly significant differences were noted and this significance generally increased with wavelength, reaching p = 0.05 by 718 nm and p = 0.001 by 839 nm (Figure 3).

3.3. Variation in Indices between Treatments

3.3.1. Thermal Indices

There were no significant differences between the two thermal indices during the pre-treatment measurements or those taken at 3 DAI (Figure 4). However, by 4 DAI values of T_c − T_a for plants in the MR treatment were significantly lower than those of the control by 0.16 °C (−1.17 vs. −1.01 °C). This divergence in T_c − T_a increased over the remaining measurements with differences becoming highly significant by 6 DAI (Figure 4). Treatment differences in T_c − T_a were highly significant from 7–14 DAI (p range of 1.5 × 10⁻¹⁰–4.9 × 10⁻¹²), during which time they exceeded 0.49 °C and over this period these treatment differences were most significant among all 83 examined indices (

Table A2. One-way analysis of variance between treatments for hyperspectral and thermal indices used in the study over the eight infection stages. The p values are colour-coded according to the strength of the significance. Red shaded cells indicate a highly significant effect, orange cells indicate marginal significance, and yellow and green shaded cells indicate moderate to strong insignificance, respectively. Indices are sorted according to their type and follow the order shown in Table A1.

Index Type	Variable	Days after Inoculation
		Pre-Treat	3	4	5	6	7	8	14
Thermal indices	Tc − Ta	0.9339	0.9973	0.0151	0.0202	0.0001	8.0 × 10−10	1.5 × 10−10	4.9 × 10−12
	TSD	0.7341	0.7462	0.7547	0.1267	0.5007	0.0438	0.7740	0.4121
Xanthophyll	PRI570	0.1696	0.3332	0.8786	0.4486	0.6389	0.9258	0.6662	0.0001
indices	PRI515	0.6713	0.9850	0.7761	0.9359	0.8190	0.7972	0.8824	0.9736
	PRI528	0.1164	0.0105	0.0249	0.0161	0.0045	0.0159	0.2869	0.0632
	PRI550	0.7438	0.8201	0.8678	0.7875	0.5777	0.7670	0.4621	0.0119
	PRIm1	0.5321	0.9617	0.8583	0.8523	0.8947	0.6965	0.6662	0.7227
	PRIm2	0.8027	0.5339	0.5010	0.9919	0.9516	0.7548	0.7028	0.0016
	PRIm3	0.4172	0.9760	0.4142	0.7622	0.7877	0.9519	0.8988	0.1006
	PRIm4	0.2592	0.6061	0.4235	0.5254	0.9872	0.9792	0.8863	0.4564
	PRIn	0.1495	0.3196	0.2400	0.5044	0.6155	0.9495	0.6615	4.2 × 10−5
	DRI PRI	0.1120	0.1363	0.0659	0.2706	0.3367	0.9111	0.3130	3.3 × 10−5
	PRI CI	0.2810	0.4728	0.8524	0.3062	0.7850	0.9060	0.6340	0.0252
R/G/B Indices	R	0.2691	0.6516	0.0221	0.5573	0.7884	0.6937	0.4057	0.3709
	G	0.3597	0.9302	0.4649	0.6874	0.8223	0.9647	0.9066	0.4573
	GI	0.4279	0.9213	0.8493	0.7929	0.5656	0.7646	0.7260	0.3613
	B	0.3526	0.4589	0.0027	0.0156	0.5287	0.2028	0.6033	0.0084
	BGI1	0.8036	0.5767	0.0114	0.4078	0.8960	0.4431	0.7473	0.0561
	BGI2	0.6555	0.7916	0.2865	0.4779	0.9907	0.8035	0.8828	0.0849
	BRI1	0.8712	0.7311	0.0330	0.4418	0.5525	0.6060	0.8989	0.0245
	BRI2	0.6972	0.9057	0.9400	0.4729	0.3743	0.8033	0.2886	0.0078
	BF1	0.9467	0.1314	0.0001	0.3741	0.4419	0.1247	0.2108	0.0128
	BF2	0.9729	0.4207	0.0010	0.3413	0.8676	0.1249	0.6432	0.0309
	BF3	0.9967	0.3856	0.0018	0.3317	0.5423	0.2856	0.7265	0.0588
	BF4	0.8910	0.4000	0.0020	0.2861	0.6901	0.2605	0.5967	0.0458
	BF5	0.9584	0.5672	0.0019	0.5733	0.6533	0.3933	0.6790	0.1121
	RGI	0.7144	0.6174	0.1288	0.7800	0.2476	0.5355	0.3252	0.4593
	RARS	0.6215	0.8175	0.1998	0.7759	0.5264	0.6026	0.5225	0.8905
	LIC1	0.9935	0.8217	0.3202	0.9349	0.6021	0.6564	0.9279	0.8995
	LIC2	0.6516	0.6119	0.8630	0.9425	0.3314	0.5478	0.2442	0.0339
	LIC3	0.8452	0.9957	0.9865	0.8679	0.7182	0.8787	0.5647	0.2464
Plant disease	CLS	0.9084	0.5916	0.1020	0.2650	0.5584	0.3057	0.1284	0.0038
indices	HI	0.2138	0.3411	0.0033	0.5509	0.1879	0.3993	0.1174	0.1151
	PMI	0.9385	0.5872	0.6885	0.1753	0.6217	0.3406	0.2644	0.0961
	SBRI	0.3296	0.6373	0.6180	0.7500	0.5532	0.6695	0.7972	0.1207
Water indices	fWBI	0.6527	0.2858	0.0027	0.6518	0.3417	0.9034	0.1109	0.1139
	WBI	0.9572	0.3024	0.0041	0.2341	0.6796	0.0348	0.1291	0.4479
	WI	0.9567	0.2952	0.0039	0.2300	0.6722	0.0340	0.1273	0.4415
Curvature index	CUR	0.1417	0.8552	0.4686	0.8227	0.0479	0.1769	0.4135	0.2775
Index type	Variable	Days after inoculation
		Pre-treat	3	4	5	6	7	8	14
Structural	NDVI	0.9698	0.8664	0.1813	0.9312	0.9171	0.9239	0.8942	0.7090
Indices	RDVI	0.9671	0.7938	0.6754	0.2942	0.8220	0.4784	0.2238	0.0020
	OSAVI	0.9256	0.7527	0.2883	0.4071	0.8518	0.6110	0.4268	0.0307
	MSAVI	0.8501	0.7599	0.3139	0.4204	0.9937	0.7376	0.5517	0.0291
	TVI	0.8454	0.8632	0.9636	0.3120	0.9335	0.5332	0.2683	0.0031
	MTVI1	0.7978	0.8614	0.9223	0.3258	0.8473	0.4990	0.2281	0.0020
	MTVI2	0.5411	0.8222	0.2969	0.4758	0.9762	0.7379	0.4683	0.0342
	CARI	0.2392	0.6730	0.4107	0.9809	0.7583	0.8778	0.9167	0.4113
	MCARI	0.2241	0.6262	0.1932	0.8480	0.7474	0.8398	0.7488	0.7026
	MCARI1	0.7978	0.8614	0.9223	0.3258	0.8473	0.4990	0.2281	0.0020
	MCARI2	0.5411	0.8222	0.2969	0.4758	0.9762	0.7379	0.4683	0.0342
	MCARI3	0.4349	0.8673	0.8378	0.5370	0.9843	0.7256	0.4463	0.1168
	SR	0.9745	0.9173	0.1201	0.7887	0.8755	0.8529	0.7028	0.5369
	MSR	0.9737	0.9753	0.1325	0.8257	0.9306	0.9124	0.7517	0.5790
	EVI	0.8383	0.7887	0.8961	0.3333	0.9987	0.6844	0.3790	0.0030
Pigment indices	VOG1	0.5178	0.6816	0.5840	0.7091	0.4973	0.5356	0.7814	0.9248
	VOG2	0.5644	0.4978	0.3041	0.5141	0.3558	0.3667	0.7560	0.8244
	VOG3	0.5705	0.5047	0.3254	0.5310	0.3678	0.3828	0.7456	0.8214
	GM1	0.5232	0.8457	0.3251	0.8766	0.7383	0.8431	0.6565	0.9276
	GM2	0.3531	0.8788	0.9008	0.9224	0.9619	0.9567	0.9495	0.8588
	GM4	0.5139	0.8536	0.3366	0.8726	0.7511	0.8240	0.6782	0.9259
	TCARI	0.2916	0.7881	0.6687	0.9170	0.8386	0.9381	0.9517	0.3950
	TCARI/OSAVI	0.2856	0.7514	0.7335	0.9951	0.8127	0.8762	0.9781	0.5465
	CI1	0.4346	0.8919	0.8206	0.9602	0.7672	0.8840	0.9508	0.9723
	SRPI	0.5624	0.9834	0.0301	0.5735	0.4467	0.8662	0.6757	0.0166
	NPQI	0.5400	0.7747	0.0682	0.5994	0.2253	0.1706	0.9208	0.6934
	NPCI	0.5349	0.9904	0.0308	0.6467	0.3982	0.8710	0.6789	0.0193
	CTR1	0.3898	0.8225	0.6454	0.9615	0.3569	0.7955	0.5274	0.0895
	CAR	0.2471	0.6085	0.7877	0.7402	0.7532	0.7448	0.7507	0.9543
	DCabCxc	0.4016	0.8245	0.9726	0.6721	0.7013	0.7063	0.5185	0.2347
	DNIRCabCxc	0.6054	0.6545	0.3816	0.5901	0.6013	0.6187	0.4315	0.2525
	SIPI	0.4444	0.8275	0.0383	0.4006	0.3902	0.9822	0.6538	0.0054
	CRI550	0.9078	0.9725	0.5531	0.6046	0.9243	0.7201	0.2542	0.2296
	CRI700	0.9602	0.8092	0.1635	0.4966	0.7319	0.5870	0.1478	0.2726
	CRI550 515	0.7129	0.8853	0.3813	0.5831	0.6905	0.5869	0.3172	0.1489
	CRI700 515	0.9113	0.6749	0.0726	0.4598	0.5230	0.4612	0.1597	0.2026
	RNIR CRI550	0.9417	0.8506	0.3818	0.9991	0.9523	0.9022	0.4103	0.4427
	RNIR CRI700	0.8965	0.8554	0.0674	0.7858	0.7020	0.7004	0.2222	0.4876
	PSRI	0.3947	0.8987	0.9464	0.9815	0.0838	0.2245	0.0703	0.0101
	PSSRa	0.9999	0.9173	0.1518	0.7881	0.9974	0.9851	0.7758	0.5944
	PSSRb	0.7417	0.9129	0.1274	0.8405	0.7811	0.7951	0.5879	0.5271
	PSSRc	0.7563	0.9371	0.2229	0.9127	0.8068	0.8923	0.5529	0.6924
	PSNDc	0.9207	0.8700	0.4106	0.9046	0.9841	0.8072	0.7608	0.4588
	RR	0.4156	0.9581	0.9406	0.6472	0.8838	0.7913	0.5692	0.0741

). Visual changes in T_c − T_a over the course of the experiment (Figure A1) for an example MR plant, show the appearance of cooler regions by 6 DAI, which extended markedly and became far cooler by 14 DAI. These reductions in T_c − T_a were most pronounced for new foliage that is symptomatic (Figure A1) but smaller reductions in T_c − T_a also occurred for older foliage that did not develop symptoms (Figure A1).

In contrast, little systematic variation was noted between treatments for T_SD across the measurement range. Values of T_SD in the control did exceed those in the MR treatment between 5 and 8 DAI, with these differences becoming significant on 7 DAI. However, by 14 DAI, values of T_SD in the MR treatment were higher than those of the control, but this difference was not significant (Figure 4).

3.3.2. Hyperspectral Indices

Although there was wide variation within each grouping, the xanthophyll indices were generally most sensitive to infection by the pathogen (Table A2). Of all the hyperspectral indices, PRI528 significantly differed the most frequently between treatments. Values of PRI528 were significantly higher within the MR treatment than the control during five measurements dates from 3 to 7 DAI, but converged after this time, with values in the control exceeding those in the MR treatment, but not significantly so, by 14 DAI (Figure 5). Although the other xanthophyll indices did not significantly differ from 3 to 8 DAI, many of these indices exhibited highly significant treatment differences at 14 DAI (Table A2). Among these, treatment differences for DRI PRI and PRIn were most significant (p ≤ 4.2 × 10⁻⁵) among all hyperspectral indices, with values of PRIn in the MR treatment exceeding values in the control at this time (Figure 5).

Most of the indices in other hyperspectral categories were relatively insensitive to treatment from 3 to 8 DAI. However, there were quite a few indices, split across a number of categories, that exhibited highly significant differences between treatments at 4 DAI (Table A2). Among these, the indices that differed most significantly included all three water indices (fWBI, WBI, WI), healthy index (HI), and indices that included reflectance in the blue band (B, BGI, BRI, BF indices and to a lesser extent SRPI, NPCI, SIPI). For a selection of these key indices values in the MR treatment were elevated above control values for SIPI, WI, HI, but reduced below control values for the blue indices, B, BF1 and BRI1 (Figure 5). By 14 DAI, many hyperspectral indices in these categories significantly differed between treatments. Although there were no significant treatment differences for any of the structural indices before this time (Table A2), there were significant differences between treatments for over half (8/15) of these indices at 14 DAI, with the most significant differences occurring for RDVI, TVI, MTVI1, MCARI1 and EVI (p < 0.01). Within other categories, the most significant treatment differences were noted for B, BRI2, CLS and SIPI (p < 0.01), but as previously mentioned, none of these differences were as strong as those noted for xanthophyll indices at 14 DAI.

3.4. Model Predictions

3.4.1. Models Using Thermal Indices

The developed models had low to moderate performance for data obtained during 3 and 4 DAI, with accuracy < 62% and an F1 score of 0.75 over these two measurement dates (

Table 3. Confusion matrix and classification statistics categorised into groupings for models that use thermal indices, narrow-band hyperspectral indices (NBHIs) or a combination of both types (NBHIs + Thermal). Abbreviations for the confusion matrix are: true positive (TP), false positive (FP), false negative (FN) and true negative (TN) and the classification statistics include precision (Prec.), recall, accuracy (Acc.) and the F1 Score. Models with an F1 score of 0.7–0.8, 0.8–0.9 and 0.9–1.0 are classified (and coloured) as having, respectively, acceptable (blue cells), excellent (green cells) and outstanding (gold cells) classification. The highest F1 score for each day after inoculation between the three sets of indices is bolded and where the same value occurs for two sets of indices both values are bolded.

Index	Days after	Confusion Matrix (%)				Classification Statistics
	Inoculation	TN	FP	FN	TP	Prec.	Recall	Acc. (%)	F1 Score
Thermal indices	3	4.8	28.5	9.3	57.3	0.67	0.86	62	0.75
	4	2.1	31.2	8.3	58.4	0.65	0.88	61	0.75
	5	6.4	26.9	6.5	60.1	0.69	0.90	67	0.78
	6	18.8	14.5	13.1	53.6	0.79	0.80	72	0.80
	7	23.2	10.1	5.9	60.8	0.86	0.91	84	0.88
	8	27.2	6.1	3.5	63.2	0.91	0.95	90	0.93
	14	24.8	8.5	4.9	61.7	0.88	0.93	87	0.90
NBHIs	3	17.5	15.9	7.2	59.5	0.79	0.89	77	0.84
	4	22.1	11.2	6.9	59.7	0.84	0.90	82	0.87
	5	17.3	16.0	7.7	58.9	0.79	0.88	76	0.83
	6	22.4	10.9	11.3	55.3	0.84	0.83	78	0.83
	7	17.3	16.0	9.5	57.2	0.78	0.86	75	0.82
	8	12.3	21.1	11.6	55.1	0.72	0.83	67	0.77
	14	21.7	11.6	9.3	57.3	0.83	0.86	79	0.85
NBHIs + Thermal indices	3	13.3	20.0	7.1	59.6	0.75	0.89	73	0.81
	4	22.5	10.8	8.1	58.5	0.84	0.88	81	0.86
	5	16.0	17.3	7.1	59.6	0.77	0.89	76	0.83
	6	22.0	11.3	9.1	57.6	0.84	0.86	80	0.85
	7	25.5	7.9	3.6	63.1	0.89	0.95	89	0.92
	8	27.2	6.1	3.5	63.2	0.91	0.95	90	0.93
	14	30.0	3.3	5.7	60.9	0.95	0.91	91	0.93

). However, by 5 DAI, there was a substantial improvement in model performance (accuracy = 0.67; F1 score = 0.78), as reflected by reductions in T_c − T_a for MR plants (Figure 4). Model performance improved substantially over the next three measurements and reached the highest accuracy and F1 score, of, respectively, 90% and 0.93 at 8 DAI. Model performance was slightly lower at 14 DAI (accuracy = 87%; F1 score = 0.90), but nonetheless the predictions for both 8 and 14 DAI were both classed as outstanding classification (Table 3) and the effectiveness of T_c − T_a at separating treatments during these measurements is clear in Figure 5 and Figure 6. In general, recall was higher than precision for all measurements, and these differences were greatest from 3 to 5 DAI. This difference reflected the higher number of false positive predictions compared with incorrect classification of MR plants (Table 3).

3.4.2. Models using Hyperspectral Indices

In contrast to thermal indices, models created using NBHIs demonstrated excellent classification during the pre-visual period with high F1 scores of, respectively, 0.84 and 0.87 for 3 DAI and 4 DAI (Table 3). The F1 values declined slightly to 0.83 from 5–6 DAI, but these values were still slightly higher than the F1 values for models created from thermal indices for these measurement days. The performance of models based on NBHIs was excellent for 7 (F1 score = 0.82) and 14 DAI (F1 score = 0.85), but was not as high for 8 DAI (0.77). However, over this period the F1 scores for NBHI models were markedly lower than that of models based on thermal indices. The recall exceeded precision for all but one of the measurement days (6 DAI), highlighting the prevalence of false positives among the misclassified plants (Table 3).

3.4.3. Models Using Both Thermal and Hyperspectral Indices

The models using both thermal indices and NBHIs generally included the best features of both datasets which in some cases resulted in synergies in model performance. As with the models with NBHIs only, excellent classification was achieved from 3 to 6 DAI, with F1 scores and accuracies ranging from, respectively, 0.81–0.86 and 73–81% over this period (Table 3). Similar, to the models with only thermal indices, the classification statistics for the models with the combined dataset, had outstanding performance (i.e., F1 score > 0.9) during 7, 8 and 14 DAI, with the F1 score and accuracies ranging, respectively, from 0.92–0.93 and 89–91%.

The DAI when model performance using the combined dataset exceeded that of the separate datasets included 6, 7, 8 and 14 DAI, and these models had F1 scores of, respectively, 0.85, 0.92, 0.93 and 0.93. As with models created using the separate datasets, recall generally exceeded precision for the combined dataset, highlighting that most misclassifications were attributable to false positives rather than incorrectly classified MR plants.

Variables that were selected for models created from the combined dataset were consistent with the best models created from the separate datasets for each DAI. The variable PRI528 appeared in five of the seven models and was among the two most important variables for four of these models (Table 1). As time since inoculation progressed T_c − T_a assumed more importance in the models, improving from 8^th place at 5 DAI, to 3rd place at 6 DAI and was the most important variable for 7, 8 and 14 DAI (Table 3).

Plots of the two most important variables, in each model for selected DAI, highlight the utility of using a combined dataset to improve treatment separation (Figure 6). During the pre-visual stage at 4 DAI there was quite clear separation between most plants using HI and BF1. Similarly, a plot of PRI528 against CUR at 6 DAI, illustrates quite distinct groupings for the two treatments, with respect to these two variables. By 7 and 14 DAI, T_c − T_a was the most important variable. Plots of T_c − T_a against PRI528 at 7 DAI showed excellent treatment separation, while a plot of T_c − T_a against DRI PRI at 14 DAI, showed almost complete separation of the treatments (Figure 6).

3.5. Tree Physiology and Relationships with Thermal Indices and NBHIs

Prior to inoculation, there were no significant differences between the three physiological variables and values for g_s; A and E varied by ≤5% between treatments. In contrast, at 9 DAI all three physiological variables significantly differed (p < 0.05) between treatments. Mean values of these variables for MR plants exceeded those of control plants by 38% for g_s (74.8 vs. 54.4 mmol m⁻² s⁻¹), 33% for E (0.743 vs. 0.557 mmol m⁻² s⁻¹) and 32% for A (4.74 vs. 3.59 µmol m⁻² s⁻¹). The linear regressions of these three variables against the 83 NBHIs and thermal indices detected only five relationships that were significant, that all included one of the two thermal indices. For all three physiological variables, T_c − T_a was the strongest significant predictor (Figure 7) which was negatively correlated with E (p = 0.0006; R² = 0.19), g_s (p = 0.0006; R² = 0.19) and A (p = 0.005; R² = 0.14). There were also significant positive relationships between T_SD and both E (p = 0.03; R² = 0.08) and g_s (p = 0.03; R² = 0.08) (Figure 7).

4. Discussion

This study highlighted the utility of thermal and VNIR hyperspectral close-range imagery for detection of myrtle rust on pōhutukawa, which is an iconic species under threat from the disease in New Zealand. Using the combined dataset, model performance, as described by the F1 score, ranged from excellent between 3–6 DAI to outstanding (i.e., F1 score > 0.9) from 7–14 DAI. Recent research has successfully used thermal imagery and leaf-level hyperspectral data from a hand held spectroradiometer to detect myrtle rust on rose apple, which is one of the most susceptible hosts of the disease [27]. Findings from the study presented here are consistent with this previous research and extend this approach to a more important target species. The collection of canopy-level imagery, rather than leaf-level hyperspectral data, as used in the previous study [27], allows more direct comparison of model performance with the canopy-level thermal imagery. The collection of this canopy-level hyperspectral imagery is also an important step towards scaling up the detection model using a camera that is often used for acquisition of canopy-level hyperspectral data from a UAV.

This study also demonstrates that robust canopy-level hyperspectral data can be collected under artificial lights. In contrast to outside acquisitions, that rely on clear sky conditions, this novel methodology allows rapid crown level changes in NBHIs resulting from pathogen infection to be characterised at a fine temporal scale. Use of such a method is particularly useful when the pre-visual window for detection is relatively narrow as was the case in this study. The set-up used here does not require the geometric or atmospheric corrections typically necessary for UAV-based captures. On the other hand, achieving uniform illumination of heterogeneous targets such as whole plants using artificial lights is quite challenging.

Hyperspectral indices were found to be more effective than thermal indices for pre-visual detection of infection. As the scanned plants comprised a far greater proportion of older resistant leaves than emerging susceptible leaves, the detected changes mostly reflected shifts in this asymptomatic mature leaf cohort. Pre-visual differences in hyperspectral indices were evident from as early as 3 DAI and by 4 DAI there were widespread differences in xanthophyll indices, the blue band, disease indices and water indices. Within these groupings, the pre-visual indices that were most significantly different between treatments, included, respectively, PRI528, BF1, HI and fWBI.

Normalised canopy temperature was the most useful variable for early disease detection and T_c − T_a showed marked reductions in response to infection from 5 DAI onwards. Following energy balance theory, lower leaf temperature results from increased transpiration which is in turn a result of higher stomatal conductance [41] and observed changes in all these three variables conformed to this theory. Austropuccinia psidii is a biotrophic pathogen [42] and as with other rusts it is likely that these physiological changes were the result of pustules rupturing the cuticle [43]. However, it is also worth noting that the physiological and thermal response of plants to infection also varies according to the infection stage. Development by biotrophic pathogens during later stages is characterised by an increase in the proportion of necrotic tissue which results in reductions in transpiration and increases in leaf temperature [43].

There was strong similarity between results reported here and those from the previous experiment that used hyperspectral and thermal indices to detect myrtle rust on rose apple [27]. Hyperspectral indices derived from symptomless older leaves were the most useful indices for pre-visual detection of myrtle rust on rose apple, while thermal indices were the most effective for later detection immediately before symptom expression and in post-symptomatic plants. As data were collected from a spectroradiometer, indices were derived from both the VNIR (400–1000 nm) and shortwave infrared (SWIR; 1000–2500 nm) ranges within this previous experiment. The most important pre-visual hyperspectral indices derived from the VNIR range on the asymptomatic older leaves were almost identical to those identified in the current study and included PRI528, BF1, HI and WBI [27]. Within the SWIR range, four water indices (Ratio975; NDWI1; SRWI, SRWI2) exhibited strong significant treatment differences two and three days before symptoms appeared, which preceded the differences in WBI, WI and fWBI, extracted from the VNIR range, that occurred the day before symptoms [27].

The key indices used here were also well aligned with previous studies that have used hyperspectral and thermal imagery for pre-visual and early detection of stress. Many studies have shown that the Photochemical Reflectance Index (PRI) is a useful pre-visual indicator of stress [44,45,46] and PRI has been previously used for both early disease detection [47,48] and characterisation of disease severity [49]. Indices derived from the blue part of the spectrum often feature prominently in pre-visual detection studies. Chlorophyll has high absorption in this range [50] and consequently many indices that include wavelengths in the blue region (e.g., NPQI, BGI2) are proxies for degradation of chlorophyll [51,52,53] which is a likely impact of infection by A. psidii [27]. Water stress is a common outcome from disease infection [54] and a range of water indices have proven useful for detection of disease [55,56].

In studies where it has been measured, normalised canopy temperature, as expressed through indices such as T_c − T_a or crop water stress index (CSWI), is frequently an important pre-visual and early indicator of disease [25,26,48]. In contrast to our study, increased canopy temperature is often associated with increased disease severity as previous research has often studied pathogens that cause early reductions in stomatal conductance and transpiration [25,26,47,48]. Consistent with our findings normalised canopy temperature was found to detect almond red leaf blotch later than hyperspectral indices [48]. However, research also shows canopy temperature to be an effective early detector of infection by Xylella fastidiosa [25] and Verticillium wilt [26,47] on olive trees and T_c − T_a was the most effective early discriminator of holm oak decline [21]. The relative efficacy of hyperspectral and thermal imagery depends on the physiological response induced by pathogen infection which likely accounts for this variation [48]. A potential issue when working with thermal data is the separation of temperature and emissivity. While there are emissivity differences between healthy and severely stressed vegetation, differences do not usually occur during the early stress stages examined in this study [57]. As a consequence it was assumed that use of a single emissivity value did not negatively influence the results.

Despite the small sample size the presented results are likely to be robust as a thorough cross validation method was used and there was considerable alignment between predictors identified here and in the earlier study on rose apple [27]. As found previously [25,26,47,48], results reinforced the utility of collecting both thermal and hyperspectral imagery. The use of both hyperspectral and thermal indices provided significant synergies in model performance and importantly, bridged the pre-visual with the early detection period. Ideally, hyperspectral imagery should be collected from a camera that captures the full VNIR/SWIR range as many water indices previously shown to be useful in pre-visual detection of myrtle rust occur within the SWIR range [27]. Previous research has also highlighted the utility of SWIR indices for pre-visual and early detection of disease [21].

Detection of myrtle rust on pōhutukawa is likely to be useful within a nursery setting or for biosecurity purposes in areas where there are important amenity trees or trees of significant cultural and ecological value. As nursery settings are relatively controlled environments, these could provide an opportunity to scale up and fine tune the described methods and examine if they are sufficiently accurate and cost effective for disease detection at a larger scale. Within field settings, previous research has developed a deep learning algorithm using regional aerial RGB imagery that can detect pōhutukawa with high accuracy [58]. The resulting boundaries from this detection algorithm could be used to identify pōhutukawa within aerial acquisitions of thermal and hyperspectral imagery.

Acqusitions of thermal and hyperspectral imagery from UAV or fixed wing aircraft will provide the most appropriate spatial resolution for detection of myrtle rust in pōhutukawa within field settings. Almost all field grown pōhutukawa have a canopy radius exceeding 1.5 m [58] which is most aligned to UAV data collection, where the possible resolutions of 15 cm for thermal imagery [59] and 5 cm for hyperspectral imagery [49,59], enables collection of many pixels per tree. Most VNIR cameras that have been developed for UAVs are less than 5 kg [60]. However, SWIR cameras usually are substantially heavier, as are integrated cameras that can obtain imagery across the VNIR/SWIR range [60] and thus captures across the useful SWIR range will likely need heavy lift UAVs. Data acquisition from a fixed wing aircraft flown at a low altitude could also provide adequate resolution hyperspectral and thermal imagery and it is possible to collect these data at resolutions as fine as 30 cm [60,61]. Although there are a number of satellite platforms that collect thermal and hyperspectral imagery these typically acquire data at a spatial resolution of 30 m for hyperspectral imagery [60] and 60–2000 m for thermal imagery [62], which is too coarse for identification of disease on individual tree canopies.

Scaling up is often complex due to a number of variables that influence the response of hyperspectral and thermal imagery, including canopy structure, sun-sensor geometry, atmospheric conditions and changing environmental conditions [63,64]. The CWSI has been widely adopted as a method of normalising thermal imagery to changing environmental conditions and methods for determining this index have been documented [20]. Radiative transfer models that can account for variability in canopy architecture and image acquisition conditions can be used to accurately simulate canopy thermal properties [65] and describe plant traits useful for predicting disease from hyperspectral imagery [49]. Although scaling up has challenges, a number of studies have successfully detected disease during pre-visual or early stages at broad scales using thermal and hyperspectral imagery [21,25,47,48] demonstrating that robust acquisition and processing pipelines can be developed.

5. Conclusions

Pōhutukawa plants that were inoculated with A. psidii showed the first symptoms by 3 DAI, with all plants displaying symptoms by 8 DAI. Indices extracted from regular captures of plant level thermal and hyperspectral imagery were able to pre-visually distinguish between infected and control plants with considerable accuracy. The NBHIs were most effective at pre-visually classifying treatments and models developed between 3 and 6 DAI had excellent discriminatory power. The normalised canopy temperature, T_c − T_a, which was derived from the thermal imagery was more effective than NBHIs at classifying treatments from 7–14 DAI. There were synergies when indices derived from hyperspectral and thermal imagery were combined, and models using a combination of indices, had excellent classification from 3–6 DAI (F1 score range: 0.81–0.85) and outstanding classification from 7–14 DAI (F1 score range: 0.92–0.93). Further research should extend this method to a nursery setting and examine the potential of thermal and hyperspectral imagery for detecting myrtle rust on field grown Myrtaceae.