Invasive Shrub Mapping in an Urban Environment from Hyperspectral and LiDAR-Derived Attributes

Proactive management of invasive species in urban areas is critical to restricting their overall distribution. The objective of this work is to determine whether advanced remote sensing technologies can help to detect invasions effectively and efficiently in complex urban ecosystems such as parks. In Surrey, BC, Canada, Himalayan blackberry (Rubus armeniacus) and English ivy (Hedera helix) are two invasive shrub species that can negatively affect native ecosystems in cities and managed urban parks. Random forest (RF) models were created to detect these two species using a combination of hyperspectral imagery, and light detection and ranging (LiDAR) data. LiDAR-derived predictor variables included irradiance models, canopy structural characteristics, and orographic variables. RF detection accuracy ranged from 77.8 to 87.8% for Himalayan blackberry and 81.9 to 82.1% for English ivy, with open areas classified more accurately than areas under canopy cover. English ivy was predicted to occur across a greater area than Himalayan blackberry both within parks and across the entire city. Both Himalayan blackberry and English ivy were mostly located in clusters according to a Local Moran’s I analysis. The occurrence of both species decreased as the distance from roads increased. This study shows the feasibility of producing highly accurate detection maps of plant invasions in urban environments using a fusion of remotely sensed data, as well as the ability to use these products to guide management decisions.


