Detecting Long-Term Urban Forest Cover Change and Impacts of Natural Disasters Using High-Resolution Aerial Images and LiDAR Data

Urban forests provide ecosystem services; tree canopy cover is the basic quantification of ecosystem services. Ground assessment of the urban forest is limited; with continued refinement, remote sensing can become an essential tool for analyzing the urban forest. This study addresses three research questions that are essential for urban forest management using remote sensing: (1) Can object-based image analysis (OBIA) and non-image classification methods (such as random point-based evaluation) accurately determine urban canopy coverage using high-spatial-resolution aerial images? (2) Is it possible to assess the impact of natural disturbances in addition to other factors (such as urban development) on urban canopy changes in the classification map created by OBIA? (3) How can we use Light Detection and Ranging (LiDAR) data and technology to extract urban canopy metrics accurately and effectively? The urban forest canopy area and location within the City of St Peter, Minnesota (MN) boundary between 1938 and 2019 were defined using both OBIA and random-point-based methods with high-spatial-resolution aerial images. Impacts of natural disasters, such as the 1998 tornado and tree diseases, on the urban canopy cover area, were examined. Finally, LiDAR data was used to determine the height, density, crown area, diameter, and volume of the urban forest canopy. Both OBIA and random-point methods gave accurate results of canopy coverages. The OBIA is relatively more time-consuming and requires specialist knowledge, whereas the random-point-based method only shows the total coverage of the classes without locational information. Canopy change caused by tornado was discernible in the canopy OBIA-based classification maps while the change due to diseases was undetectable. To accurately exact urban canopy metrics besides tree locations, dense LiDAR point cloud data collected at the leaf-on season as well as algorithms or software developed specifically for urban forest analysis using LiDAR data are needed.


Introduction
Ecosystem services measure the urban forest's benefits and economic value in four categories: supporting (e.g., biodiversity), regulating (e.g., air quality), cultural (e.g., health), and provisioning (e.g., fresh water) [1]. Measuring urban tree canopy cover and other metrics, such as size and height, provides the most basic quantification of the urban forest's potential ecosystem services. Currently, the most accurate way to collect tree attribute data is by trained personnel collecting individual tree metrics using ground assessments within an urban forest. However, ground assessments are both time-consuming and expensive, and it is not physically possible or legal to access every tree within an urban environment.
The coordinates of its central location are 44.323°N and 93.958°W. The total area of the City is 14.76 km 2 with 11,400 residents [42]. The total boulevard tree population is 4824 trees. The majority of the population consists of 1414 maple (Acer spp.), 1094 ash (Fraxinus spp.) and other species such as 404 basswood (Tilia spp) and 301 hackberry (Celtis spp.) [43]. On March 29, 1998, between 17:18hrs and 17:36hrs, a tornado struck the city, and many of the town's houses, businesses, civil buildings and urban forests were damaged. Therefore, one aspect of this study examined the impact of the 1998 tornado on the urban canopy cover area. The possibility of detecting forest canopy change due to the tree diseases Dutch Elm Disease (Ophiostoma ulmi and O. nova-ulmi) and Butternut Canker (Sirococcus clavigigenti-juglandacearum) was also explored, as hundreds of millions of trees have been infected since the 1960s in the USA.

Data and Preprocessing
The complete methodological process, described in Sections 2.2-2.6, is shown as a schematic in Figure 2. The research data consisted of digital airborne ortho-images and aerial photographs,

