Stress Detection in New Zealand Kauri Canopies with WorldView-2 Satellite and LiDAR Data

: New Zealand kauri trees are threatened by the kauri dieback disease ( Phytophthora agathidicida (PA)). In this study, we investigate the use of pan-sharpened WorldView-2 (WV2) satellite and Light Detection and Ranging (LiDAR) data for detecting stress symptoms in the canopy of kauri trees. A total of 1089 reference crowns were located in the Waitakere Ranges west of Auckland and assessed by ﬁeldwork and the interpretation of aerial images. Canopy stress symptoms were graded based on ﬁve basic stress levels and further reﬁned for the ﬁrst symptom stages. The crown polygons were manually edited on a LiDAR crown height model. Crowns with a mean diameter smaller than 4 m caused most outliers with the 1.8 m pixel size of the WV2 multispectral bands, especially at the more advanced stress levels of dying and dead trees. The exclusion of crowns with a diameter smaller than 4 m increased the correlation in an object-based random forest regression from 0.85 to 0.89 with only WV2 attributes (root mean squared error (RMSE) of 0.48, mean absolute error (MAE) of 0.34). Additional LiDAR attributes increased the correlation to 0.92 (RMSE of 0.43, MAE of 0.31). A red / near-infrared (NIR) normalised di ﬀ erence vegetation index (NDVI) and a ratio of the red and green bands were the most important indices for an assessment of the full range of stress symptoms. For detection of the ﬁrst stress symptoms, an NDVI on the red-edge and green bands increased the performance. This study is the ﬁrst to analyse the use of spaceborne images for monitoring canopy stress symptoms in native New Zealand kauri forest. The method presented shows promising results for a cost-e ﬃ cient stress monitoring of kauri crowns over large areas. It will be tested in a full processing chain with automatic kauri identiﬁcation and crown segmentation.


Introduction
The deadly kauri dieback disease (Phytophthora agathidicida (PA)) was first officially confirmed as a new pathogen by Beever in 2008 [1] in the Waitakere Ranges west of Auckland and later taxonomically described by Weir et al. [2]. Meanwhile, it was detected over most of the natural distribution range of New Zealand kauri [3]. The endemic kauri trees (Agathis australis) are a key species for the forest ecosystems on the North Island [4], a famous tourist attraction, and are of high cultural significance for Māori [5,6]. Regular monitoring of stress symptoms in tree canopies requires a cost-efficient, objective method that is suitable for covering large areas and able to detect the first signs of stress. Multispectral band and index combinations for the detection of kauri and canopy stress symptoms have been identified by Meiforth et al. [7,8]. Airborne Light Detection and Ranging (LiDAR) data recently became available for the northern kauri forests, so that the main part of the distribution range of kauri trees in New Zealand is now covered. In this study, we analyse the use of WorldView-2 (WV2) satellite data in combination with LiDAR data to detect stress symptoms in kauri trees in an object-based approach on manually segmented kauri crowns.