INTRODUCTION
Many human-dominated landscapes are invaded by non-native species, causing adverse impacts for native flora and fauna (Wilcove et al., 1998;Ehrenfeld, 2003), public health (Mack et al., 2000;Laaidi et al., 2003;Jones et al., 2009), and on ecosystem services including agricultural production, water filtration, recreation and tourism, flood mitigation, and cultural services (Pimentel et al., 2000). The decreased resilience and compromised integrity of these ecosystem services caused by invasive plant species costs at least US$34 billion annually in the United States alone (Pimentel et al., 2005;Pejchar and Mooney, 2009). Reducing the negative effects of plant invasions by restricting their spread is difficult as invasive plant species can distribute quickly through various mechanisms such as greater dispersal abilities (Kowarik, 1995;Pysek and Hulme, 2005), heightened fitness from hybridization (Ellstrand and Schierenbeck, 2000), or quicker germination (Gundale et al., 2008) than native species.
Due to the detrimental impacts of invasive plant species, city authorities are required to initiate programs to control or eradicate the plants. Notably, managers in urban areas are particularly interested in controlling invasive plants (Pysek, 1998;Pimentel et al., 2000;Pysek and Hulme, 2005), as novel habitats and increased habitat disturbance associated with urban environments provide areas where certain invasive species can thrive (Gundale et al., 2008;Lampinen et al., 2015;Hawthorne et al., 2015). In addition to being a hotspot for many plant invasions, urban ecological areas are often the most affected by them. As more people migrate to urban areas and existing urban areas are densified, ecosystem services that urban forests and other urban natural areas provide, such as increasing recreational opportunities (Arnberger, 2006), reducing the urban heat island effect and various types of pollution (Akbari et al., 2001;Escobedo et al., 2011), improving air quality (Akbari et al., 2001), reducing storm water runoff (Mcpherson et al., 2005), improving esthetics (Price, 2003), and increasing opportunities for social interaction (Kuo, 2003) will become increasingly important (Hough, 2014). It is therefore critical that the integrity of these urban natural areas is maintained, partially by controlling the spread of invasive plant species.
Detecting invasive species is typically undertaken using field surveys when field crews are available, which may be costly and have inconsistent methodologies over space and time. Additionally, field crews cannot survey large areas in sufficient detail to develop accurate maps of invasion locations. As such, different approaches are needed to monitor and map the distributions of invasive species. Remote sensing technology, specifically, light detection and ranging (LiDAR) data and hyperspectral imagery, can produce maps of the distribution of invasive plant species, augmenting previous detection approaches (Huang and Asner, 2009;He et al., 2011;Rocchini et al., 2015). LiDAR sensors are active remote sensing systems, meaning that they require and use their own energy rather than passively detecting solar energy. These sensors send pulses of electromagnetic radiation toward objects and record the time it takes for the pulse to reflect off objects and back to the sensor. Because the pulse travels at the speed of light, the distance to the object can be calculated from the return time (Bachman, 1979). In contrast, hyperspectral imagery is a passive remote sensing technology that uses the objects' reflections of solar energy as its inputs. This type of imagery records reflectance values of an object over numerous narrow spectral channels, providing a detailed spectral signature of an object (Shaw and Burke, 2003).
Light detection and ranging data alone has been used for species habitat modeling and detection. Singh et al. (2015) mapped the distribution of an understory invasive species in urban forests using various LiDAR-derived orographic predictors as well as LiDAR-derived forest structural characteristics, such as return height variance, standard deviation, and mean to gather information about the various height strata in a stand. Some of these vegetation characteristics were found to be predictors of invasive species presence. In addition to its applications toward mapping plant distributions, LiDAR data has also been shown to be able to identify tree species and quantify the density of shrub cover using the intensity values of the LiDAR returns (Kim et al., 2009;Ørka et al., 2009;Wing et al., 2012). However, these analyses require extensive calibration (Kashani et al., 2015) and have not identified individual shrub species, thus additional information, such as detailed spectral information, may be warranted when available and when mapping distributions of shrub species.
Hyperspectral imagery excels where LiDAR data is often hindered, by offering a number of spectral channels in the visible (400-700 nm), near-infrared (700-1350 nm), and shortwave infrared (1350-2500 nm) wavelength, and narrow spectral channel widths (Adam et al., 2010). Despite the wealth of information contained in hyperspectral datasets, the applicability of using such imagery is highly dependent on the spatial resolution of the imagery and the target of analysis, as at coarse resolutions spectral signals mix and may impede finescale species detection (Roth et al., 2015). Peña et al. (2013) found that increasing the pixel size of hyperspectral imagery from 1.2 to 2.4 m markedly decreased the spectral separability between tree species and decreased a species classification's accuracy by up to 25%. Early studies investigating hyperspectral detection of invasive plant species used instruments covering the visible through shortwave infrared wavelengths, yet were restricted to mapping with large (up to 20 m 2 ) pixels or to mapping communities of invasive shrubs rather than specific plants (Parker Williams and Hunt, 2002;Underwood et al., 2003). More recently Underwood et al. (2006) built upon this previous research and used hyperspectral imagery with an increased spatial resolution to detect individual plants in a scrubland delta in California, USA. Hyperspectral imagery has also been applied for invasive species detection in forested areas. Asner and Vitousek (2005) mapped tree invasions by analyzing hyperspectral imagery for nitrogen concentration and relating this concentration to invasive species presence across their study area in tropical forests of Hawaii. Recent studies mapping invasive plants with hyperspectral imagery include Narumalani et al. (2009) which used airborne hyperspectral imagery focused on the visible and near-infrared wavelengths to map various species along riparian areas, Bentivegna et al. (2012) which also used airborne hyperspectral imagery, but to map an invasive shrub along highways, Ishii and Washitani (2013) and Mirik et al. (2013) which both used airborne hyperspectral imagery to map invasive herbs in grasslands, Calviño-Cancela et al. (2014) which mapped invasive plants of various lifeforms in a coastal park, and Große-Stoltenberg et al. (2016) used imagery in the visible infrared to shortwave infrared wavelengths to map an invasive tree species in dune ecosystems. While hyperspectral imagery alone has shown success detecting invasive plants, knowledge of habitat characteristics may be important in some cases for adequate mapping (Andrew and Ustin, 2008).
Fusion of LiDAR and hyperspectral technologies has also been successful in mapping invasive species, as when combined they can be used to mask certain areas or provide additional contextual environmental information. Asner et al. (2008) used hyperspectral imagery to identify and map tree species in speciesrich Hawaiian rainforests and masked shaded pixels using LiDAR data. Chance et al. (2016) undertook a similar approach by masking densely forested and shaded areas using LiDAR data, to then map the distributions of two understory invasive shrub species on high spatial resolution hyperspectral images. Andrew and Ustin (2009) used hyperspectral imagery to detect an invasive plant and LiDAR data to create spatially contiguous potential habitat models from the sources of invasions, as detected by the hyperspectral imagery. This study combined LiDAR-derived metrics such as surface elevation, slope, aspect, and curvature with training field plots to produce the models. Gu et al. (2015) incorporated environmental information, and mapped gradients of species composition, including invasive species, in urbanrural gradients using hyperspectral imagery and LiDAR-derived structural variables, however, locations of individual species were not mapped.
Despite these advances in modeling and mapping invasive plant species distributions, the use of LiDAR data has been restricted to masks or structural and orographic variables when combined with hyperspectral imagery. Techniques developed for modeling direct and diffuse light regimes from LiDAR data in urban areas (see Tooke et al., 2012) may increase model accuracy by providing high resolution, relevant environmental context, in particular for plants that are sensitive to irradiance conditions (Metcalfe, 2005). Specifically, as light is a highly heterogeneous environmental condition that effects plant growth and interactions (Valladares, 2003), and as many invasive plants show preferences for certain light conditions (e.g., Metcalfe, 2005;Jansen et al., 2011), it may hold high predictive power in modeling invasive species distributions.
Another consideration when developing models with remotely sensed to detect plant distributions is deciding which type of model to utilize. While many ecological models can be developed using parametric or non-parametric regressions, they may not perform well when there are many relatively weak explanatory variables with complex interactions, as is often the case with remotely sensed data (Cutler et al., 2007). Random forests (RF) models, which are a collection of classification trees developed by selecting random subsets of the data in each, may solve this issue, as they are asserted to still be powerful when using many variables of various types (Breiman, 2001) and have been shown to outperform the majority of other classifiers in many modeling operations (Fernández-Delgado et al., 2014). Additionally, they have been used to map various invasive species with higher classification accuracies achieved compared to other models (e.g., Lawrence et al., 2006;Singh et al., 2015).
The aim of this study is to analyze the use of a combination of LiDAR data and hyperspectral imagery to map the distributions of two invasive plant species, Himalayan blackberry (Rubus armeniacus) and English ivy (Hedera helix), in the urban area of Surrey, British Columbia, Canada. First, a spectral detection algorithm using hyperspectral imagery was performed across open areas of the city to map the spectral closeness of each pixel compared to the two species (Figure 1). Second, LiDAR data was utilized to create a map of local solar irradiance conditions during the growing season in vegetated areas across the city.
Lastly, a map from a RF model was produced with LiDAR-derived environmental variables, the results of the spectral classification, and the irradiance models as inputs to map Himalayan blackberry and English ivy across the entire city at a 1.0 m pixel resolution. The accuracies of the maps were assessed, the importance of model variables and spatial information were analyzed, and implications were discussed.

STUDY AREA AND SPECIES
Surrey (49 • 11 N, 122 • 51 W) is located in the greater Vancouver area of British Columbia, Canada (Figure 2), covers 316 km 2 , and contains a mosaic of land cover patches including urban, agricultural, forested, wetland, and grassland area. Vegetated areas often contain dense, shrubby and herbaceous understories, as is commonly found in the temperate coastal climate of the area. Additionally, the city contains a mixture of public and private lands. On public lands, city management maintains an active parks and natural areas (Figure 2) system that contains many of the aforementioned land cover types. These land cover types are also present on private lands.
Himalayan blackberry (Rubus armeniacus) and English ivy (Hedera helix) are two invasive shrub species that land managers aim to control, due to their harmful effects on biodiversity, city residents, and users of natural areas. Himalayan blackberry displaces native vegetation and birds in both disturbed and undisturbed areas (Amor, 1973;Caplan and Yeakley, 2006;Astley, 2010), and alters successional patterns along streams (Fierke and Kauffman, 2006). Individual Himalayan blackberry plants grow and reproduce quickly, producing large amounts of seeds, and outcompete native plants (Schwartz et al., 1996;Caplan and Yeakley, 2006). Plants overgrow walking paths and roadways, creating a nuisance for city residents and managers. Himalayan blackberry generally grows near sea level in moist sites, in a variety of conditions from pastures to riparian areas and forest edges (Gaire et al., 2015). English ivy may produce allergens (Jones et al., 2009;Paulsen et al., 2010), has competitive advantages over native plants (Thomas, 1980;Dlugosch, 2005), and endangers users of natural areas by weakening trees (Thomas, 1980). English ivy also decreases light levels in the understory, resulting in reduced cover of native herbs and shrubs in the understory (Thomas, 1980;Harmer et al., 2001;Dlugosch, 2005), and may be a deterrent to wildlife (Barnea et al., 1993). It is challenging to control due to the many various life forms it can assume: herbaceous vines, climbers, herbs, woody shrubs, and sometimes small trees (Metcalfe, 2005). It is mostly a forestdwelling species, but is common in open areas as well (Clergeau, 1992).