Data and Preprocessing
The complete methodological process, described in Sections 2.2-2.6, is shown as a schematic in Figure 2. The research data consisted of digital airborne ortho-images and aerial photographs, hardcopy historical aerial photographs, scanned images and historical maps, vector data, and LiDAR data. Tables 1 and 2 provide detailed information regarding the multi-source datasets. Seven selection criteria determined the aerial images used: (1) high image resolution, i.e., 600 Dots Per Inch (DPI) or pixel resolution ≥ 0.93 m; (2) minimum map scale of 1:20,000; (3) trees with leaves on; (4) historical photographic images taken during the same months; (5) images with the full areas within the city boundary of 1928 and 2017; (6) ease of access to and availability of images; (7) images just prior to and post the 1998 tornado. The images selected were for 1938, 1951, 1964, 1995, 2008, and 2017 (Table 1).
Remote Sens. 2020, 12, x FOR PEER REVIEW 4 of 21 data. Table 1 and Table 2 provide detailed information regarding the multi-source datasets. Seven selection criteria determined the aerial images used: (1) high image resolution, i.e., 600 Dots Per Inch (DPI) or pixel resolution ≥ 0.93 m; (2) minimum map scale of 1:20,000; (3) trees with leaves on; (4) historical photographic images taken during the same months; (5) images with the full areas within the city boundary of 1928 and 2017; (6) ease of access to and availability of images; (7) images just prior to and post the 1998 tornado. The images selected were for 1938, 1951, 1964, 1995, 2008, and 2017 (Table 1).    * Data collected in leaf-off periods, April 8-May 5 and November 2-19, 2010, via fixed-wing aircraft equipped with LiDAR system (Optech Gemini) including differential GPS unit and inertial measurement system. Area data horizontal positional accuracy: acquired at or below 1700 m above mean terrain with a horizontal accuracy of 2.54 cm. Vertical positional accuracy values (RMSE): better than 10 cm. Fundamental vertical accuracy of the classified bare earth: 0.08 cm at 95% confidence level in the 'open terrain' land cover category. Diminished accuracy expected in areas of dense vegetation due to fewer points defining bare earth in those areas [44].
To prepare the images, first, the 1938 and 1964 Black and White (B+W) panchromatic images were cropped to remove black borders and distorted sections of the image deemed unsuitable for analysis, e.g., incorrect tone. Second, the 1995 aerial image and the 1928 engineering map were transferred from hard copy to soft copy by digital scan (3-band, natural color at 600 DPI). Third, the 1938Third, the , 1951Third, the , 1964, and 1995 aerial photographic images lacked camera parameters and exterior orientation parameters [2] and could not be orthorectified. Therefore, the photographic images and the 1928 engineering map were georeferenced [2] within ArcGIS 10.6.
The root-mean-square errors (RMSEs) range from 1.40 m to 113.60 m; on average, 17 GCPs were used. The 1951 south image and 1928 engineering map (Table 3) had relatively high RMSEs; both were georeferenced with a first-order polynomial transformation. The 1928 GCPs were based on a reference photograph from 1951, due to difficulty in determining viable GCPs in 1938. The high RMSE is likely due to the limited GCPs available, which were based on the road layout; it is also probable that the road layout changed between 1928 and 1951. For example, the Saint Peter Broadway Bridge (highway 99) was remodeled three years after the 1928 engineering map [45]; therefore, the road layout detailed in 1951 may have changed accordingly. It should be noted that the process for accuracy of individual images was based on a compromise between both RMSE and visual interpretation due to the disparity in image type, location of suitable GCPs for overlapping and adjacent images, limited availability of GCP sites, and temporal variation of the photographic images when georeferencing. Therefore, low RMSEs do not guarantee high positional accuracy [46]. After georeferencing, the city limits were digitized and a polygon feature class created from the 1928 engineering map. The 1928 and 2017 feature class boundaries were buffered to 60.96 m for two reasons. First, the city boundary is a human construct and the urban forest does not conform to a linear boundary. Using visual estimation, it was determined that no trees in the boundary areas had a canopy width greater than 121.92 m; therefore, the 60.96-m buffer contained trees within the city boundary. Second, the 60.96-m boundary and the remaining minimum bounding rectangle with the 1928 and 2017 city boundary was a large enough area to perform validation results on the aerial photographic images. The 1998 tornado boundary was created by buffering two polyline shapefiles to a distance of 2.01 km, clipping them to the 1928 and 2017 boundary, reflecting the geographic extent of the tornado as described by the NOAA Storm Prediction Center's (SPC) United States severe report database [47].
The aerial photographic images were clipped to the 1928 and 2017 city feature classes and a single raster dataset was created by mosaicking the clipped aerial photographic images. To aid distinction between land use classes for the B+W panchromatic and 1995 three-band color image, a 3 × 3 variance (second-order) texture layer was created from the mosaicked photographic images within ERDAS IMAGINE 2018 [48]. Each texture layer was stacked with the mosaicked photographic image.