Kauri Trees and Kauri Dieback Symptoms
Kauri trees show a large phenological variety, from young conical trees with dense foliage to large open dome-shaped crowns (Figure 1). The mature trees develop a massive trunk and large side branches and reach heights of up to 40 m in the study areas. The colour of the foliage in non-symptomatic kauri varies from light green to a darker blue green, depending on the exposition to light and thickness of the wax coating. The flat lanceolate needle-like leaves create a spiky foliage surface. The foliage forms clusters in the medium growth stages and mature crowns show a more open canopy [5], with visible branch material and gaps from an aerial viewpoint. Stress symptoms of PA first become visible as yellowing and thinning of foliage, which exposes small top branches. Later stages of stress lead to structural changes caused by an ongoing loss of foliage, with an increase of visible bare branch material and gaps, until the tree dies [1,9]. The higher gap fraction increases the influence of undergrowth on the canopy reflectance. Small trees with a conical shape first develop a bare top on the main stem, while the foliage in the lower part of the canopy is affected in later stages. However, these types of symptoms can also have other causes, such as continuing drought periods [10], storm damage, or other diseases [1]. Figure 1. Kauri crown size classes: Non-symptomatic small, medium, and large kauri crowns, as a profile on a Light Detection and Ranging (LiDAR) point cloud (a-c) [11] and on 7.5 cm aerial images below (d-f) [12].
Amongst the most widely used vegetation indices for canopy stress detection with broader bandwidths of satellite images are combinations of near-infrared (NIR) and red bands, which are responsive to changes in canopy chlorophyll and leaf biomass. Other widely used indices are based on red and green bands sensitive to stress-related changes in leaf pigment compositions and bands in the NIR spectrum, which are influenced by the amount of leaf and canopy water and thus also structural changes due to foliage loss [13,[24][25][26][27].
A new generation of commercial high-resolution satellites, such as WorldView (WV) and RapidEye (RE), launched in 2009 and 2008, respectively, has provided additional red-edge bands (WV, 705-745 nm; RE, 690-730 nm). In green vegetation, the red-edge region covers a steep increase in reflectance values between the red and NIR1 spectra and shows a high sensitivity to changes in chlorophyll, which causes a blue shift of the central red-edge region. Several studies have successfully applied spectral indices with red-edge bands from satellite data for stress detection in tree canopies. Oumar and Mutanga [28] utilised indices on the NIR and red-edge bands from WV2 data to detect a decline in Eucalypt plantations in South Africa. Ortiz et al. [17] detected early stages of bark beetle attack by combining TerraSAR-X and RapidEye data. Eitel et al. [29] proved that a normalised difference vegetation index (NDVI) with a red-edge band is suitable for detecting early stress symptoms in conifer woodland trees in New Mexico.
Textural information derived from satellite imagery has successfully been used to assess structural features in tree canopies [30][31][32]. It also has improved the analysis of stress symptoms in tree canopies, especially in later stages of decline with structural changes in the canopy [19,23].
Stress symptoms in kauri trees were analysed with airborne hyperspectral data by Meiforth et al. [8]. According to their study, a NIR/red NDVI, followed by indices with bands in the NIR and red-edge range are the most important for describing the full range of stress symptoms. However, pigment-sensitive indices with narrow green and red bands had a higher importance for the detection of stress symptoms in smaller crowns with denser foliage. Indices with NIR bands, such as the water band index (WBI) at 900 and 970 nm, were more relevant to describe stress symptoms in larger crown sizes with an open crown structure. The NIR bands are sensitive to structural changes and changes in the leaf water content. The paper recommends a stratified approach according to different growth stages of kauri.
An appropriate spatial resolution for the target objects is important for the successful use of satellite imagery in order to avoid mixed pixels with shadow, neighbouring vegetation, and soil [33]. In Meiforth et al. [8], a minimum crown diameter of 3 m was defined for stress detection in kauri trees for the 1 m resolution of a hyperspectral image. Fassnacht et al. [34] found that the 5 m pixel resolution of their HyMAP data (Integrated Spectronics) was too coarse to evaluate bark beetle infections at crown level. However, dead trees could still be mapped with a high overall accuracy (84-96%). Immitzer et al. [35] and Pu and Shawn [36] only used tree crowns that could be manually identified on WV2 respective Ikonos images. Immitzer et al. [35] defined at least two pixels (=4 m) per object as the minimum target size for an analysis of tree canopies with a WV2 image. Other studies analysed larger forest stands. Meddens et al. [37] discovered that a spatial resolution of 2.4 m performed best when mapping the infestation of forest stands with mountain pine beetle infection in Colorado. They demonstrated that a higher spatial resolution also increases the variability within the same level of crown damage and can cause a problematic "salt and pepper" effect in a pixel-based analysis. Lottering [38] and Ismail et al. [39] showed that higher-resolution data (1.25 and <1.75 m) are better suited to detecting early stages of stress in eucalyptus stands and a pine plantation, while later stages of stress are best described with lower spatial resolutions of 2.5 and 2.3 m. To achieve a sub-pixel horizontal accuracy in the processing of WV2 data, Aguilar et al. [40] recommended avoiding off-nadir angles higher than 20 • . They achieved the highest spatial accuracies by using ortho-ready products with third order three-dimensional (3D) rational functions with the supplied rational polynomial coefficients, an accurate digital elevation model, and a polynomial adjustment with well distributed ground control points.
An object-based approach is well-suited to handling a high variance in target objects [41][42][43][44]. Several studies have found that object-based methods are superior to pixel-based methods for single tree crown analysis [35,[45][46][47]. They have also been successfully applied for stress analysis and change detections in tree canopies [48][49][50]. For use in forest areas, an object-based approach requires a prior segmentation of crowns with respect to homogenous forest stands. An object-based analysis simplifies the combination of different data sources under the requirement that these are accurately aligned [51][52][53]. However, an essential condition for utilising the advantages of an object-based approach is that the pixel size is significantly smaller than the average size of the objects of interest [43].
The fusion of multispectral bands with a higher-resolution panchromatic channel, so-called pan-sharpening, allows the creation of a higher spatial resolution multispectral image. It can improve the use of textural attributes in an object-based analysis [54]. While some pan-sharpening methods can alter the spectral information, certain techniques, such as the Gram-Schmidt and Ehlers pan-sharpening methods, have been tested and shown to be superior for preserving the spectral values [55,56].
Stress symptoms in tree canopies caused by soil pathogens are expressed as a multifactor gradient with different states of chlorosis and defoliation until only bare branches are left [26]. While classification algorithms only distinguish between discrete classes, regression algorithms give a more detailed description of symptom stages with a continuous output. A regression approach also facilitates a more differentiated evaluation of changes in a time series. Training based on the correlation coefficient allows regression models to handle a high variance in the data. Several studies have successfully used regression approaches to analyse stress symptoms in tree canopies [19,57], predict tree mortality [58], and evaluate foliar moisture [59]. However, a non-linear relationship in the numeric range of reference values causes a shift between the scales of the reference and predicted values. Meiforth et al. [8] enhanced the linearity by rescaling the value range and also suggested testing a two-step approach, in which the dead and dying trees are first identified in a binary classification. The remaining stress levels can then be analysed in a regression approach.
Light Detection and Ranging (LiDAR) systems generate an accurate, high-resolution, three-dimensional point cloud, which enables to segment single tree crowns and homogenous units of forest stands [60][61][62][63][64][65]. LiDAR data allows estimating crown parameters, such as the tree height, crown shape, canopy density, foliage texture, and structural crown characteristics [66,67]. The LiDAR intensity values add a spectral component as the amount of reflected energy at the peak amplitude of the returned LiDAR signal in the near-infrared region, usually at 1064 or 1550 nm. Intensity data have been successfully used to describe tree canopies and detect tree species [61,[68][69][70]. The suitability of LiDAR attributes for assessing structural crown characteristics as a measure of defoliation has been proven in Scots pine forests at both stand level [71][72][73] and single tree level [72]. In combination with biochemical and biophysical information from passive optical sensors, LiDAR data can demonstrably enhance individual tree analysis [35,52,61,74,75]. The combination of different data sources requires accurate spatial alignment. Data from airborne LiDAR systems (ALS) are usually the remote sensing data with the highest spatial accuracy compared, for example, to photogrammetric data [76,77]. However, the accuracy of terrain models derived from ALS ground data is impaired in terrain with steep slopes and vegetation with dense foliage and undergrowth [76,78].

Objectives
Spectral indicators with a sensitivity to stress symptoms in kauri trees have already been defined by Meiforth et al. [8] for use with an airborne multispectral system. This study analyses whether these findings can be implemented in an operational detection strategy based on WV2 satellite data. With LiDAR data available for most of the distribution range of New Zealand kauri, LiDAR-based crown attributes were integrated for crown stratification and to test their use for stress detection. The study pays special attention to differentiating early signs of stress. It has the following objectives: 1.
Test the performance of WV2 attributes for crown-based stress detection in kauri trees and define the recommended minimum crown size; 2.
Test a two-step method on WV2 attributes by first identifying dead and dying trees in a classification and then applying a regression for the remaining stress symptom levels;

3.
Test the performance of LiDAR attributes in combination with WV2 data for canopy stress detection.
A minimum requirement for the performance of the developed method for a value range of symptom levels from 1 to 5 is a root mean square error (RMSE) smaller than 0.5, so that the predicted symptom values stay mainly within one reference symptom level. While the analysis in this study is based on manually edited crowns, the application for monitoring requires additional steps to automatically pre-segment crowns with respect to homogenous stand units for the existing LiDAR data and to define kauri trees. To prepare for this full workflow, we edited the reference crowns for this study on the LiDAR data.

Study Area
The Waitakere Ranges on the west coast of Auckland are amongst the highest PA-affected forest areas, and the first officially confirmed site in New Zealand with a verified PA infection (Figure 2a, [1]). The Ranges are characterised by a warm temperate climate and a rough terrain, with steep slopes and elevations from sea level to 474 m [79]. The three selected study sites ( Figure 2b) cover trees with the full range of symptom stages in both homogenous and mixed stands, with mainly second-growth kauri forests in the Maungaroa area, larger trees in the Kauri Grove area, and all growth stages of kauri in the Cascades.

LiDAR Data and Aerial Images
Airborne LiDAR data were obtained for the study areas on 30 January 2016 by AAM New Zealand, with a Q1560 LiDAR sensor (400 kHz, 58 • field of view) [11]. The sensor achieved five pulses/m 2 on average, which resulted in 35 average returns/m 2 and 0.17 m average point spacing. The ground return, however, was significantly less, with 0.8-1.5 ground returns per m 2 in the dense forest stands. The sensor utilises the 1064 nm wavelength to record intensity. The spatial accuracy was stated as +/−0.1 m standard error for vertical information and +/− <0.5 m standard error for horizontal information at a 68% confidence level. A digital terrain model (DTM), digital surface model (DSM), and ground normalised crown height model (CHM) were generated in LAStools [83] with the "spike-free" method, according to Khosravipoura et al. [84]. Outliers that were classified as "noise" were removed. The method uses all returns to create a triangular irregular network (TIN) from the highest returns downwards. The user defines a "freeze distance" for the maximum length of the triangle sides, and an "insertion buffer" in a vertical downward direction, which defines the points that are included in the TIN creation, before the freezing process starts. All triangles with sides shorter than the freeze distance are preserved and underlying points are ignored in further steps of the TIN creation for the remaining LiDAR points. This method prevents the creation of downward spikes. In the final step, the surface model is created through a linear interpolation of the spike-free TIN. As a rule of thumb, the lengths of the triangle sides should be around three times the average pulse spacing. We calculated two versions of height models with freeze distances of 0.6 m and 1 m and a pixel size of 0.25 m, which is larger than the average pulse spacing of the last returns of 0.23 m for this dataset.
A three-band red-green-blue (RGB) aerial image was acquired during the same flight as the LiDAR data collection, with a 15 cm pixel resolution. It was delivered in two versions: orthorectified on the surface model and the terrain model [11]. Another aerial image was obtained by order of the Auckland Council in May 2017, with three RGB bands at a 7.5 cm resolution and orthorectified on a DTM [12].

WorldView-2 Image
The WV2 image was captured on 15 March 2017 at 11:33 a.m. (local time) during cloud-free conditions with 45.7 • sun elevation and 18.2 • off-nadir angle. It covers 100 km 2 of the Waitakere Ranges ( Figure 2b) and was delivered as standard imagery with nearest neighbour resampling in an "ortho ready" format. The eight multispectral bands have a resolution of 1.8 m, and the pan channel features a 0.45 m resolution (Table 1, Figure 3). They reach up to 1043 nm in wavelength, but do not reach the characteristic spectral pattern of kauri trees in the NIR2 range ( Figure 3).
Atmospheric correction of the multispectral bands was carried out with the Fast Line-of-sight Atmospheric Analysis of Hypercubes (FLAASH) module in ENVI, which incorporates the MODerate resolution atmospheric TRANsmission (MODTRAN) radiation transfer model [85]. The eight multispectral bands were pan-sharpened to a 0.5 m resolution with the Gram-Schmidt method [86] using ENVI's Spectral Processing Exploitation and Analysis Resource (SPEAR) tools. The Gram-Schmidt method performs well in retaining spectral information [55,87] and has been used in several studies with satellite images for the analysis of tree canopies [88,89]. The image was orthorectified using ENVI's orthorectification tool for WV2 images with the rational polynomial coefficients (RPCs) provided by Digital Globe with the image ( Figure A1, Appendix A). The orthorectification included a 1 m DEM (spatial accuracy at 95% confidence: vertical +/−0.2 m, horizontal +/−0.6 m) for the entire area [90]. Remaining deviations were addressed with a shift in the X and Y direction and a georeferencing with a second-order polynomial transformation based on 1102 ground control points (GCPs) for the three areas. The mean control point error for the WV2 image in the three study areas was an RMSE of 0.87 m (X 0.57 m, Y 0.66 m). Shadow areas were defined with a threshold value of 300 on a brightness layer, to match manually identified outer crown shadows while keeping internal crown shadows. The brightness layer was calculated as the average of the blue, green, red, and NIR1 bands [91]. Table 1. WorldView-2 image: wavelengths, bandwidths, and spatial resolution [92]. NIR: near-infrared.

Reference Crown Set
The reference crown set with 1089 crowns of kauri and dead/dying trees is based on the crown locations and attributes acquired during fieldwork in the summer of 2015/16 and complemented and updated during the summer season of 2016/2017 [7,8] (Table 2). The recorded crown attributes in the field include basic metrics, such as the stem diameter at breast height, the cardinal crown spread, and the canopy base height, as well as a detailed description of canopy characteristics, such as an estimated crown density, coverage in percentage classes [93], and the amount of dead branches [93]. In addition, other signs of stress were recorded such as decolouration of the foliage and the number of epiphytes. The crowns used in this analysis have a mean diameter of at least 3 m and were stratified into three growth stages using the mean crown diameter ( Figure 1). An outer buffer of 10% of the mean crown diameter was removed before the analysis to reduce edge effects and effects of a spatial misalignment between the LiDAR CHM and the WV2 image. After removing the shadow areas, only crowns with sunlit areas larger than 50% compared to the full crown size were considered for the analysis. Table 2. Distribution of crown size classes according to their mean diameter per stress symptom level for the refined 7-step reference system with symptom levels from 1 = non-symptomatic to 5 = dead.

Reference Values for Stress Assessment
An overall canopy score for stress symptoms was assessed for each crown at five levels, from 1 for non-symptomatic dense foliage, to 5 for "dead tree". This score was based on the fieldwork and the assessment of aerial images, according to the method described by Meiforth et al. [8]. For a more detailed evaluation of the first stress symptoms, the stress levels 1 to 2 were refined in half-steps ( Figure 4) for this study. Level 1 describes crowns with non-symptomatic dense green foliage, which is mainly found in small and some medium crowns. Level 1.5 applies to crowns with still green, but more open foliage, with exposed gaps and single patches of visible branch material from an aerial viewpoint. This level is typically found in non-symptomatic medium and larger crowns. The stress level 2 features visible smaller branches over the full crown extent and the partial yellowing of foliage. Further increasing gaps and visible small to medium bare branches are the main visible characteristics of level 2.5 in the aerial images, as well as a more intense yellowing in some trees. Level 3 features large branches that become visible as linear structures. Level 4 marks a degree of foliage loss that deteriorates the overall crown architecture and shows a decolouration from yellow-green to brown in the remaining foliage. Level 5 describes dead crowns and also includes different amounts of photosynthetic active undergrowth and epiphytes. Figure 5 gives an impression of the crown polygons for different forest stands on the WV2 image with an indication of the crown stress symptom values.  To account for the fact that the LiDAR data were acquired in 2016, one season before the WV2 image, the aerial images from both years were compared, and crowns that showed visible changes or where the status of the crown could not be identified on both aerial images were removed. Since both summers were relatively moist, the stress symptoms had not changed as rapidly as in a drought situation.

Attribute Calculation
A total of 203 LiDAR attributes was calculated for crown stratification and stress detection. For the stratification of crown sizes, basic crown metrics were derived from the coarser CHM with the 1 m spike-free threshold, such as the maximum crown height and mean diameter. Figure 6 gives an overview of the workflow for the attribute calculation. For stress detection, LiDAR attributes were calculated on both versions of the DSM and CHM (0.6 m and 1 m spike-free threshold). They describe the shape, crown structure, and texture. The range of attributes covers height ratios, percentiles and bincentiles, point density and coverage, shape, roundness, slope, and curvatures, as well as first-and second-order texture measures, such as variance, range, kurtosis, homogeneity, and skewness (see Table A1, Appendix C). These measures were either calculated directly on the LiDAR point cloud for the crown polygons or first calculated as rasters with different kernel sizes. We also calculated statistics from the intensity values of the first returns. Height cut-offs of 6 and 8 m for calculation of the density and coverage on the LiDAR point cloud for low and high stands were derived from the field measurements of the canopy base height. Table A1 in Appendix C gives an overview of the attributes and tools.
Spectral (43) and textural (111) attributes were calculated on the eight pan-sharpened multispectral WV2 bands with the band math tool in ENVI. The original set of 31 calculated vegetation indices is based on literature research. They also include the best performing indices from the hyperspectral analysis presented in Meiforth et al. [8], as far as the wavelengths and bandwidths of the WV2 image allowed for a calculation. In addition to the indices, bands from a minimum noise fracture (MNF) transformation, a three-band Munsell hue, saturation, and value transformation (HSVM), and a brightness layer [91] were calculated. In total, 110 texture rasters were based on the PAN channel and band 7 with occurrence and co-occurrence matrices on different kernel sizes. For more details, see Table A2 in Appendix C.
The rasters were scaled to the 1-99 or 2-98 percentiles, according to the histogram distribution, to remove outliers. The resulting attribute rasters were aggregated on the crown polygons with the zonal statistic tool in QGIS by the mean crown value and the standard deviation for the spatial attributes.

Attribute Selection and Regression Approach
The aim of the attribute selection was to identify the subset of attributes with the highest correlation to the stress symptom levels. First, highly correlated and redundant attributes were removed. The attributes were clustered according to their origin (WV2, LiDAR height, and LiDAR intensity) and function as spectral or structural. In each of the clusters, pre-selection was performed with a Wrapper Subset Evaluator [94,95], which outperformed other subset selection methods, such as the correlation-based feature selection (Cfs) and classifier subset evaluator. We used the open source data mining software WEKA, which was developed by Waikato University in New Zealand [96].
The wrapper method evaluates attribute combinations with a specified machine learning algorithm based on a defined evaluation criterion. We used a random forest (RF) regression in 200 iterations, with the Pearson correlation coefficient (correlation) as the evaluation measure. The RF regression also corresponds to the algorithm used for the analysis in this paper and a former analysis of kauri stress symptoms [8]. RF is a non-parametric algorithm that does not require a certain distribution of data [97]. It can handle a large number of input variables with a high variance and is less sensitive to overfitting, unless the data is noisy [98]. To avoid overfitting, we tested and pruned the RF tree to a maximal depth of eight. We tested several search methods and "best first" with bidirectional selection, and a search termination of eight consecutive non-improving nodes gave the best results. The "best first" method combines a backward elimination with a forward selection by finding the best attribute subset from the Wrapper Subset Evaluator. The pre-selected subsets were complemented with other attributes highly correlated to the stress levels, which were identified with the correlation attribute evaluator in WEKA. The subset selection of all attributes was repeated with different seed values. Depending on the number of attributes selected, the resulting attributes from all subsets were again combined, and the subset selection was repeated with a different range of seed values, in order to reduce the selection to the most relevant attributes. Figure 7 shows four raster images of the selected attributes. The performance of each subset from the final selection round was then tested in an RF regression with a three-fold split and 1000 repetitions, to define the best performing subset. The correlation, the root mean squared error (RMSE), and the mean absolute error (MAE) were calculated as evaluation metrics. When two subsets performed equally well, the one with fewer attributes was chosen. The importance of each attribute in the final selection was calculated as the average impurity decrease for the RF regression.
We also tested an approach that was discussed in Meiforth et al. [8], by first identifying the dead and dying trees in a binary RF classification, before applying a regression for the lower symptom levels.

Results Objective 1: Test the Performance of WV2 Attributes for Crown-Based Stress Detection in Kauri Trees and Define the Recommended Minimum Crown Size
A selection of eight WV2 attributes and the maximum crown height resulted in a correlation of 0.85 and an RMSE of 0.59 for all 1089 crowns with a minimum diameter of 3 m (Table 3, Table A3 in Appendix D). For a minimum crown diameter of 4 m, the correlation increased to 0.89, while the RMSE decreased to 0.48 (Table 3). Figures 8 and 9 show that mainly the small crown sizes and higher stress levels caused outliers (error >1). Table 3. Test of different minimum crown diameters for WV2 attributes for all crowns for seven stress levels. In addition to eight WV2 attributes, the maximum crown height value was included for crown stratification. The performance was tested with a random forest (RF) regression with three random folds in 1000 repetitions and a depth of 8.

Crown Diameter
No.   To analyse the effect of spatial misalignment between the WV2 image and the CHM, we measured the spatial deviation between both datasets for 100 randomly selected crown polygons ( Figure A2, Appendix B). The test resulted in an overall RMSE of 0.87 m, with 60 crowns having no position error. After subtracting a buffer of 10% of the mean crown diameter, 90% of the remaining crown polygons (edited on the CHM) were within 0.5 m of the crown position on the WV2 image. However, the horizontal inaccuracy of six of the 100 crowns on the WV2 image accounted for more than 20% of their mean crown diameter (mean crown size 6.3 m, standard deviation 3.0 m).
The boxplot diagram in Figure 10a confirms the overall good match of the model, with a correlation of 0.89 for all crowns with a diameter larger than 4 m. However, the first quartiles of the half-step levels (1 to 1.5 and 2 to 2.5) and the levels 3 and 4 overlap, while the first quartiles of the other levels are more separated. The diagram also shows a slightly non-linear relationship between the actual and predicted values, especially towards the dead and dying trees. After a rescaling of the value for dead trees from 5 to 7, a more linear relationship could be established, facilitating the interpretation of the results (Figure 10b). The correlation for the rescaled value range remained unchanged, while the RMSE and MAE values increased with the higher value range. The most important attributes are a combination of red/NIR1 bands in a normalised vegetation index (NDVI_75), followed by a ratio of the red and green bands (red-green ratio index (RGRI)) and the standard deviation of the first band of an MNF transformation (Table 4). For mean crown values of a green NDVI (NDVIg) with red-edge and green bands, the red-edge band and a brightness layer were also selected. Further attributes include the mean value of a seven-pixel kernel of the PAN band and the range of a three-pixel kernel on the NIR1 band. Figure 11 illustrates the predicted versus the actual stress values and their match with an NDVI_75 raster.  A comparison with a reference scheme on the basic five stress levels showed a lower performance (correlation of 0.88, RMSE of 0.55) than the seven-level reference scheme with a more detailed assessment of the first symptoms (Tables 5 and A3 in Appendix D). While the NDVIg and red-edge band were not selected for the five stress levels, a ratio of the two NIR bands (modified normalised difference water index (mNDWI)) improved the performance. Table 5. Performance of an RF regression for the basic five-level reference scheme and the refined reference scheme with seven stress levels based on WV2 attributes for crowns with a diameter larger than 4 m. The correlation (corr.), root mean squared error (RMSE), and mean absolute error (MAE) were calculated for an RF regression (depth 8, 200 iterations) in a three-fold random split with 1000 repetitions.

Results Objective 2: Classification to Identify Dead and Dying Trees
Dying and dead trees (level 4 and 5) with a diameter larger than 4 m could be identified in an RF classification based on WV2 attributes, with a user's accuracy of 84.3% and a producer's accuracy of 88.7% ( Table 6). The remaining stress levels from 1 (non-symptomatic) to 3 (medium symptoms) could be described for crowns larger than 4 m, with an RMSE value of 0.37 (MAE of 0.28, correlation of 0.72). The main differences in the attribute selection to the full symptom range are the higher importance of the NDVIg with red-edge and green bands and the selection of the mNDWI on the NIR bands for the first stress levels from 1 to 3 (Table A3 in Appendix D).

Results Objective 3: Test the Performance of LiDAR Attributes for Stress Detection
The combination with LiDAR attributes improved the stress detection compared to only WV2 attributes for all tested reference schemes (Table 7). For crowns larger than 4 m and seven stress levels, the correlation increased from 0.89 to 0.92 and the RMSE was lowered from 0.49 to 0.43. While the removal of very small crowns under 4 m diameter notably improved the correlation and RMSE for WV2 attributes, the correlation for only LiDAR attributes did not change and the RMSE only improved slightly ( Figure 12).
The most important attributes for a combination of LiDAR and WV2 data (Table 8) are the NDVI on the red/NIR1 bands, followed by the ratio between the maximum height and the 50th percentile crown height (R_max_P50) and the average intensity values. The standard deviation of the first MNF band was also selected with a high importance amongst other spatial attributes for crowns with a diameter larger than 4 m.
Both the identification of dead and dying trees ( Table 9, Table A4 in Appendix E), as well as the detection of first stress symptoms (level 1 to 3), improved with additional LiDAR attributes. In the attribute selection for the first stress symptom, the intensity values had a higher importance, while the R_max_P50 and brightness attributes were not selected (Table 8).    Table 9. Accuracies (Acc.) for an RF classification to identify dead and dying trees (level 4 and 5) with WV2 and LiDAR attributes for crowns with a diameter larger than 3 m (1089 total) and larger than 4 m (895 total). The RF classifier was set up with 1000 iterations in a 10-fold cross-validation.

Discussion
The high overall performance of WV2 attributes with a correlation of 0.89 (RMSE of 0.48, MAE of 0.34) for the seven stress symptom levels and crowns with a diameter larger than 4 m diameter showed that WV2 data is well-suited to describe stress symptoms in kauri canopies. The seven-level reference scheme includes a refined description of the first stress symptoms, which is important for an early identification of potentially infected stands.

Discussion: Minimum Crown Size, Spatial Accuracy and Stratified Approach
The high number of outliers in crowns with a diameter smaller than 4 m, especially for dead and dying trees, is an indicator that the multispectral resolution of 1.8 m in the WV2 image causes mixed pixels in small crowns and is too coarse to detect the small bare top branch of severely affected small crowns. For the use of only WV2 attributes on the full range of stress symptoms, a minimum crown size of at least 4 m is recommended, which improves the correlation from 0.85 to 0.89 and the RMSE from 0.59 to 0.48. Stands with kauri of smaller crown sizes should be pre-segmented and analysed in homogenous stand units.
Spatial misalignment between the WV2 image and the CHM also caused errors. This cannot be completely avoided, especially in steep terrain with high crown heights [106]. The subtraction of an outer buffer of 10% of the mean crown diameter helped to reduce the effects of spatial inaccuracy. However, still six of the 100 tested crown polygons showed displacements of more than 20% of the remaining mean crown diameter. To improve the spatial accuracy, the use of a DSM instead of a DEM should be tested for orthorectification. The different appearance of kauri growth stages requires a stratified approach [8]. If only WV2 data are available for regular monitoring, the inclusion of the maximum crown height derived from LiDAR data of former years allows such a stratification.

Discussion: Attribute Performance
According to the findings presented in Meiforth et al. [8], an NDVI with a red and NIR1 band combination (NDVI_75) and a red-green ratio (RGRI) are the most important attributes for stress detection in the full seven-level symptom range (32.5% and 25% importance, respectively), including dead and dying trees. While the NDVI_75 saturates for high chlorophyll and biomass levels, it is able to capture the loss of photosynthetic activity and structural changes caused by foliage loss [99,100]. The also selected RGRI combination, the green NDVI (9.2% importance), and the mean of the red-edge band (3.7%) are more sensitive to pigment changes [103,107]. The standard deviation of the first band of an MNF transformation (12.6% importance) and the range of the NIR1 band (4.3%) provide information on structural and textural crown characteristics, which are influenced by a higher number of large bare branches and internal shadows. The same features are also captured by a selected brightness layer (3.9%). These attributes have a higher importance in a five-level reference scheme with a greater emphasis on advanced stress levels (Table A3, Appendix D).
For the detection of the first stress symptoms in trees with a still intact crown architecture (level 1 to 3), the green NDVI has higher importance. It detects the blue shift of the red-edge point caused by a decline in photosynthetic activity and changes in the foliage pigment concentrations [107]. A modified normalised difference water index (mNDWI) for the NIR1 and NIR2 bands was also selected for the first stress levels, and is sensitive to a reduction of leaf-cell wall scattering and reduced canopy water caused by foliage loss [108].

Discussion: Dealing with Non-Linearity: Rescaling and a Two-Step Approach
While the RF regression is well suited to handle non-linear relationships, the boxplot in Figure 10a shows that the medians of the resulting values for dead and dying crowns are slightly lower than the actual values. This effect of non-linearity is mainly caused by the more extreme reflectance values of dead branches and shadows in the later stages of decline.
A rescaling of the reference values for the dead crowns (i.e., subdividing the range for this crown state in several levels of severity) suggested by Meiforth et al. 2020 [8] leads to a better match between the medians of the results compared to the actual values (Figure 10b), while the overall correlation coefficients remain unchanged. Merging the mapped severity classes again after the prediction process allows to maintain the original value range from 1 (non-symptomatic) to 5 (dead), which also corresponds to the field reference scheme used by the Auckland Council. A change detection however should be based on the original resulting values to utilise the detailed information in the continuous output range of the regression analysis.
The test of a two-step approach showed that the identification of dead and dying crowns in a binary RF classification requires additional LiDAR data for a stress description of crown sizes smaller than 4 m. If only WV2 attributes are available, the minimum crown size should be 4 m mean diameter. Although, even for crowns larger than 4 m, ca. 15% of the dead and dying crowns were still misclassified as less symptomatic. The remaining lower stress levels from 1 to 3 can be described in a regression with the pre-selected attribute sets with an RMSE of 0.37. However, the mistakes from prior classification need to be added to these results, which makes the interpretation more complicated, especially when this method is used in a change detection. For an easier interpretation, we recommend using only a regression approach.

Discussion: Performance of Additional LiDAR Attributes for Stress Detection
The improved performance with additional LiDAR attributes for stress detection can be explained by the additional 3D information in higher spatial resolution, the high spatial accuracy and the spectral information in the intensity values. Moreover, there is no misalignment with the crown polygons, since the reference crowns were manually edited on the LiDAR data.
The ratio between the maximum crown height and 50th percentile height has high importance for reference schemes that include dead and dying crowns, but it was not selected to assess the first stress symptoms. The advanced loss of foliage is characterised by a lower 50th percentile height of the LiDAR returns, while the maximum height in the form of branches still remains high, even in the later stress stages, caused by the remaining (dead) branches. Further selected LiDAR height attributes describe textural characteristics, as an indicator of foliage loss that exposes gaps and bare branches, such as the variance, standard deviation, and cross curvature in CHM and DSM models. The performance of these attributes depends on an accurate spatial alignment, since their values increase on the edge of the crowns (Figure 7c). For the time being, the acquisition of airborne LiDAR data is still too expensive for regular monitoring of stress symptoms. However, the inclusion of the maximum crown height enables a size stratification and improves the performance of the stress detection with the WV2 data. For the stratification, the crown height can also be derived from older LiDAR data.

Discussion: Recommendations for Further Studies
The resulting stress symptom values of the remote sensing analysis need to be interpreted to translate them into health categories by local experts considering the growth stage, the stand situation, the environmental conditions, and other possible causes of stress. Spatiotemporal dynamics of stress symptoms over several years analysed in change detections might also help to distinguish between PA infection and other causes of stress [16,109,110]. We also recommend testing the use of the original, not pan-sharpened multispectral bands, especially for stress detection in medium and larger crown sizes.
Further studies should investigate crown and stand segmentation in kauri forests based on LiDAR data. Zörner et al. [111] developed a method for large tree detection in native New Zealand forest in the Wellington region. Another promising method for canopy segmentation was developed by Wagner et al. [65] in an Atlantic rainforest based on a U-net convolutional network.
A stand-based analysis of stress symptoms with satellite data requires a better understanding of spectral stress responses in associated tree species. Meiforth et al. [7] showed that the effect of waxy, shiny foliage surface in kahikatea, wooden seed capsules in kānuka, and the influence of older needle cohorts in perennial conifers like rimu can lead to confusion with dead and dying crowns [7]. A spectral unmixing of forest stands with known spectra of kauri and other tree species might also help with the interpretation of stress assessment in mixed forest stands [33,112,113].
We also recommend testing the use of other datasets, which offer more cost-efficient options, such as optical height models derived from Stereo-VHR, instead of LiDAR [114]. Additionally, the newly launched constellations of micro-satellites could be tested, which offer opportunities for producing very-high spatial resolution data at a relatively low cost [115].

Conclusions
This study presents a method for detecting stress symptoms in the canopies of New Zealand kauri trees with LiDAR and WV2 data. A unique characteristic of this study is the large representative reference dataset with 1089 manually edited crowns. It covers the full range of stress symptoms and size classes of kauri in different stand situations. A reference scheme for stress symptoms in five levels was refined to seven levels with a higher differentiation of the first stages of stress. A crown diameter of 4 m was determined as the recommended minimum object size for stress detection, to avoid higher numbers of outliers and mixed pixels for smaller crown sizes. The removal of a 10% outer crown buffer helped to reduce edge effects and spatial misalignment. A selection of eight WV2 attributes in combination with the maximum crown height resulted in a correlation of 0.89 (RMSE of 0.48, MAE of 0.34) in an RF regression for crown sizes larger than 4 m. The performance could be further improved by adding LiDAR attributes (correlation of 0.92, RMSE of 0.43). For cost-efficient regular monitoring of stress symptoms, we recommend the use of WV2 satellite data on pre-segmented crowns with a diameter larger than 4 m and homogenous stand units for smaller crown sizes. The maximum crown height should be included in the analysis to allow for stratification. With the newly available LiDAR data for crown segmentation over the main distribution area of natural kauri forests, the results of this study have important implications for cost-efficient large area monitoring of stress symptoms in New Zealand's kauri forests.  Acknowledgments: Our sincere thanks go to all people and institutions who supported this project. We are especially grateful to David Norton, who helped with the administration and provided input for the fieldwork and ecological aspects of the analysis. Nick Waipara, Lee Hill, and Yue Chin Chew from Auckland Council helped to establish the project and provided data. The Kauri Dieback Programme (Planning and Intelligence Team) gave constructive feedback and support during this research. We would also like to thank the staff of the Arataki visitor centre, Fredrik Hjelm from the Living Tree Company, and Joanne Peace for their excellent support and guidance during the fieldwork. Our thanks go to Massey University for the provision of the hyperspectral data acquisition, as well as AAM NZ Ltd. for the acquisition and processing of airborne LiDAR data. Jeanette Allen, Vicki Wilton, and Nicole Gellner helped with the university administration.

Conflicts of Interest:
The authors declare no conflict of interest.
Appendix A. Illustrations for the Orthorectification of the WV2 Image Figure A1. Different steps in the orthorectification illustrated with the PAN channel of the WorldView-2 image, for a forest section in the Cascade area. The red polygons mark crowns of kauri and dead/dying trees that were used in the analysis. The 10% buffer, which was removed for the analysis, is marked with a stippled line. The orange polygons mark crowns from other species and kauri smaller than 3 m diameter that were not used in this study. Figure A2. Spatial offset of crown polygons, edited on the LiDAR CHM, to the WorldView-2 image, measured for 100 randomly selected crowns. The offset is marked in classes with (a) the distance in percent from the crown diameter and (b) the absolute distance in meter. The blue values mark the offset for the full crown polygons, and the green values mark the offset for the inner buffer that was used in the analysis, after a 10% outer buffer was removed.  Appendix D. WV2 Attribute Importance for the Basic and Refined Symptom Levels Table A3. Selected attributes and attribute importance in percent for different stress symptom references (basic 5 levels 1 , first symptoms 2 , refined 7 levels 3 ) based on crowns with a diameter larger than 4 m for WV2 attributes. The attribute importance for an RF regression was calculated as the average impurity decrease and converted to percentages. D6_k5var variance variance of a 5 × 5 kernel on the DSM (spike-free threshold 0.6 m) 13.8 0.17 [104] Cg6_k5var median variance of a 5 × 5 kernel on the CHM (spike-free threshold 0.6 m) 12.3 0.51 [104] RatDMHght ratio between the mean diameter and maximum height 13.6 0.00