Field Measurements
The presence and absence of Himalayan blackberry, English ivy, and other common plant species were recorded by City of Surrey field crews throughout 2012 and 2013 in park natural FIGURE 1 | Flow diagram showing the use of the LiDAR data, hyperspectral imagery, and field data to detect Himalayan blackberry (Rubus armeniacus) and English ivy (Hedera helix). SAM, spectral angle mapper.
areas. The municipal government sent experts to locate invasive species across the city according to a stratified design that targeted areas known to be susceptible to invasions first. These initial surveyed areas included forest edges, riparian areas, fields, edges to developed properties, and trails. Following the survey of these selected sites, other locations in parks and natural areas across the city were surveyed for invasions. Single points in large, discrete, established patches dominated by the invasive species with an area of greater than 40 m 2 separated by least 10 m from another plot were recorded as having either a presence or absence of the two species. These plots also contained a clear dominance of the invasive species over all species. An eTrex 10 GPS device (Garmin, Olathe, KS, USA) with a dGPS positional accuracy within 3 m was utilized. After the field survey, plots were randomly categorized as either training or validation, creating 646 training plots (108 Himalayan blackberry, 123 English ivy, and 415 of other species), and 496 validation (160 Himalayan blackberry, 114 English ivy, and 222 of other species) which were used for detection models.

Remotely Sensed Data
Between April 3 and April 11, 2013, airborne LiDAR data were collected for the entire City of Surrey (Figure 2) by Airborne Imaging (Calgary, AB, Canada). The Leica ALS70-HP (Leica Geosystems, Heerbrugg, Switzerland) system sent pulses with a wavelength of 1064 nm from an altitude of 1000 m with a view angle of 20 • and a scan frequency of 53 kHz resulting in a point density of 25 points/m 2 . TerraScan software (Terrasolid, Helsinki, Finland) classified points as ground or non-ground. TerraScan processing was completed using standard default parameters (Kraus and Pfeifer, 1998;Axelsson, 1999).
The Compact Airborne Spectrographic Imagery (CASI) 1500 acquired hyperspectral imagery over Surrey on May 4, 2013 at a 1.0 m pixel resolution. The imagery consisted of 72 spectral channels from 363 to 1051 nm with a constant channel width of 9.6 nm. Following the acquisition of the images, they were radiometrically corrected, georeferenced, and converted to reflectance imagery. To radiometrically correct the imagery, raw digital numbers, value representations of the contents of a pixel, were converted into spectral radiances for each spectral channel based on calibrations from ITRES Research (Calgary, AB, Canada), who calibrated the radiance imagery to known lab calibration files and corrected for dark and electronic offsets. These corrected images were georeferenced to a 1.0 m digital elevation model (DEM) produced from the LiDAR data and ground-based global positioning data (GPS) data. Lastly, atmospheric conditions and topographic and bidirectional effects were corrected to produce reflectance imagery that ranges from 0 to 100% using the ATCOR-4 procedure (Richter and Schlapfer, 2016).
Spectra of the two species of interest, Himalayan blackberry and English ivy, and other common shrub or herbaceous species in open areas in Surrey, lamium (Lamium galeobdolon) and grasses (lawn and European beachgrass [Ammophila arenaria]), were collected in April 2015 using a portable Analytical Spectral Devices (ASD) full range (FR) spectrometer (Analytical Spectral Devices, Boulder, CO, USA) within ±2 h from solar noon while the leaves were still attached. This spectrometer recorded spectra from 350 nm to 2500 nm at a 3 nm channel resolution up to 1000 nm and a 10 nm channel resolution at wavelengths between 1000 and 2500 nm. To account for differing light conditions between samples, the spectrometer was calibrated before each measurement with a Spectralon panel (Labsphere, North Sutton, NH, USA). At each measurement, the sensor was held approximately 10 cm from the surface of the plant, and leaves were manipulated to be stacked at least six deep, to achieve a standard optical depth (Datt, 1998;Jones et al., 2010). While 10 cm away from the plant, the sensor was moved around in circles to acquire spectra from various parts of the plant simultaneously, providing spectra that were representative of the entire plant from the view on the imagery. This sampling protocol controlled for variable irradiance conditions between samples by calibrating for the irradiance at each sample with the Spectralon panel and by limiting our sample area to the small area covered by the 10 cm circles. Between 10 and 15 samples were collected for each species. After measuring spectra, ASD channels were converted to the 72 CASI channels via averaging (Chance et al., 2016). Figure 3 shows the spectra of these species prior to resampling.