OBIA Urban Forest Canopy Analysis and Change Detection
For analyzing high-resolution B+W panchromatic aerial photographs, texture can be crucial to maximize extraction of land use categories, particularly forest and urban forest areas [48][49][50]. B+W Panchromatic photographs have low radiometric and spectral information and texture adds another variable for extraction [48]. Feature Analyst, the software chosen for OBIA analysis, includes texture as an object input [51]; the combination of the texture layer created in ERDAS IMAGINE 2018 and that included within Feature Analyst gave the best results.
Feature Analyst uses a contextual classifier in its segmentation process which utilizes object size, edge type, spatial context, and shape to produce vector files that can be edited [51,52]. After preliminary investigations using this software, it became apparent that the extraction classification process was most effective if supervised learning using multiple classes (Table 4) was used rather than unsupervised or a single-class approach [53]. Large shrubs were included within the urban forest, as it is virtually impossible to differentiate between trees and large shrubs in the B+W panchromatic images due to the resolution and hue/tone. It also became apparent, as documented by Yuan [48], that differentiation between grass and soil in the B+W panchromatic images was challenging; the classes were combined as grass/soil. For consistency, these classes were kept for B+W panchromatic, three-band natural color, and four-band color infrared images. All the images had some form of shadow; therefore a shadow/other class was created as detailed in Qiu, Wu [54]. Training samples were manually digitized drawing close to the edge of class features (objects) to represent the shape, size, spectral content, texture, patterns, and contextual data and to distinguish them from adjacent objects [54,55]. Due to the differences between the three types of aerial photographic images, multiple supervised learning operations were repeated using different histogram stretches, learning options, and input representations to extract each class.
Post-processing of the OBIA data was performed first in ArcGIS by visibly comparing the OBIA canopy cover class polygons to the photographic images and removing inaccuracies. Second, the OBIA polygon was converted to raster and imported into ERDAS IMAGINE 2018, where a majority function using a 3 x 3 window was used to remove single pixels incorrectly classified as canopy cover.
To quantitatively determine and compare canopy change between the years, a post-classification change detection model was created within ERDAS IMAGINE 2018 [56]. The change detection model created a thematic layer from a matrix that evaluated two historical image files. This thematic layer displays the unique difference in the values of the two overlapping original images [57].

