Mapping Physiognomic Types of Indigenous Forest using Space-Borne SAR, Optical Imagery and Air-borne LiDAR

: Indigenous forests cover 24% of New Zealand and provide valuable ecosystem services. However, a national map of forest types, that is, physiognomic types, which would beneﬁt conservation management, does not currently exist at an appropriate level of detail. While traditional forest classiﬁcation approaches from remote sensing data are based on spectral information alone, the joint use of space-based optical imagery and structural information from synthetic aperture radar (SAR) and canopy metrics from air-borne Light Detection and Ranging (LiDAR) facilitates more detailed and accurate classiﬁcations of forest structure. We present a support vector machine (SVM) classiﬁcation using data from the European Space Agency (ESA) Sentinel-1 and 2 missions, Advanced Land Orbiting Satellite (ALOS) PALSAR, and airborne LiDAR to produce a regional map of physiognomic types of indigenous forest. A ﬁve-fold cross-validation (repeated 100 times) of ground data showed that the highest classiﬁcation accuracy of 80.5% is achieved for bands 2, 3, 4, 8, 11, and 12 from Sentinel-2, the ratio of bands VH (vertical transmit and horizontal receive) and VV (vertical transmit and vertical receive) from Sentinel-1, and mean canopy height and 97th percentile canopy height from LiDAR. The classiﬁcation based on optical bands alone was 72.7% accurate and the addition of structural metrics from SAR and LiDAR increased accuracy by 7.4%. The classiﬁcation accuracy is su ﬃ cient for many management applications for indigenous forest, including biodiversity management, carbon inventory, pest control, ungulate management, and disease management.