LiDAR-Derived Variables
As orography and forest structural characteristics may partially drive understory plant invasions (Royo and Carson, 2006), and as some invasive plant species, particularly shrub species, can alter forest structural characteristics (e.g., Hartman and McCarthy, 2008;Kuebbing et al., 2014), numerous orographic and forest attribute layers were created. Summaries of the vertical LiDAR data distribution were produced across the entire extent of the city, at a 1.0 m pixel resolution. The following layers were created using FUSION software (version 3.42, United States Forest Service, Seattle, WA, USA): a digital elevation model (DEM) of the elevation of the ground, a digital surface model (DSM) of the elevation of the first returns including trees and buildings, 75th (P75), 90th (P90), and 95th (P95) height percentiles, kurtosis of height, skewness of height, coefficient of variation of height, canopy cover and penetration above a height of 2.5 m, and cover below a height of 2.5 m ( Table 1). As height percentiles are a function of the locations of points vertical profiles, increased understory cover due to plant invasions may decrease the height percentile values when comparing forests with similar structure and height otherwise. Metrics that describe vertical point distributions were FIGURE 3 | Reflectance of Himalayan blackberry (Rubus armeniacus), English ivy (Hedera helix), grasses, and lamium (Lamium galedo. . .) in Surrey, BC, Canada as obtained from an Analytical Spectral Devices spectrometer in the spectral range of the hyperspectral imagery. Mean reflectance is indicated by dark lines, and ±1 standard deviation is indicated by the lighter lines.
included, as kurtosis of height was found to be negatively related to an invasive understory plant in Singh et al. (2015). For cover and penetration, a height of 2.5 m was chosen as a threshold as Himalayan blackberry and English ivy plants that are not vines on trees mostly occur below this height in Surrey.
Orographic variable layers were created in ArcMap (version 10.3, Esri, Redlands, CA, USA). Slope and aspect were calculated from the DEM. A cosine transformation was applied to aspect to convert the value range from −1 to 1, 1 being north (360 • ), and −1 being south (180 • ). Curvature rasters, which describe the concaveness of a surface, were created from the DEM. From the penetration layer, a measure of distance to open areas was created for each pixel. Pixels with penetration values above 90% were classified as open areas. Distances to these areas were calculated using Euclidean Distance. Lastly, a topographic wetness index, which is a unitless index quantifying hydrological accumulation based on slope and contributing area (Sørensen et al., 2006), was produced (Beven and Kirkby, 1979; Table 1).