Accuracy Assessment of Urban Forest Canopy cover Extracted by OBIA
For accuracy assessment, a site-specific thematic error matrix [58] and Kappa coefficient were created, using the method proposed by Zhou and Troy [59]. Stratified random sampling enabled the canopy and non-canopy classes to be proportionately weighted; there were fewer canopy class features compared to non-canopy features and these were unevenly distributed [54]. To mitigate against potential positional inaccuracy, a small difference in cell size between photographic images of different years, and ensure a representative sample within a category polygon [58], a sample unit of a 5 × 5 pixel block was created for each year. Congalton and Green [58] state that when the research site is less than 404,686 hectares and 12 classes, a minimum of 50 samples per class is needed for accuracy assessment. Therefore, 150 samples were used: 50 per class (forest, non-forest) and 50 randomly distributed. As no higher-resolution images were available for reference data, the same high-spectral-resolution aerial photographic images were used for the accuracy assessment process. Field sample collection was possible for the 2017 dataset. The area selected for field sample data was based on the ability to legally and practically assess areas within the 2017 city boundary: the city owned land (parcel) and the rights of way (ROW) of all roads within the city boundary ( Figure 3).
The field sampling map was created in ArcGIS by buffering the Nicollet County streets shapefile to reflect the correct ROW distance for streets within the city, dependent on municipal code. This dataset was merged with the city tax parcel polygon shapefile. Any significant difference between the 2017 photographic image and 2019 ground truthing was mitigated by the researcher's in-depth knowledge as the city's urban forester.
Remote Sens. 2020, 12, x FOR PEER REVIEW 8 of 21 Remote Sens. 2020, 12, x; doi: FOR PEER REVIEW www.mdpi.com/journal/remotesensing based on the ability to legally and practically assess areas within the 2017 city boundary: the city owned land (parcel) and the rights of way (ROW) of all roads within the city boundary ( Figure 3). The field sampling map was created in ArcGIS by buffering the Nicollet County streets shapefile to reflect the correct ROW distance for streets within the city, dependent on municipal code. This dataset was merged with the city tax parcel polygon shapefile. Any significant difference between the 2017 photographic image and 2019 ground truthing was mitigated by the researcher's in-depth knowledge as the city's urban forester.

Urban Forest Canopy Cover Assessment Using i-Tree Canopy
i-Tree is a freely available suite of urban forest ecosystem service evaluation software, created by the U.S. Forestry Service and other collaborators [60,61]. i-Tree Canopy is one tool within i-Tree that allows users to define land cover categories (e.g., water, trees) using a web browser application (Google Maps and Google Earth) to enable photo interpretation of current and historical aerial imagery using randomly selected points [62]. The Google Maps 2019 image and the Google Earth 2008 image were selected for their relevance.
Both 1928 and 2017 polygon boundary shapefiles were imported into i-Tree Canopy and five classes, as per OBIA, were created. Employing the 2019 Google photographic image in Google Maps;

Urban Forest Canopy Cover Assessment Using i-Tree Canopy
i-Tree is a freely available suite of urban forest ecosystem service evaluation software, created by the U.S. Forestry Service and other collaborators [60,61]. i-Tree Canopy is one tool within i-Tree that allows users to define land cover categories (e.g., water, trees) using a web browser application (Google Maps and Google Earth) to enable photo interpretation of current and historical aerial imagery using randomly selected points [62]. The Google Maps 2019 image and the Google Earth 2008 image were selected for their relevance.
Both 1928 and 2017 polygon boundary shapefiles were imported into i-Tree Canopy and five classes, as per OBIA, were created. Employing the 2019 Google photographic image in Google Maps; one thousand randomly generated survey points were produced, and, through photo-interpretation, each point was assigned a class, as recommended by the i-Tree Canopy technical notes [63]. To photo-interpret the 2008 Google Earth image, a Keyhole Markup Language (KML) file, representing the randomly generated survey points produced for the 2019 images, was imported into Google Earth. Using the 2008 Google Earth photographic image, each imported point was photo-interpreted and assigned a class.