Introduction
Indigenous forests cover 24% of New Zealand, primarily in mountainous and hilly terrain, representing a 70% reduction from the pre-human state (ca. 800 years ago) due to fire, forest clearance, and logging [1]. These forests provide highly valued ecosystem services, including recreation, sense of belonging, tourism, and important provisioning and regulating services, such as wild foods, fresh water, and regulation of both climate and soil erosion [2]. They contain much of the unique biodiversity in New Zealand where 78% of native plants are endemic and 91% of terrestrial Animalia are endemic [3]. The evergreen forests are warm and temperate in the far north and progress through to cool and temperate in the south [4]. Unfortunately, many plants and animals introduced to New Zealand have become pests in indigenous forests. Some browsing animal pests (such as Australian brushtail possums, pigs, and deer) impact forest plants directly, whereas animal predators (such as cats, rats and stoats) impact the indigenous fauna in the forests. Recently, two tree diseases-kauri dieback 2011 with 24-cm wavelength and acquiring polarimetric data (horizontal and vertical transmit and receive) [20].
Another remote sensing technique able to measure canopy structure is airborne Light Detection And Ranging (LiDAR), also referred to as Airborne Laser Scanning (ALS). In ALS, laser pulses are emitted downwards to the ground and reflected by canopy features back to the sensor, thereby capturing vertical information on canopy structure [21]. In 2013 and 2014, the Wellington Regional Government GIS Group conducted an extensive LiDAR survey with a minimum point density of 1.3 points/m 2 and a vertical accuracy of ±0.15 m over the Wellington region (812,000 ha). While their main intention was to build new, large-scale topographic maps of the region [22], the survey also provides an unprecedented opportunity to study forest composition and individual tree parameters on a large scale.
In addition to imagery quality and fitness for purpose [23], the choice of classification method is also important. Nonparametric methods are proving to be effective for classifying forest types [24], especially Random Forests and Support Vector Machine (SVM) classification, due to their ability to synthesize classification functions on discrete and continuous data sets, and to deal with unbalanced data, and also to their insensitivity to noise or overtraining [25,26]. These properties enable the combination of spectral data with other types of data to improve the accuracy of classifying forest types. However, research into combining satellite optical data with other remotely sensed data relating to vegetation structure, such as LiDAR and SAR, is incipient [27,28]. Classification of objects or segments rather than pixels, termed GEOBIA (Geographic-Object-Based Image Analysis), has improved classification accuracy in forest applications [29][30][31][32][33]. Temporal analysis has also proved useful [12,[34][35][36][37].
In this paper, we investigate how the classification of forest physiognomic types using Sentinel-2 may be improved with the addition of other imagery relating to vegetation structure. We consider radar backscatter from Sentinel-1 and PALSAR (Phased Array L-band Synthetic Aperture Radar), and canopy height metrics (mean and 97th percentile) from airborne LiDAR. These data are available for the whole of the Wellington region, New Zealand-the study area. Fusion of the data is made possible by calculating mean values within segments of typical size 5 ha, that is, using a GEOBIA approach. We use SVM classification because of its ability to deal with unbalanced data and superior insensitivity to noise and overtraining [25,26]. With SVM, we explore the best combinations of Sentinel-1, Sentinel-2, PALSAR, and CHM (canopy height model) data by five-fold cross-validation. Using the best performing classifier, we produce a map of forest physiognomic types in the Wellington region along with summary statistics.

Study Area and Forest Physiognomic Types
The Wellington region, in the lower part of the North Island (Figure 1), New Zealand, has 23% of its 812,000 ha of land, covered in indigenous forest. In pre-Māori times, nearly all the Wellington region was covered by indigenous forest, but after their arrival (ca. 1000 A.D.), a significant proportion was burned for shifting agriculture, and after the arrival of Europeans (ca. 1840 A.D.) an even larger proportion was burned and cleared for pastoral agriculture. Most of the indigenous forest is now confined to the protected mountainous areas of the Tararua, Rimutaka, and Aorangi ranges. Many remnants of indigenous forest are spread throughout the rest of the region, and there are large blocks of shrubland reverting to indigenous forest on the east coast. The climate of the Wellington region is mild and wet: annual rainfall varies between 800 mm in the eastern Wairarapa and 7000 mm at the tops of the Tararua range (ca. 1500 m), and mean monthly temperatures at Wellington city range between 8 • C in July and 16 • C in January.
The forests of the region are dominated by various mixtures of species from three groups: conifers, all from the Podocarpaceae family; broad-leaved evergreen species from a wide range of families; and Southern beechs (Nothofagaceae). Forests dominated by broad-leaved trees with scattered emergent podocarps (Broadleaved-podocarp forest) occur on most lowland sites (podocarp abundance has been reduced by historical selective logging in many areas). Broadleaved-podocarp forest is usually luxuriant with small trees, shrubs, and ferns forming the understoreys with abundant mosses, lichens, lianes, and epiphytes. There is a shift to beech dominance at higher elevations, or on sites with dry climates or infertile soils. Beech forest tends to be dominated by just one or more of the southern-beech species and usually has lower diversity of associated plant species. In this paper, we focus on distinguishing forests based on the relative proportions of tree species of differing physiognomy in the canopy. We focus on the five combinations that occur in the Wellington region: Broadleaved-podocarp forest; Beech-broadleaved-podocarp forest; Beech-broadleaved forest; Broadleaved forest; and Beech forest. Podocarp forest is uncommon in the Wellington region, but common in other regions.

Processing of Sentinel-2 Imagery to Standardised Reflectance
The radiance (i.e., apparent brightness) of vegetation as seen by a satellite is the product of irradiance and the vegetation reflectance, attenuated through the atmosphere. Both the irradiance and reflectance are influenced by the position of the sun relative to the slope normal of the vegetation surface. As reflectance is also influenced by the viewing direction of the satellite relative to the slope normal, in a satellite image of hilly or mountainous terrain, the apparent radiance of vegetation is strongly influenced by topography. For use in vegetation classification it is desirable for a satellite image to represent only the properties of the vegetation, or other reflecting surfaces. To achieve this, it is necessary to process the satellite image to standardised spectral reflectance, which is the spectral reflectance the vegetation would have on flat terrain with the sun and satellite in standard positions (nadir view for the satellite, and a solar elevation of 50 degrees). We processed the satellite imagery to standardised reflectance using the method of Dymond and Shepherd [38], which requires physical modelling of radiation from the sun and sky through the atmosphere, reflection of the light from the vegetation canopy, and the transmission of the reflected light through the atmosphere to the satellite sensor.

Segmentation of Sentinel-2 imagery
A cloud-free mosaic of the Wellington region (99.5% cloud free) was produced by extraction of highest priority cloud-free pixels from the Sentinel-2 images (Section 3.1), with highest priority given to dates nearest the middle of November 2017. The mosaic was then segmented using an iterative elimination process. An initial segmentation was produced from k-means clustering (k = 60), with small clumps being eliminated iteratively, smallest first, until a minimum mapping unit was achieved. The simple method produces results comparable with other well-known algorithms [39] but has the advantages of scalability to large datasets and a target minimum mapping unit. The method can retain smaller segments than the minimum mapping unit where they are significantly spectrally distinct from neighbours, such as small water bodies. The algorithm has been made available in the open source software library RSGISLib [40]. The segments are used as a framework for estimating robust means of Sentinel-1, Sentinel-2, and PALSAR imagery, and of canopy height from LiDAR data.

Processing of Sentinel-1 SAR imagery
We use the Sentinel-1 Level-1 Ground Range Detected (GRD) product consisting of focused SAR data that has been detected, multi-looked and projected to ground range. Individual orbit slices are downloaded for the study region and stored on the NeSI high-performance computer (HPC) facility [41]. A pre-processing workflow has been developed to prepare the uncalibrated SAR data for use in this study. The following steps are performed using routines built into the ESA Sentinel Application Platform (SNAP) version 6 and run on the HPC in parallel for individual orbits: (1) thermal noise removal; (2) satellite orbit correction; (3) radiometric calibration; (4) mosaic of individual orbit slices; (5) speckle filtering of radar noise using a single-product Lee-filter (window size: 7 × 7, sigma: 0.9, target window size: 3 × 3); (6) radiometric terrain flattening using NASA's SRTM elevation model; (7) geometric terrain correction and re-projection to NZTM (New Zealand Transverse Mercator) with 10 × 10 m pixels; and (8) conversion to decibel (dB) units and clipping off noise along the near and far range edges (using a custom polygon boundary for each orbit). For radiometric and geometric terrain correction NASA's SRTM 1 Sec product [42] is used.

Processing of PALSAR Imagery
The PALSAR Global Mosaic [24,43] is a high-resolution slope-corrected and ortho-rectified L-band gamma-naught map. The data used for this study were acquired at 10 m resolution in a fine-beam dual-polarization mode (HH and HV) (H=horizontal and V=vertical) for the year 2007. The data were resampled to the NZTM map projection, mosaicked, and calibrated to backscatter power using published calibration factors [32]: where DN is the 16-bit unsigned value from the ALOS-PALSAR product. Backscatter was then converted to decibels.

Production of Canopy Height Model
A LiDAR survey was conducted over the entire Wellington region [22], primarily in early 2013 but with some additional aircraft flights later in 2013 and 2014, depending on weather. The LiDAR scanner was an Optech Airborne Laser Terrain Mapper ALTM 3100EA flown at a nominal height of 1000 m above the ground. Target point density was 1.7 points per square meter (ppsm) with 50% swath overlap to ensure the minimum specification of 1.3 ppsm and a vertical accuracy of ± 0.15 m. While these were the minimum requirements, the actual point density of the dataset ended up higher, ranging between 4 and 6 ppsm on average and > 6 ppsm in regions with overlapping flight paths. Coincident with the LiDAR survey was an aerial photographic survey with blue, green, red and near infrared bands, which was ortho-rectified to 30 cm pixels.
The 1261 LiDAR flight lines of point cloud data were merged, tiled and further processed using the open-source LiDAR processing library, Sorted Pulse Data Library (SPDLib) [44] as described in [45]. A canopy height model (CHM) at 1 × 1 m pixel resolution was generated with outliers being removed by applying a 5 × 5 m Gaussian median filter and screening out pixels below 0.5 m (undesired low vegetation) and above 60 m (erroneous high trees).

Ground Data
Data from 580 vegetation plots in the Wellington region having tree species observations were extracted from the National Vegetation Survey (NVS) database [46]. Where plots are permanent, only the most recent measurement was taken. These spanned the years 1962 to 2017, with 87% of the plots being measured since 1980. Plots subject to experimental treatments (e.g., exclosures) were excluded.
Plots are of size 20 × 20 m and contain measurements of diameters at 1.35 m for all trees of diameter over 20 mm. On 276 of the plots, relevé data (species abundance in cover classes within height tiers) are available [47]. These data were summarised by listing (i) the five most abundant tree species by basal area and cover; (ii) relative basal area and cover of trees classed as podocarps, beech and broadleaved species; and (iii) the vegetation alliance to which the plot belonged. Summaries informed assignment of the plot to a physiognomic type. A polygon-shapefile of physiognomic forest type was generated from this summary data by overlaying visually, using ArcGIS, the plot location on the orthophotography (i.e., the natural colour and colour infrared images acquired together with the LiDAR measurements) and on a colour rendition, according to the height of the canopy height model. The plots typically are located in sets, usually along altitudinal transects located within catchments using a two-stage random/systematic sampling scheme. It was possible to determine the boundaries of the physiognomic type associated with these sets of plots from the visually apparent relationships with the canopy height model and orthophotography. Where the physiognomic type derived from the plot data was inconsistent with the canopy height model and orthophotography (which could arise owing to successional change since the plot was measured or inaccuracies in location), the physiognomic type was inferred from the imagery. We drew polygons of homogeneous physiognomic types of typical size 5-10 ha. These types (hereafter termed 'ground data' were assigned to the segments (i.e., the raster attribute table of the segment image) derived in Section 3.3 if a segment was wholly or primarily contained (>90%) in a polygon.

Support Vector Machine Classification
The mean standardised reflectance of each Sentinel-2 band was calculated for each segment and stored in the raster attribute table of the segment image. Likewise, the mean backscatter in each Sentinel-1 band and in each PALSAR band was calculated for each segment and stored in the raster attribute table, along with the mean and the 97th percentile (below which 97% of heights are below) of the canopy height model. Together with the ground data (449 segments in total), these data formed the training and test data. We explored the best combinations of Sentinel-1, Sentinel-2, PALSAR, and CHM data by performing five-fold cross-validation, where the ground data is randomly partitioned into 5 equal sized subsamples and each subsample used to test the remaining 4 subsamples assigned to training. The five-fold cross-validation was repeated 100 times with different random selections of subsamples to get an average overall classification accuracy (which we seek to maximise). For the SVM classification, the R-package e1071 is used as interface to the SVM implementation in the LIBSVM software [48]. A linear kernel and default parameters for cost (= 1) and gamma (= 0.1) were used. Table 1 shows all the optical bands, radar bands, and CHM metrics that make a significant contribution to cross-validation accuracy when included in a SVM classification with a linear kernel (significant meaning greater than 1% contribution when combined with other like and significant data, that is, with other optical, with other canopy height, or with other radar data). When all significant optical and radar bands and CHM metrics are included, the accuracy of the SVM classification is 80.05%. The relative importance of each band and metric is shown in the difference between using all bands and metrics and with leaving one out. The most important band/metric is CHM97, the 97th percentile of canopy heights in a segment, with a 4.27% contribution to accuracy. The most important optical bands are B02, the blue band, with a 1.85% contribution to accuracy, and B12, the 2nd short-wave infrared band, with a 1.85% contribution to accuracy also. The most important radar band is the Sentinel-1 VH/VV with a 1.63% contribution to accuracy.

Results
When all the bands and metrics are combined, the red edge band, B05, and PALSAR HH become less significant and make a small contribution to the final accuracy, with B05 making a negative contribution. Table 2 repeats the analysis of Table 1 but omits B05 and PALSAR HH. When only these bands and metrics are combined, the accuracy of the SVM classification increases to 80.48% from 80.05%. The most important band/metric is still CHM97, with an increased 5.19% contribution to accuracy. The most important optical band is now solely B02, with an increased 2.15% contribution to accuracy. The most important radar band is Sentinel-1 VH/VV with an increased 1.80% contribution to accuracy. Table 3 shows the cross-validation accuracy when using all the significant optical bands, then adding the significant CHM metrics to the significant optical bands, and then adding the significant radar bands to the optical and CHM metrics. The optical bands alone achieve 72.67% accuracy. With the addition of the CHM metrics this is increased to 78.48%, which is a significant increase of 5.81%. With the addition of the significant radar bands to the optical and CHM metrics the accuracy is increased to 80.05%, another significant increase of 1.57%. Table 1. Five-fold cross-validation accuracy (100 repeats) of predicted forest physiognomic types from support vector machine classification by leaving out one from all the significant bands (i) Sentinel-2 spectral bands (2,3,4,5,8,11,12), (ii) mean and 97th percentile of canopy height model, and (iii) Sentinel-1 VH/VV and PALSAR HH bands. Grey shade shows which bands were used. The SVM classifier established above was trained on all locations for which ground data were available and then applied to all the image segments of indigenous forest in the Wellington region (according to the Land Cover Database [11]) to produce a regional map of forest physiognomic types. Figure 2 shows a subset of this map in a transect across the Tararua ranges. The spatial distribution of physiognomic types is consistent with common observations. Broadleaved-podocarp forest is often in valley floors and neighbouring slopes. Progression upwards in elevation corresponds with a transition to Beech-broadleaved-podocarp forest and then to Beech-broadleaved forest. Near the tree line, the forest is usually pure Beech. On the western flanks of the Tararua ranges, near the coast, Broadleaved forest predominates. Table 3. Five-fold cross-validation accuracy of predicted forest physiognomic types from support vector machine classification from three distinct groups, the significant optical bands, the significant LiDAR metrics, and the significant radar bands. Grey shade shows which bands were used.  Table 3. Five-fold cross-validation accuracy of predicted forest physiognomic types from support vector machine classification from three distinct groups, the significant optical bands, the significant LiDAR metrics, and the significant radar bands. Grey shade shows which bands were used.    Table 4 shows the areas of each forest physiognomic type as mapped in the Wellington region. The areas of the physiognomic types by height class (according to the 97th percentile of the CHM), short (< 20 m), medium (20-30 m), and tall (>30 m). Beech-broadleaved forests are the most common, with 63 thousand ha, and Broadleaved-podocarp forests are the second most common, with 35 thousand ha. Beech forests and Beech-broadleaved-podocarp forests are less common, both with areas of 16 thousand ha. Broadleaved forests are uncommon, with only 8 thousand ha. Podocarp forests do not exist in the Wellington region, (but do in other regions of New Zealand). Table 5 shows the common forest alliances within the physiognomic types in the Wellington region according to the forest plot data (note that no forest plots occur on the relatively uncommon broadleaved forest).  Height classes differ in their distribution among the physiognomic types. In Broadleaved-podocarp forest, medium is the most common height class, but short and tall are also common. In Beech-broadleaved-podocarp, medium and tall are common height classes, but there is no short. In Beech-broadleaved forest, medium is the most common height class, but short is also common, and there is very little tall. In Broadleaved forest, short and medium are common height classes, and there is very little tall. Likewise, in Beech forest, short and medium are common height classes, and there is very little tall.

Discussion
The 80.5% accuracy (rounded up from 80.48%) of the Wellington map of forest physiognomic types is a significant improvement over the 50.0% accuracy of a current national map of forest physiognomic types, EcoSat Forests [10] (published at a scale of 1:750,000), whose accuracy in the Wellington region we assessed using the NVS forest plots. EcoSat Forests were derived from a combination of Landsat imagery and the vegetation component of the New Zealand Land Resource Inventory (NZLRI). The Landsat imagery was used to bring the NZLRI vegetation up to date (i.e., to 2002/2003 from the 1970s) and to define forest boundaries within land units. The forest codes in the NZLRI do not have the mixed classes, Beech-broadleaved-podocarp forest and Beech-broadleaved forest, which are common in the Wellington region, except if they occupy separate parts of a land unit. It is not surprising, therefore, that there is a significant improvement in the method presented in this paper, which maps mixed forest physiognomic types directly. This was one of the reasons that EcoSat Forests was published at a scale of 1:750,000. However, it has been used at a range of scales, including 1:50,000, because it is the only national digital map of forest physiognomic types with reasonable currency (~2000).
While the accuracy of mapping forest physiognomic types is much improved with the method presented in this paper, it is still not high enough for some forest management applications, such as disease management. In Table 6, we list major management applications for indigenous forest mapping and assess (by canvassing expert opinion) whether the 80.5% classification accuracy would make the map of forest physiognomic types fit for purpose. Fitness of purpose is assigned a 'Yes' if a significant contribution is made to the application and 'No' if a minor contribution is made. Fitness for purpose is assigned 'Yes/no' if a significant contribution is made to some aspects of the application while a minor contribution is made to other aspects. The method presented in this paper (labelled 'Forest Physiognomic Types') is fit for purpose for four of the management applications, which compares favourably with one for EcoSat forests and none for LCDB. However, the method presented here has been applied to only one region, Wellington, and so EcoSat forests and LCDB will continue to be used for national applications until a national map of forest physiognomic types is available.
The Wellington region was chosen as the study area because it has a large proportion of indigenous forest and is one of the two regions in New Zealand with complete LiDAR coverage. Several other regional councils are currently acquiring regional LiDAR. The New Zealand government has recently recognised the utility of LiDAR data by announcing a nineteen-million-dollar funding boost for regional acquisition of LiDAR. All New Zealand may have LiDAR coverage within a few years, making a national map of forest physiognomic types possible using the method presented in this paper. There may be regional differences in the SVM classifier due to a gradual change in climate from north to south and also different soils and topography, but there is a good spread of vegetation plots throughout New Zealand that may be used to refine regional calibrations of the classifier.
The addition of structural information to optical data improves the classification of forest types. Indeed, even with just the ratio of Sentinel-1 VH over VV band and the 97% percentile of canopy height band, there is nearly complete discrimination between Broadleaved-podocarp and Beech forests (Figure 3). This is due to the very different canopy structure of these two physiognomic types. Broadleaved-podocarp forest commonly has tall podocarp trees in a sparse upper canopy (typically over 25 m tall), with a lower, 15-m tall secondary canopy. The canopy surface appears rough and heterogeneous and varies markedly in height. In contrast, Beech forest commonly has a smooth upper canopy with most trees reaching to the same upper canopy level at 20 m. Therefore, we would expect the ratio of volume scattering to surface scattering to be higher for the Beech forest than for the Broadleaved-podocarp forest, as seen in Figure 3.    Support Vector Machine classification has handled the addition of structural metrics to optical data satisfactorily and robustly, with the accuracy of training and test data being the same. It compares favourably to other classifiers which achieved lower classification accuracies in the R Caret package, such as Random Forests (71.5%), Extreme Random Trees (78.0%), Neural Networks (73.8%), Classification Trees (C5.0) (73.8%), Classification Trees (CART) (62.5%), Boosted Regression Trees (73.8%), Linear Discriminant Analysis (77.3%), K-Nearest Neighbour (70.5%), and Parallelepiped Classifier (53.3%). While Extreme Random Trees handled the addition of structural metrics well, as did the Support Vector Machine classifier, it was not robust -the accuracy of training data (100%) and test data (78.0%) was very different. Other classifiers exhibiting a lack of robustness were Random Forests and Classification Trees (both C5.0 and CART). The robustness of the Support Vector Machine is desirable when it is applied in regions where training data are sparse.
It would be possible to further increase the mapping accuracy of our method by adding environmental ancillary information, such as elevation from a DEM (Digital Elevation Model), or a soil map [49]. However, we have chosen to limit the method to using direct measures of vegetation, rather than indirect environmental drivers. This is because accuracy is not everything-the ability to apply the method in other regions with minimum effort is also important, as is the ability to interpret results. When ancillary information is used, the complexity of the classifier increases, and overfitting becomes more likely when ground data are limited. Extension of the classifier to other areas without extensive ground data then becomes risky. Artefacts can also appear when forest boundaries are created by environmental thresholds, such as an elevation contour. While we prefer not to use environmental ancillary information, because we are aiming to produce a nationally applicable method, for applications in specific localities, it might be useful, even necessary. We will however, investigate the use of temporal spectral signatures, once the length of record of Sentinel-2 is sufficiently long-to further increase mapping accuracy and associated contribution to management applications.

Conclusions
Radar and canopy height metrics, which are related to vegetation structure, may be combined with optical imagery to improve the accuracy of forest classification. We combined C-band Sentinel-1 (VH/VV) radar, LiDAR derived mean canopy height, and 97th percentile of canopy height, with Sentinel-2 optical imagery to improve the mapping accuracy of an SVM classifier for mapping forest physiognomic types from 73.1% to 80.5%. The optical bands found useful for mapping these physiognomic types were blue (490 ± 50 nm), green (560 ± 25 nm), red (665 ± 20 nm), near infrared (835 ±65), short-wave infrared (1610 ± 70 nm), and short-wave infrared (2185 ± 120 nm). The 80.5% classification accuracy of forest physiognomic types achieved by the combination is sufficiently high to contribute to several management applications for indigenous forest including biodiversity management, carbon inventory, pest control, ungulate management, and disease management.