Invasive Species and Land Cover Classifications
Rule images from a spectral angle mapper algorithm (SAM; Kruse et al., 1993) for Himalayan blackberry and English ivy were produced using the hyperspectral imagery in open areas (all areas with less than 20% canopy cover above a height of 2.5 m). SAM algorithms are well-established spectral matching algorithms that classify pixels based on the angle between spectral vectors. Rule images indicate the degree to which the spectral signature of a pixel matches the spectral signature of the target object. Smaller differences between angles represent more closely aligned spectra. Spectra from 73 of the field plots in open areas (38 of Himalayan blackberry and 35 of English ivy) were extracted from the hyperspectral imagery. As hyperspectral imagery can be processing intensive and contain redundancy in the spectral signals across channels (Clark et al., 2005), a spectral channel selection was applied to test whether or not all channels were needed for accurate classification. Chance et al. (2016) identified key spectral channels for differentiating Himalayan blackberry and English ivy from lamium and grasses, other common species in Surrey, and each other, using the instability index (ISI; Somers et al., 2010) channel selection and spectra from the ground-based acquisition. Table 2 describes these spectral channels. Prior to producing SAM rule images, forested areas were omitted by eliminating pixels with more than 20% cover above a height of 2.5 m according to the LiDAR-derived canopy cover layer. SAM rule images were produced by classifying the hyperspectral imagery using both the channel subsets and all channels of the spectra for Himalayan blackberry and English ivy across the open areas of the city in ENVI software (version 5.1, Exelis, McLean, VA, USA).
As land cover type may be an indicator for plant invasions (see Roy et al., 1999;Pearson et al., 2004), a land cover map of Surrey established by Plowright et al. (in review) was used as a categorical input in the detection models (Figure 1). Seven classes, coniferous forest, deciduous forest, grass, bare earth, paved areas, buildings, and water were classified at a 1.0 m pixel resolution using LiDAR data and hyperspectral imagery. The overall accuracy of this land cover classification was 88.6% (Plowright et al., in review). The resulting total surface area of pervious surfaces (coniferous forest, deciduous forest, grass, and bare earth) was 215 km 2 , 68% of the total area of Surrey. Topographic wetness index Wetness based on hydrological accumulation as modeled by the slope and contributing area, and as described by Beven and Kirkby (1979 All variables are rasters with a pixel size of 1.0 m except for the irradiance layers, which have pixel sizes of 3.0 m.

Irradiance Model
Because invasive plant distributions may be related to light levels, especially in forests in western North America (Parendes and Jones, 2000), an irradiance model indicating the light regime across the city was developed. Methods established by Tooke et al. (2012) were used to map irradiance at 3 m by 3 m pixel resolution. In this procedure, the LiDAR returns, the DSM, and the DEM were used in conjunction with solar angles and typical atmospheric conditions to calculate direct and diffuse irradiance on the 15th day of each month between March and October in 2013. Direct irradiance at each hour from 5 am to 10 pm local time was calculated according to the following three steps: atmospheric transmission, viewshed calculation, and vegetation transmission. Prior to taking these steps, the solar position was calculated according to the ENEA solar position algorithm (Grena, 2008) for each time interval.
The atmospheric transmission step calculates the irradiance that penetrates the atmosphere, a function of turbidity and Second, a viewshed of the effective horizon from the point 2.5 m above the DEM was produced for 36 viewing azimuths. At each pixel and azimuth, the angle of the tallest obstruction within 100 m, whether tree, building, or otherwise, was obtained from the DSM, relative to the orientation and slope of the pixel. This information was used to calculate the sky view factor (SVF), a measure of the directly visible proportion of the sky.
Within vegetated pixels, an extinction coefficient, calculated from the proportion of LiDAR returns at regular height intervals, was determined. A two-parameter Weibull distribution function, previously shown to be capable of characterizing vegetation structure from LiDAR returns (Coops et al., 2007;  2 | Key CASI (Compact Airborne Spectrographic Imagery) spectral channels and their corresponding wavelengths for differentiating Himalayan blackberry (Rubus armeniacus) and English ivy (Hedera helix) from each other and other common species in Surrey, BC, Canada, and the possible causes of the responses shown at these wavelengths and discussed in Chance et al. (2016).

Species
Wavelength at channel center (nm) Corresponding CASI channel number Possible reason for response as discussed in Chance et al. (2016) Himalayan blackberry 512 16 Nitrogen, phosphorus, potassium or a combination of these (Mutanga et al., 2004) 559 21 Lack of nitrogen or potassium (Mutanga et al., 2004) 655 31 Electron transition response to Chlorophyll a (Curran, 1989) Tooke et al., 2012), was produced at each using a vertical scaling parameter, a shape parameter, and the extinction coefficient as inputs. The vertical scaling parameter and shape parameter were estimated using a Levenberg-Marquardt leastsquares analysis for LiDAR returns in each pixel (Tooke et al., 2012). The extinction coefficient was the proportion of LiDAR returns above a certain height (Tooke et al., 2012). These parameters varied across the study area as they were determined by the vertical distribution of LiDAR returns at each pixel. The diffuse irradiance component was calculated using the diffuse transmissivity function and solar altitude function from Hofierka et al. (2002), which depends on atmospheric turbidity and solar elevation, in conjunction with the SVF (see Tooke et al., 2012, Eq. 9). The two output layers were the average daily diffuse and direct irradiance during the growing season from March to October with units of MJ m −2 day −1 .

Development
Four binary RF detection models were developed for each species in R programming environment (version 3.1.2, R Core Team) using the randomForest package (Liaw and Wiener, 2002) with 2000 trees and default values for the number of variables selected at each node and all other input parameters. Two models detected Himalayan blackberry and English ivy respectively in open areas using SAM rule images and LiDAR-derived variables relevant to open areas, and the other models detected the species in areas of canopy cover using only LiDAR-derived variables. In this present study, the model development processes were iterative, as initially all relevant remote sensing predictor variables ( Table 1) were considered. Using the training data, correlated and unimportant variables were iteratively removed until only uncorrelated variables with high predictive power remained, using mean decrease accuracy provided by the RF models as the measure of predictive power. Of variables that were correlated, the variable with higher predictive power was included in the model. Variables with negative mean decrease accuracies were excluded. Prior to fitting the model, two masks were created. Pixels with less than 10% cover below a height of 2.5 m were omitted and dirt, paved areas/rock, buildings, and water were masked using the land cover classification. The RF models were applied to the unmasked areas only, and produced maps of presence and absence of Himalayan blackberry and English ivy across the study area. For both species, the resulting maps of open areas and areas of closed canopy were combined to produce one map for each species or its distribution across the entire city.
According to Breiman (2001), RF classifications do not require independent accuracy assessments as some training data is inherently used as validation data to produce unbiased error. However, using this implicit error estimate may inflate the accuracy of the classification in some cases (Millard and Richardson, 2015). For detecting Himalayan blackberry and English ivy, the validation data from field plots were used to conduct an independent accuracy assessment. A minimum mapping unit (MMU) of 3 m radius circles around the validation plots was established since the two irradiance layers had pixels of 3 m resolution. Detection of the species of interest, either Himalayan blackberry or English ivy, within one of these plots counted as a presence in the accuracy assessment, indicating a true positive. Detection in an absence plot was recorded as a presence, and thus a false positive. The same logic applied for true and false negatives. Accuracies, kappa coefficients, and true skill statistics (Allouche et al., 2006) were assessed accordingly from the output map.
The uncertainty associated with accuracy was assessed as well at each validation point. RF models can provide spatial information about the probability of a pixel being classified as a certain class (either presence or absence in this present study; Cutler et al., 2007). This RF probability metric can be utilized to derive information about the classification uncertainty (Loosvelt et al., 2012). In this present study, U = 1 -p max was calculated for each pixel, where p max was the probability that a pixel was most likely to occur in its class based on the RF model. A higher value indicated class confusion and a lower value indicated greater certainty about the result (Loosvelt et al., 2012). In this present study, the maximum probability within the MMU that corresponded to the same class as the validation plot was used to assess the uncertainty.

Analysis
Variable importance was provided by the RF algorithm. The map was used to quantify the distributions of Himalayan blackberry and English ivy across the extent of the city and within parks. The total areas of Himalayan blackberry and English ivy were determined across the entire city and in parks according to the number of 1 m by 1 m pixels classified as the respective species. For analyzing the distribution of the species across the parks system, the number of parks containing the species was tabulated. Additionally, ecological relationships were analyzed by extracting information from the LiDAR-derived layers at all the field plots, training and validation. For this analysis, field plots were one of six categories: Himalayan blackberry open, Himalayan blackberry closed, English ivy open, English ivy closed, absence open, or absence closed. Values of the variables were compared between open and closed areas, as well as species presence and absence.
As the distribution of invasions across the city has management implications, clusters and outliers of Himalayan blackberry and English ivy were identified using Anselin Local Moran's I. Moran's I is an indicator used to quantify spatial autocorrelation (Moran, 1950;Anselin, 1995). Anselin Local Moran's I, which measures local autocorrelation, has proven to be an effective indicator of spatial clustering in ecological studies (see Fu et al., 2014;Swetnam et al., 2015). This metric ranks the abundance of species in each cell of a grid and determines which areas are significantly clustered or separated. A grid of cells 100 m by 100 m was applied across the entire study. At each cell in this grid, the number of pixels with Himalayan blackberry or English ivy invasion was tabulated, creating a measure of abundance. Anselin Local Moran's I was calculated across the city with α = 0.05 for both Himalayan blackberry and English ivy.

Models
The overall accuracies of the invasive species detection models ranged from 77.8 to 87.8% (Table 3). Himalayan blackberry in open areas was detected best, followed by English ivy in open areas (areas with canopy cover less than 20% above a height of 2.5 m), English ivy in areas with closed canopies (areas with greater than 20% canopy cover above a height of 2.5 m), and Himalayan blackberry in areas with closed canopies ( Table 3). Detection of Himalayan blackberry in open areas had the highest kappa coefficient (0.75) followed by English ivy in open areas (0.56), Himalayan blackberry in areas with closed canopy (0.55), and English ivy in areas with closed canopy (0.53; Table 3). True skill statistics took a different order, with Himalayan blackberry in open areas receiving the highest (0.78), followed by Himalayan blackberry in areas with closed canopy (0.55), the English ivy in open areas (0.52) and areas with closed canopy (0.48; Table 3). Himalayan blackberry was detected with higher true positive rates (sometimes referred to as sensitivity) than English ivy in both open areas and areas with closed canopies ( Table 3). English ivy was detected with higher true negative rates (sometimes referred to as specificity) than Himalayan blackberry in both open areas and areas with closed canopies (Table 3). False negative rates (miss rates) were higher for English ivy than for Himalayan blackberry and false positive rates (fall-out) were higher for Himalayan blackberry than English ivy ( Table 3). All classifications had low uncertainty, as the majority of validation points showed uncertainty values below 0.25 (Figure 4). Himalayan blackberry in areas with closed canopies had the highest uncertainty (Figure 4).
In open areas, the SAM rule images were ranked as highly important for both Himalayan blackberry and English ivy (Figure 5). For both species, the SAM rule image produced using the subset of spectral channels was chosen over the SAM rule image from all spectral channels (Figure 5). Direct irradiance was more important to English ivy detection than to Himalayan blackberry detection in open areas (Figure 5). In areas with closed canopies, the most important predictor variables -75th and 95th percentile height, skewness, and the coefficient of variation of height -were related to forest structural characteristics (Figure 5). Aspect was an important variable in both open areas and areas with closed canopies for both species (Figure 5). Land cover classification was ranked as the least important variable in all classifications, however, its product, distance to impervious surface, and was ranked highly for detecting Himalayan blackberry in open areas (Figure 5). Figure 6 shows examples of the final classification results.

Distribution of Invasive Species
Across the city, Himalayan blackberry covered 1.18 km 2 and English ivy covered 1.51 km 2 , corresponding to 0.5 and 0.7% of the pervious areas respectively ( Table 4). Anselin Local Moran's I results showed that Himalayan blackberry and English ivy were similarly clustered, with 98.5 and 99.3% of 100 m by 100 m grid cells containing either of the species located in clusters ( Table 5). 1.5% of grid cells with Himalayan blackberry and 0.7% of grid cells with English ivy were spatial outliers ( Table 5). Within parks, Himalayan blackberry covered 0.16 km 2 and English ivy covered 0.35 km 2 , corresponding to 0.8 and 1.8% of the pervious surface respectively ( Table 4). Although the area covered by English ivy was greater than Himalayan blackberry in parks, Himalayan blackberry was present in more parks (Table 4).  In addition to differing in their spread across the city and parks systems, Himalayan blackberry and English ivy invasions differed in their relationships to forest types, irradiance, and roads. Himalayan blackberry was found more often in deciduous areas with closed canopies than in coniferous areas with closed canopies (Figure 7). The opposite was true of English ivy. In both open areas and areas with closed canopies, Himalayan blackberry did not have a strong relationship to direct irradiance compared to absence plots, whereas English ivy occurred at areas with less direct irradiance (Figure 8). Himalayan blackberry and English ivy in areas with closed canopies followed similar trends in regards to the cosine of the aspect (Figure 8). Both species had higher relative frequencies at most south facing aspects than absence plots, and vice versa for north facing aspects (Figure 8). For both species in open areas, there were relatively fewer occurrences at extreme north and south aspects compared to absence plots, yet there were relatively more occurrences at east or west aspects when compared to absence plots (Figure 8). For both species in both areas, the majority of profile curvature values occurred between 0 and 20, which corresponds to slightly concave surfaces (Figure 8). Topographic wetness for both species in both areas were skewed toward drier values (Figure 8). Slope was skewed toward lower values; however, English ivy in open areas and areas with closed canopies occurred most frequently at slopes between 5 and 10 degrees, whereas Himalayan blackberry and absences occurred mostly at slopes between 0 and 5 • (Figure 8). Himalayan blackberry had a peak coefficient of variation of height at 0.8, whereas English ivy had a peak at 0.7 (Figure 8). English ivy had a coefficient of variation of height more skewed toward smaller values than Himalayan blackberry (Figure 8). Both presences and absences for both species in open areas were closer to impervious surfaces than in areas with closed canopies (Figure 8).
The absolute amount of Himalayan blackberry and English ivy decreased as the distance from roads increased, yet when this area was considered proportionally to the land area at various distances from roads, there were no relationships between invasive species occurrence and the distance to a road at areas less than 500 m from a road (Figure 9). Areas further than 500 m away from road had proportionally higher rates of invasion than areas within 500 m of a road (Figure 9).

DISCUSSION
Few studies have mapped the distribution of invasive shrub species in complex urban areas using remotely sensed products. In one of these studies, Singh et al. (2015) mapped an invasive shrub, Chinese privet (Ligustrum sinense), in Charlotte, North Carolina at a 5 m by 5 m resolution using LiDAR data and a RF classifier and achieved overall accuracies between 81 and 89%. The results of this present study were similar to those found in Singh et al. (2015) with overall classification accuracies of 77.8 and 87.8% ( Table 3). The true skill statistic, a metric of classification performance that is independent of prevalence, which the kappa coefficient and overall accuracy are not (Allouche et al., 2006), showed that Himalayan blackberry was actually better classified than English ivy in all cases due to the low true positive rates in the English ivy classifications ( Table 3). The relatively low true positive rates for English ivy in both open areas and areas with closed canopy indicate that the RF classifier may have underestimated the area covered by English ivy. The low true positive rates may also indicate that the classifier confused English ivy with other species surveyed in Surrey, and the estimated area may not be inaccurate. However, English ivy is more specialist than Himalayan blackberry (Metcalfe, 2005;Gaire et al., 2015), thus the second possibility is less likely considering the high accuracies of the Himalayan blackberry classifications. Additionally, the uncertainty analysis showed that Himalayan blackberry classifications actually had higher uncertainties that the English ivy classifications (Figure 4), indicating that instead of confusing English ivy with other species the RF classifier was likely to underestimate the species range for other reasons. As English ivy covered more areas in parks and across the city than Himalayan blackberry, but was present in a proportionately smaller number of parks (Table 4), English ivy was likely more clustered than Himalayan blackberry. Anselin Local Moran's I results did not corroborate this, however, as the number of 100 m by 100 m grid cells that were considered parts of Himalayan blackberry or English ivy clusters were similar for the two species (Table 4). It is logical that English ivy was clustered when considering its habitat preferences. As English ivy prefers coniferous forests (Clergeau, 1992 ; Figure 7), its potential range was limited to a few patches across the city.
Open areas were classified with higher accuracies than areas with closed canopies for both Himalayan blackberry and English ivy. This is likely due to the inclusion of the rule images from the hyperspectral imagery in the model. The SAM rule images provided a continuous variable on which the RF model could make decisions. Other SAM classifications, such as that classifying Himalayan blackberry and English ivy in Chance et al. (2016), transform the continuous rule image into a binary classification image by defining a threshold SAM value, below which is classified as presence and above which is classified as absence. By keeping the classification results as continuous, the RF model had greater flexibility in defining this threshold based on the values of the other predictor variables. For both Himalayan blackberry and English ivy, the SAM rule images were highly ranked as important (Figure 5). More specifically, the SAM rule images produced from the subset of spectral channels were chosen for the final models rather than those produced with all channels as they were ranked as more important. This is consistent with previous results comparing detection accuracies of different channel selections. For example, Cho et al. (2010) used an SAM-based spectral channel selection process to increase overall detection accuracies of savanna tree species. Naidoo et al. (2012) used RF-and SAM-based spectral channel selections and showed that structural information combined with spectral information increases classification accuracy of tree species in RF models. More recent studies utilizing a subset of spectral channels for classification include Somers and Asner (2013) and Kuo et al. (2014), both of which showed an increase in classification accuracies using channel subsets. Additional considerations for the hyperspectral imagery are the phenology of deciduous trees in the study area and the timing of the imagery, as leaf-off conditions for deciduous trees while understory invasive plants have leaves present increases the area for spectral analyses. However, as hyperspectral imagery is often collected to serve multiple purposes, catering the acquisition timing to one or even several species may be difficult. Thus, satellite imagery programs have been well-suited to mapping invasive understory plants in deciduous forests due to their multiple acquisition dates (Bradley, 2014). For example, Resasco et al. (2007), Wilfong et al. (2009), andShouse et al. (2013) all used Landsat sensors during fall or winter conditions to map the understory invasive plant, honeysuckle (Lonicera maackii) across different environments. Kimothi et al. (2010) map Lantana camara during leaf-off conditions in India. While these studies successfully mapped plant invasions, the canopies in the forests were dominated by deciduous trees. Such an approach would not cover as wide as breadth in coniferousdominated forests. However, there is an opportunity to explore utilizing hyperspectral imagery over the deciduous portions of this forest in leaf-off conditions to detect understory invasive plants.  In this present study, LiDAR data effectively modeled the species distributions even without combining it with hyperspectral data, highlighting the applicability of LiDAR data to be a source of information for ecological models. While this has been shown in the past using orographic and structural information derived from LiDAR data, this present study added LiDAR-derived irradiance to the ecological models. In terms of the use of specific LiDAR products, direct irradiance was more highly ranked than diffuse irradiance, hence it being included in the final models. Direct irradiance was also highly ranked for English ivy detection in open areas (Figure 5). Consistent with the known habitat preferences of the English ivy for cooler, shaded areas, and of Himalayan blackberry for open areas (Clergeau, 1992;Gaire et al., 2015), Himalayan blackberry preferred areas with more direct irradiance than English ivy (Figure 8), indicating that when in open areas, English ivy still established in shaded locations. Direct irradiance was likely important for determining the areas of English ivy in open areas due to the limited sites at which direct irradiance conditions suit English ivy in open areas. Forest structural characteristics were more important than orographic variables in areas with closed canopies for both species (Figure 5). Specifically, measures of the vertical distribution of LiDAR points, the coefficient of variation of height and skewness, were important predictors ( Figure 5). Based on Figure 8, English ivy had a different preference for coefficient of variation than absent plots. This indicates that English ivy either established in sites with certain structural characteristics, influences the structure, or both. Previous research indicates that English ivy plants prefer shaded forests with trees with large trunks (Schnitzler and Heuzé, 2006), which in Surrey may be forest stands with relatively tall, older trees of similar ages. Additionally, as English ivy can decrease the amount of understory cover by other species (Biggerstaff and Beck, 2007), many of which in Surrey are taller shrubs; the standard deviation of height of these stands may decrease.
While Himalayan blackberry and English ivy distributions may be more related to forest structural characteristics than orographic ones, the pixel resolution may influence this signal. Gu et al. (2015) used a pixel size of 20 m to map forest composition with LiDAR-derived forest structural variables and hyperspectral imagery, and found that the imagery was more important for their purposes, perhaps due to the limitations caused by the large pixel size. Numerous studies have shown that increasing the coarseness of the spatial resolution of hyperspectral imagery can decrease its applicability for detecting species (e.g., Underwood et al., 2007;Peña et al., 2013), however, few studies have quantified the effects of LiDAR footprint sizes or point densities for ecological mapping of single understory plant species. A review of some relevant studies indicates that point densities of 25 points per m 2 sampled to 1 m pixels, as are present in this current study, may be at a finer scale than needed for modeling or management purposes as follows. One study focusing on estimating biomass of one understory species from LiDAR data found that pixel sizes between 5 and 30 m produced better models than those at 1 m (Li et al., 2015). As previously mentioned Singh et al. (2015) mapped understory plant invasions at a 5 m pixel size. Barber et al. (2016) used LiDAR data with 1.4 points per m 2 sampled to 25 m pixels to estimate understory shrub abundance based on ecological relationships. Furthermore, Brubaker et al. (2013) tested the accuracies of LiDAR-derived DEMs and differences between derived products such as curvature, slope, and roughness and found no substantial differences between DEMs produced from LiDAR data with more than 10 points per m 2 and less than 1 point per m 2 . As evidenced by these varying results, point densities, and pixel sizes, more research is needed to quantify ideal LiDAR resolutions for different applications. Moreover, in urban areas where LiDAR data is often collected for multiple purposes (e.g., engineering, flood control, forestry), the data may need to be downgraded after acquisition to suit certain purposes.
The land cover classification, a product of both the hyperspectral imagery and the LiDAR data, was ranked as the variable with the lowest importance in all models (Figure 5). Because water, paved areas, and buildings were masked, only coniferous forests, deciduous forests, and areas of grass were considered in the RF model. The spatial breadth of these classes may have contributed to the low importance of the land cover classification, as other variables could likely better explain the spatial distributions of the invasive species. However, analysis of the locations of the invasions found that Himalayan blackberry and English ivy differed in abundance between deciduous and coniferous forests (Figure 7), with Himalayan blackberry showing a preference for deciduous forests and English ivy for coniferous forests, consistent with the habitat preferences of the two species (Clergeau, 1992;Gaire et al., 2015).
Analysis of the distance between roads and invasion showed that both absolute coverage of Himalayan blackberry and English ivy was inversely related to the distance from roads (Figure 9). These results have management implications, as they indicate target areas for city managers to direct resources toward curbing the spread of invasion. Additionally, targeting areas near roads may be relatively low effort as field crews can access invasion patches without traveling far into urban natural areas. However, when considered proportionally, the area covered by Himalayan blackberry and English ivy was not related to the distances from roads (Figure 9). This contradicts the general trend shown in previous research on English ivy ecology, which shows decreasing proportional cover of English ivy with increasing distance from roads (Arévalo et al., 2008). However, this previous study differs from the present study in that it did not occur in an urban area; the forest was made up of mostly deciduous trees, and only analyzed species composition up to 60 m from roads. The distance to impervious surfaces, while not directly analogous to distance to roads, is also indicative of human activity and is somewhat important for predicting Himalayan blackberry distributions in open areas and areas with closed canopies according to Figure 5. Figure 8 shows that occurrences in open areas are more frequently closer to impervious surface than in areas with closed canopies, but this is likely due to there being more open area closer to impervious areas. One reason that impervious surfaces may not be a strong indicator of invasions in this study are that impervious surfaces relate better to the beginning of invasions, but once invasions are established, the relationship weakens (Beauséjour et al., 2015). Another reason may be that because this study is in an urban area, all areas are close enough to roads and impervious surfaces not enough variation exists to observe a substantial effect. Additionally, both roads and impervious surfaces aid in dispersal, but there may be other variables that ultimately lead to the invasion remaining on the landscape (Beauséjour et al., 2015).
This present study highlights a trend in ecological and remote sensing research of combining the two disciplines. Often remote sensing scientists create data products without consideration for how the products are going to be used. In these cases, ecologists must find ways to cater these products to their research. This present study shows that processing remotely sensed data specifically for ecological modeling of certain plant species is effective. Therefore, land managers and researchers modeling invasive species distributions should consider the goals of the modeling in conjunction with the spectral and spatial resolutions of remotely sensed data before processing it.

CONCLUSION
Current municipal operational approaches for detecting plant invasions rely on field crews that cannot produce spatially contiguous information across a large area about the distribution of plant species. Remote sensing technologies can augment this approach by detecting plants across large urban areas with a single methodology. This study produced detection models for Himalayan blackberry and English ivy across Surrey, BBC, Canada. RF models were used in conjunction with variables from hyperspectral imagery and LiDAR data to detect Himalayan blackberry and English ivy invasions across open and forested areas of the city. Spatial relationships on the resulting maps were quantified by analyzing clusters and relationships with environmental variables. The resulting RF detection models classified Himalayan blackberry with 78.8 and 87.8% accuracies and English ivy with 81.9 and 82.1% accuracies, with plants in open areas being better detected than those in areas with closed canopies. Spatial analysis of the resulting maps showed that Himalayan blackberry occurred more often in deciduous forests whereas English ivy tended to occur in coniferous forests. Both species were found to have a negative relationship with distance from roads when considering absolute area. This study highlights the applicability of LiDAR data and hyperspectral imagery in mapping plant species, specifically invasive plant species distributions in open areas and areas with closed canopies. This study also demonstrates that creating LiDAR-and hyperspectral-derived products specific to species to is an effective way to model their distributions. Furthermore, this present study shows that novel LiDARderived products, such as irradiance models, may be useful for ecological models over large extents. Some limitations of the study include its reliance on data from public land and its applicability to only mature plants. Future research could address these limitations by broadening the training datasets to include private property and plants of different life stages.

AUTHOR CONTRIBUTIONS
CC: Wrote and compiled the manuscript, performed analyses and data processing. NC: Wrote and edited parts of the manuscript. advised on which analyses to FIGURE 9 | (A) Proportion covered by Himalayan blackberry (Rubus armeniacus) and English ivy (Hedera helix) and the total area covered in the city related to the distance from roads, and (B) total area covered by Himalayan blackberry (Rubus armeniacus) and English ivy (Hedera helix) related to the distance from roads in Surrey, BC, Canada.
perform and how to process to the data. AP: Processed data for the land cover classification section and wrote parts of the results and methods that corresponded to the land cover classifications. TT: Provided code, advice, and edits related to the irradiance models. AC: Provided edits and suggested analyses related to the irradiance models and random forest classification results. NA: Provided data used in this paper and manuscript edits.

ACKNOWLEDGMENTS
The City of Surrey and the Natural Sciences and Engineering Research Council of Canada (NSERC) through Engage and Engage + grants supported the research. The authors thank assistants Yuhao Lu, Joanna Lee, and Ignacio San-Miguel. Additional thanks to Dr. Jeanine Rhemtulla and Ken Crosby for their insightful suggestions and comments. Further thanks to the three reviewers for their knowledgeable and informative advice.