Urban Forest Metrics Detection Using LiDAR
LiDAR data was analyzed with LiDAR360 software [64]. To improve data's quality, LiDAR data preprocessing included reclassifying the point cloud and removing outliers. LiDAR360 can re-classify by machine learning, using a small training sample that was manually corrected to reclassify points and then using the random forest method along with the training sample to edit the entire dataset in batches [64]. Initial results were poor; therefore, the following three LiDAR point classification processes were applied and proved successful. First, bare earth ground points were classified/reclassified using the improved progressive TIN densification filtering algorithm proposed by GreenValley International Ltd [64], Zhao, Guo [65]. The maximum building size was set to 160 m, the actual length of the largest building within the study site. Second, the bare earth ground points were refined by fitting quadratic surfaces using a concoid filter [64]. The point cloud was classified, paying particular attention to the vegetation points using the interactive editing module, which allows the user to manually edit points or groups of points with a profile tool. Finally, to remove topographic relief elevation effects, the LiDAR points were normalized by subtracting the closest classified bare earth ground point elevation from other classified points' z value [64].
To create a canopy height model (CHM), a Spike Free TIN method [64,65] was used to interpolate a Digital Surface Model (DSM) and a Digital Elevation Model (DEM). Medium vegetation returns were used instead of first returns to improve accuracy by removing buildings, etc., which in an urban environment can be higher than the urban canopy. The CHM was created by subtracting the DSM from the DEM. Canopy cover density was determined by calculating the ratio of medium vegetation returns to the total number of LiDAR first returns for a pixel [64]. Vegetation below 2 m was removed, and the pixel size was determined by measuring the width of the largest single tree canopy, 33 m.
Two different tree segmentation models-CHM segmentation and the point cloud segmentation model-were used to determine individual tree attributes, e.g., tree location, height, crown diameter, crown area and volume. In CHM segmentation, watershed segmentation recognizes and demarcates individual tree crowns using the CHM created in the previous step as the input layer. The watershed segmentation algorithm inverts individual tree canopy points to represent a catchment basin. The level at which the water would start to overflow the basin represents the bottom of the individual tree canopy [64,66]. From this segmentation, tree location, height, crown diameter, and crown area can be calculated. Point cloud segmentation uses the relative spacing between trees. As spacing at the top of trees is greater than at the bottom, LiDAR points are included or excluded based on their relative distance to each other, starting at the top of the tree and working down. Using this segmentation process, the same tree attributes were created [12,64]. To extract tree metrics, a normalized point cloud consisting of vegetation points was created, and, due to the close proximity of trees in areas within the research site, the spacing threshold was set to 0.5 m.

Accuracy Assessment of OBIA and Stratified Random Sampling Results
In general, for OBIA (Table 5), the overall classification accuracy was excellent: over 90% except 1995 (89.33%). The overall Kappa statistics show strong agreement (>80%) [58] between the years, excepting 1995 (0.78%). In general, the error matrix shows that the producer's accuracy for canopy (89%-97%) is higher than the user's accuracy (84%-94%), excepting the 2017 OBIA (1928 boundary). The error matrices for the entire period 1938-2017 show that the producer's accuracy (89%-97%) (omission error) of canopy cover ranges from 3% to 11%; slightly higher than user's accuracy (84%-94%) (commission error) of 16% to 6%, excepting the 2017 OBIA (1928 boundary) (Table 5). However, different producer's and user's accuracies were obtained using the 2017 city boundary. Therefore, accuracy assessment results appear sensitive to the study site boundary. To enhance the OBIA accuracy assessment validation, ground truthing was conducted. Due to only being able to access public property combined with local flooding, the area available to ground truth within the city was approximately 1 /4 of the total 2017 boundary. Despite these restrictions, the error matrix is comparable to those within photointerpretation (Table 5).

Trend of Urban Forestry Canopy Change and the Impacts of the 1998 Tornado and Tree Diseases
For the 1928 City boundary, the OBIA analytical results (Table 7)   i-Tree's stratified random sampling confirms that from 2008 to 2019 the percentage canopy cover within the 1928 boundary is~30%-35% (Table 6). From 1938 to 1995 (57 years), canopy cover increased as a proportion of city area by 5%. From 1995 to 2017 (22 years), canopy cover increased by 10% of the city area. In fact, the 10% increase occurred between 1995 and 2008 and canopy cover remained stable from 2008 to 2019.
The canopy change map (Figure 4) shows that between 1938-1995 and 1995-2017 the majority of canopy increase took place on the east boundary of the city along the Minnesota River flood plain, within the central urban areas, e.g., streets and backyards and towards the north, west, and south boundary edges of the city. The increase in canopy area along the river flood plain is due to a change in land management; areas shifted from farmland to naturalized woodland. The tornado response is primarily responsible for the increase in the central urban reforestation. Canopy cover increased towards the edges of the city boundary, as new development of previous farmland created new areas for planting trees. While overall canopy cover increased, many areas remained non-canopy areas, or changed to non-canopy, primarily before 1995. The change to non-canopy is particularly evident along the Minnesota River flood plain between 1938-1951 and 1964-1995, where the shift in the river's position and subsequent flooding have changed the land use. Within the city center, areas have changed to non-canopy, likely due to redevelopment, e.g., dividing a large parcel containing one property and trees into smaller parcels containing fewer trees. increased as a proportion of city area by 5%. From 1995 to 2017 (22 years), canopy cover increased by 10% of the city area. In fact, the 10% increase occurred between 1995 and 2008 and canopy cover remained stable from 2008 to 2019. The canopy change map (Figure 4) shows that between 1938-1995 and 1995-2017 the majority of canopy increase took place on the east boundary of the city along the Minnesota River flood plain, within the central urban areas, e.g., streets and backyards and towards the north, west, and south boundary edges of the city. The increase in canopy area along the river flood plain is due to a change in land management; areas shifted from farmland to naturalized woodland. The tornado response is primarily responsible for the increase in the central urban reforestation. Canopy cover increased towards the edges of the city boundary, as new development of previous farmland created new areas for planting trees. While overall canopy cover increased, many areas remained non-canopy areas, or changed to non-canopy, primarily before 1995. The change to non-canopy is particularly evident along the Minnesota River flood plain between 1938-1951 and 1964-1995, where the shift in the river's position and subsequent flooding have changed the land use. Within the city center, areas have changed to non-canopy, likely due to redevelopment, e.g., dividing a large parcel containing one property and trees into smaller parcels containing fewer trees. The canopy change between 1995 and 2008 within the tornado tract ( Figure 5) shows increased canopy cover on streets and backyards, as well as along the Minnesota River flood plain. The canopy to non-canopy area change during this time span is primarily due to the tornado, as these areas have changed due to infrastructure conversion or redevelopment, e.g., houses, parking lots. The results of The canopy change between 1995 and 2008 within the tornado tract ( Figure 5) shows increased canopy cover on streets and backyards, as well as along the Minnesota River flood plain. The canopy to non-canopy area change during this time span is primarily due to the tornado, as these areas have changed due to infrastructure conversion or redevelopment, e.g., houses, parking lots. The results of this study reveal that more than 45% of the canopy cover within the defined F3 tornado (158-206mph) boundary was lost.
The butternut (Juglans cinerea L.) and, particularly, the American elm (ulmus americana L) trees were both highly visible within the urban forest and naturalized areas prior to the 1980s, but due to Dutch Elm Disease and Butternut Canker they have disappeared, excepting a few disease-resistant "survivor" trees. Due to a probable combination of image resolution and the length of time between images it was not possible to determine any correlation between the loss of these tree species and canopy cover. However, up until 1995, canopy cover remained relatively constant, therefore replanting using other tree species, e.g., maple (Acer spp.), or natural regrowth mitigated the loss.
The results primarily focused on the 1928 city boundary. However, comparing the 1928 (941 hectares) and 2017 (170 hectares) city boundary data for canopy cover for both OBIA and random sampling for the years 2008 and 2017/19, the data shows that the areas have a comparable canopy cover of~34%, even though the 2017 area is almost double that of the 1928 boundary.
Remote Sens. 2020, 12, x FOR PEER REVIEW 12 of 21 Remote Sens. 2020, 12, x; doi: FOR PEER REVIEW www.mdpi.com/journal/remotesensing this study reveal that more than 45% of the canopy cover within the defined F3 tornado (158-206mph) boundary was lost. The butternut (Juglans cinerea L.) and, particularly, the American elm (ulmus americana L) trees were both highly visible within the urban forest and naturalized areas prior to the 1980s, but due to Dutch Elm Disease and Butternut Canker they have disappeared, excepting a few disease-resistant "survivor" trees. Due to a probable combination of image resolution and the length of time between images it was not possible to determine any correlation between the loss of these tree species and canopy cover. However, up until 1995, canopy cover remained relatively constant, therefore replanting using other tree species, e.g., maple (Acer spp.), or natural regrowth mitigated the loss.
The results primarily focused on the 1928 city boundary. However, comparing the 1928 (941 hectares) and 2017 (1705 hectares) city boundary data for canopy cover for both OBIA and random sampling for the years 2008 and 2017/19, the data shows that the areas have a comparable canopy cover of ~34%, even though the 2017 area is almost double that of the 1928 boundary.

Urban Forest Canopy Density, Height Assessment and Individual Tree Metrics using LiDAR
LiDAR was used to create a canopy cover density model, CHM, and tree metrics for 2010. The canopy cover model shows the total canopy coverage within a pixel area of 10.06 m 2 . Although the resolution of canopy coverage is low compared to the OBIA 2008 data (0.28 m 2 ), the two datasets corroborate each other with respect to canopy cover ( Figure 6).

Urban Forest Canopy Density, Height Assessment and Individual Tree Metrics using LiDAR
LiDAR was used to create a canopy cover density model, CHM, and tree metrics for 2010. The canopy cover model shows the total canopy coverage within a pixel area of 10.06 m 2 . Although the resolution of canopy coverage is low compared to the OBIA 2008 data (0.28 m 2 ), the two datasets corroborate each other with respect to canopy cover ( Figure 6). Figure 7 shows canopy cover height. The majority of the highest canopy and tallest trees are in naturalized areas within the flood plain (>12.19m), which is expected due to environmental conditions, e.g., soil, nutrients, water allocation, etc., allowing extended growth of pioneer species specific to the environment, such as poplar (e.g. Populus deltoides). Isolated areas of high canopy within the city center are trees that survived the tornado. The path of the tornado through the city (Figure 7) is visible, as is the lower height of the canopy within the tornado tract compared to the canopy height north and south of the tornado boundary. The high canopy located on the west boundary in the flood plain is due to the density of trees protecting them from wind exposure as compared to isolated individual trees within the urban environment.
The tree attributes determined for both the Canopy Height Model (CHM) and point cloud (PC) segmentation processes were tree location, tree height, crown diameter, and crown area. The point cloud segmentation process also provided crown volume. However, both processes extracted a different number of trees. Table 8 shows a selection of tree attributes for trees located along the south west corner of Minnesota Square Park (Figure 8). Both the x and y tree locations are accurate. However, the remaining metrics each have accuracy discrepancies and cannot be rectified without performing regression analysis, not included within this study.
Remote Sens. 2020, 12, x FOR PEER REVIEW 13 of 21 Figure 6. LiDAR canopy density. Figure 7 shows canopy cover height. The majority of the highest canopy and tallest trees are in naturalized areas within the flood plain (>12.19m), which is expected due to environmental conditions, e.g., soil, nutrients, water allocation, etc., allowing extended growth of pioneer species specific to the environment, such as poplar (e.g. Populus deltoides). Isolated areas of high canopy within the city center are trees that survived the tornado. The path of the tornado through the city (Figure 7) is visible, as is the lower height of the canopy within the tornado tract compared to the canopy height north and south of the tornado boundary. The high canopy located on the west boundary in the flood plain is due to the density of trees protecting them from wind exposure as compared to isolated individual trees within the urban environment.  The tree attributes determined for both the Canopy Height Model (CHM) and point cloud (PC) segmentation processes were tree location, tree height, crown diameter, and crown area. The point cloud segmentation process also provided crown volume. However, both processes extracted a different number of trees. Table 8 shows a selection of tree attributes for trees located along the south west corner of Minnesota Square Park (Figure 8). Both the x and y tree locations are accurate. However, the remaining metrics each have accuracy discrepancies and cannot be rectified without performing regression analysis, not included within this study.  to the increase in canopy area along the river flood plain. Urban canopy change to non-canopy within the city center areas could also be attributed to urban redevelopment. The shift in the river's position and subsequent flooding have caused an evident change to non-canopy along the Minnesota River flood plain between 1938-1951 and 1964-1995. Nevertheless, due to the limitations of the spectral and temporal resolutions of the images used in this study, it was not possible to determine the canopy cover change caused by tree diseases.

Limitations and Potentials of Canopy Metrics Extraction using LiDAR Data Analysis
Except for tree locations, the other canopy metrics determined for both CHM and PC segmentation processes based on the LiDAR data demonstrated remarkable discrepancies (Table 8). For example, tree heights for both processes vary by 0.91-2.43 m. For crown diameter and crown area, the differences are even greater between processes, e.g., for the same tree crown diameter (CHM 25.76 m vs. point cloud 38.98 m) and crown area (CHM 521.28 m 2 vs. point cloud 1193.32 m 2 ). These canopy metrics' accuracy is limited by several factors. Firstly, the research LiDAR was collected with the leaves off; forestry research LiDAR is typically collected with the leaves on. Secondly, the LiDAR point cloud was 0.08 points/m 2 , 2.25 times less dense than the minimum commonly used within LiDAR-derived forest research, 0.18 points/m 2 [12,66,67]. Thirdly, after improvement of the LiDAR data quality, errors remained in the point cloud classification and therefore there is the potential to locate trees where none exist. More detailed information about this can be found in Blackman [68]. Fourthly, LiDAR360 offers only one algorithm for CHM segmentation, based on a paper that uses an algorithm where the research location is oak savannah woodland, not an urban forest. The only building structure was a fire lookout, removed from the data prior to the assessments [66]. Likewise, there is only one point cloud segmentation algorithm provided by LiDAR360 software based on a previous study that uses a segmentation algorithm where the research area is a mixed conifer forest, not an urban environment [12]. In conclusion, the canopy metrics analysis is limited by the quality of the publicly available LiDAR data and the limitations of the LiDAR data processing software.

Conclusions
This paper has shown the importance of remote sensing techniques to extract invaluable information from historical data sources. It demonstrated that both the OBIA and the point-based non-image classification approach could determine the temporal change in the canopy cover accurately from historical aerial images. The point-based random sampling method is an efficient method to determine the urban forest canopy cover area but lacks locational information. In contrast, the OBIA offers the ability to ascertain spatial-temporal canopy change leading to more detailed analysis, but this can be time-consuming and requires specialist knowledge. Based on this research, urban forest stakeholders can now choose whether to use OBIA or random-point-based assessment to determine urban tree canopy cover dependent on their specific purposes.
The post-classification comparison change analysis indicated that the 1998 tornado had the largest impact on canopy cover changes in the study site in addition to the other two identified factors -urban development and the shift of the Minnesota River channel. However, tree disease impact on canopy cover was undetectable in the OBIA-based classification maps derived from the historical aerial images. This might be partially attributed to a limitation of the study-only having access to hard-copy historical images acquired in some discrete years but lacking access to digital multi-spectral airborne images with high spectral and temporal resolutions.
The research also revealed the value and potential of LiDAR data and techniques to determine urban forest metrics. However, the canopy metrics analysis in this study is limited by the quality of the freely available LiDAR data and the limitations of the LiDAR data processing software. To accurately extract canopy metrics (such as crown area, diameter, and volume) in addition to locations of urban canopies, more dense LiDAR point cloud data collected during the leaf-on season, sample and independent variable data, as well as LiDAR processing software with algorithms developed specifically for urban forest analysis are needed.