Next Article in Journal
The Impact of Time Difference between Satellite Overpass and Ground Observation on Cloud Cover Performance Statistics
Previous Article in Journal
Remote Sensing of Submerged Aquatic Vegetation in a Shallow Non-Turbid River Using an Unmanned Aerial Vehicle
Previous Article in Special Issue
Application of In-Segment Multiple Sampling in Object-Based Classification
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A Versatile, Production-Oriented Approach to High-Resolution Tree-Canopy Mapping in Urban and Suburban Landscapes Using GEOBIA and Data Fusion

Spatial Analysis Laboratory, University of Vermont, Aiken Center, 81 Carrigan Dr., Burlington, VT 05405, USA
*
Author to whom correspondence should be addressed.
These authors contributed equally to this work.
Remote Sens. 2014, 6(12), 12837-12865; https://doi.org/10.3390/rs61212837
Submission received: 15 July 2014 / Revised: 8 December 2014 / Accepted: 16 December 2014 / Published: 22 December 2014
(This article belongs to the Special Issue Advances in Geographic Object-Based Image Analysis (GEOBIA))

Abstract

:
The benefits of tree canopy in urban and suburban landscapes are increasingly well known: stormwater runoff control, air-pollution mitigation, temperature regulation, carbon storage, wildlife habitat, neighborhood cohesion, and other social indicators of quality of life. However, many urban areas lack high-resolution tree canopy maps that document baseline conditions or inform tree-planting programs, limiting effective study and management. This paper describes a GEOBIA approach to tree-canopy mapping that relies on existing public investments in LiDAR, multispectral imagery, and thematic GIS layers, thus eliminating or reducing data acquisition costs. This versatile approach accommodates datasets of varying content and quality, first using LiDAR derivatives to identify aboveground features and then a combination of LiDAR and imagery to differentiate trees from buildings and other anthropogenic structures. Initial tree canopy objects are then refined through contextual analysis, morphological smoothing, and small-gap filling. Case studies from locations in the United States and Canada show how a GEOBIA approach incorporating data fusion and enterprise processing can be used for producing high-accuracy, high-resolution maps for large geographic extents. These maps are designed specifically for practical application by planning and regulatory end users who expect not only high accuracy but also high realism and visual coherence.

Graphical Abstract

1. Introduction

Tree canopy provides many essential ecological services to human populations, both at local scales and globally. These benefits include multiple environmental functions: trees intercept rain, reducing the volume of stormwater runoff and, ultimately, water pollution [1,2]; they improve air quality by filtering pollutants [3]; they moderate the heat-island effect observed in areas with large expanses of impervious surfaces [2]; and they are an important part of the global carbon cycle, serving as a carbon sink [4]. In addition to these direct environmental benefits, trees also have important effects on the social function of human communities, improving property values [5], reducing crime [6], and fostering a sense of place and safety [2,7], all factors integral to the perceived quality of life in urban and suburban communities.
Although the importance of healthy trees to environmental function and community welfare is increasingly well-known, forests worldwide are threatened by a host of anthropogenic and natural agents. Direct habitat loss from development is a primary factor of anthropogenic change in urban environments, and indeed land-cover conversions contributed to recent tree-canopy declines in multiple American cities [8]. Among natural forces, numerous insect pathogens threaten trees in many parts of the world, including the Asian long-horned beetle (ALB) [9] and emerald ash borer (EAB) [10]. For example, thousands of trees in Worcester, Massachusetts, USA were removed after discovery of ALB in 2008, entirely eliminating mature tree cover in some neighborhoods [11]. Other threats include fungal pathogens [12], ice storms [13,14], and climate-induced drought and heat stress [15]; some of these stressors could produce interactive effects that exacerbate tree decline [12]. If current projections of climate change are realized, with some parts of the world experiencing more frequent and intense storms, droughts, and heat waves, increased tree mortality could produce observable forest diebacks [16].
Fortunately, cities in the United States and Canada are responding to these threats, developing ambitious programs to maintain existing trees and to expand urban forest cover. Perhaps the most notable is New York City’s Million Trees Program, which aims to plant a million trees during a 10-year period [17], but other cities with planting programs includes Washington, DC [18], and Philadelphia, PA [19]. These proactive cities have established specific goals for enhancing tree cover, usually expressed as a desired proportion of tree canopy relative to total land area (e.g., 40% tree canopy goal for Washington, DC). Of particular emphasis for these programs is the planting of “street” trees in road medians and along road edges, which help mitigate runoff volumes from impervious surfaces. As a matter of social justice, many programs also seek to address areas where tree canopy is most lacking; these areas are often among the poorest and most underserved neighborhoods.
To be effective, tree-planting programs first need a comprehensive assessment of existing tree canopy, a baseline documentation that supports establishment of specific community goals. Traditional methods of manual photointerpretation combined with high-resolution imagery can provide informative city-wide estimates of tree cover [8], but these methods are laborious and do not highlight fine-scale patterns. Similarly, excellent field-based methods exist for assessing tree species distributions and volumes [20,21], but they are generally sample-driven and lack adequate detail for entire municipalities. Moderate-scale tree canopy maps developed from remote sensing data, such as the 30-m National Land Cover Dataset (NLCD), provide informative regional assessments of forest cover, but they too lack the needed specificity for evaluation of neighborhoods with sparse tree canopy [22]. Ideally, tree canopy maps will be accurate to the scale of individual trees, permitting analysis of municipalities at a range of scales from broad political units (e.g., city, county, or state) to individual property parcels.
In combination with high-resolution remote sensing data, geographic object-based image analysis (GEOBIA) offers a solution to this mapping challenge. It can accommodate a variety of data types, including LiDAR cloud point data, LiDAR derivatives (e.g., digital elevation models, normalized digital surface models), multispectral imagery (e.g., WorldView-2), and thematic GIS files (e.g., roads, building footprints), providing a data-fusion capability that maximizes the value of existing geospatial data investments while minimizing or eliminating the need to acquire new datasets. It also permits efficient processing of large datasets with enterprise-processing capabilities, analyzing multiple tiles simultaneously.
Recent projects have demonstrated GEOBIA’s applicability to tree-canopy mapping using rule-based expert systems in which automated workflows efficiently convert high-resolution remote sensing data into GIS-ready products [23]. Some of the most densely-developed urban areas in the United States have been mapped with this methodology, including New York City [24] and Philadelphia [25]. It has proven similarly effective with larger geographic extents encompassing a spectrum of urban, suburban, and rural land uses, including every county in the state of Maryland [26]. All of these projects have achieved classification accuracies exceeding 90%, with some as high as 98%. They have also achieved a high degree of visual realism and coherence, aesthetic factors that are very important to map end users and that often cannot be equaled with pixel-based statistical classifiers [27].
Despite these successes, GEOBIA’s relevance to production-level land-cover mapping has not been widely documented. The GEOBIA literature has traditionally focused on development, refinement, and optimal use of segmentation algorithms [27]. This emphasis is vital to continued innovation in GEOBIA techniques and broadened application to natural resources assessment and management, but it overlooks GEOBIA’s capacity for accommodating huge datasets in short timeframes, producing output that covers large geographic extents and is intended for immediate use in practical settings. GEOBIA is more than an experimental approach; it is a well-tested and adaptable set of tools that geospatial analysts can use in production environments where the emphasis is on rapid conversion of data into information.
This paper demonstrates a GEOBIA-based tree canopy mapping protocol developed by the University of Vermont Spatial Analysis Laboratory (SAL) in conjunction with the United States Department of Agriculture (USDA) Forest Service. The original impetus for this method, termed the Urban Tree Canopy (UTC) Assessment Protocol, was the growing need for high-resolution land cover maps that could accurately summarize tree canopy at the scale of individual property parcels, information that in turn could help document current conditions and prioritize tree-planting efforts in urban and suburban landscapes. It has proven to be a highly flexible method for characterizing the green infrastructure of cities, counties, and now entire states, accommodating input datasets of varying quality and quantity and providing the statistical accuracy and visual realism demanded by a diversity of stakeholders: planners, policy analysts, research scientists, elected officials, and citizens. It offers a model for producing accurate land cover data across large areas in a timely and cost-effective manner, data that will facilitate legitimate time-series analysis as tree canopy changes in response to future environmental and social conditions. To date, this protocol has been performed for more than 70 communities in North America, yielding more than 1 trillion pixels of land cover data.

2. Methodology

2.1. Data Acquisition and Processing

Limited financial resources are available to perform most UTC assessments. The SAL almost exclusively relies on existing geospatial datasets in its tree canopy mapping projects. This approach helps cities and counties demonstrate the utility of their GIS programs and remote sensing data collections, and makes many UTC assessments feasible by keeping project costs relatively low. Use of existing datasets also reduces or eliminates the time needed to acquire, process, and distribute new datasets. The input datasets for these projects are often highly variable in content and quality, reflecting the diverse motives and practices underlying data acquisition and processing. Unless datasets are systematically and irreparably flawed, however, most have some value in GEOBIA; sometimes even the approximate location of features can facilitate mapping, especially thematic GIS layers representing non-tree canopy landscape elements (e.g., road centerlines, utility transmission line rights-of-way). These layers can help refine the initial tree canopy classification by eliminating objects that are unlikely to represent actual trees.

2.1.1. Multispectral Data

Free or low-cost multispectral imagery exists for many parts of the United States. In particular, the National Agriculture Imagery Program (NAIP) provides 1-m ground sample distance (GSD), leaf-on aerial imagery for the 48 conterminous states [28]. Acquired for each state since 2002, usually at 2–3 year intervals, the standard NAIP acquisition is 4-band imagery (blue, green, red, and near-infrared). The 4-band datasets are preferred for tree canopy mapping because the near infrared band is necessary for calculating vegetation indices such as the normalized difference vegetation index (NDVI) [29]. When no higher resolution datasets are available, NAIP imagery is a useful input for discriminating vegetation from other land cover features, although locational errors can occur in areas with variable terrain [26].
Some municipalities also coordinate acquisition of high-resolution aerial orthoimagery (i.e., imagery orthorectified to accommodate relief displacement) for their planning and monitoring activities. The additional spatial detail often provided by these datasets further aids in vegetation discrimination, and the SAL has used orthoimagery with cell resolutions as high as 0.152-m GSD [24]. These data are typically acquired under leaf-off conditions to better facilitate the updating of cadastral maps. Often the local orthoimagery only contains the visible bands. In addition to these aerial imagery datasets, high-resolution satellite imagery has been available for some projects, such as WorldView-2 [30].
Most multispectral datasets require relatively little processing before use in GEOBIA automated feature extraction. NAIP imagery and other datasets are often obtained as multiple separate tiles covering the area of interest, and the usual protocol is to mosaic the tiles into a single scene. If high-resolution satellite data are used, more intensive pre-processing is required, such as pansharpening and orthorectification. These preprocessing tasks are carried out using ERDAS IMAGINE.

2.1.2. LiDAR

LiDAR point clouds in LAS format with ground point classified and a minimum of four returns are preferred for tree canopy mapping, providing maximum flexibility for examining the original LiDAR points and extracting necessary derivative products. More recent tree canopy mapping projects involve direct manipulation of LiDAR point clouds within GEOBIA software, which eliminates many of the LiDAR layer generation tasks listed below. Most projects to date, however, have relied on datasets derived from raster-based summaries of specific LiDAR components. Most notable is the normalized digital surface model (nDSM), which indicates the height of aboveground features such as trees and buildings. Various software packages are available for examining and filtering point clouds and exporting desired components to raster output, but the SAL usually uses Quick Terrain Modeler (Applied Imagery, Silver Spring, Maryland, USA) or SCALGO (Scalable Algorithmics, Aarhus, Denmark). These programs permit data extraction and a multitude of processing and output options.
Original point clouds are not always available for individual projects, and sometimes the SAL relies on LiDAR derivatives produced by other entities. This scenario eliminates processing flexibility and quality control of the final products but reduces the labor devoted to data pre-processing.

DEM

When a classified LiDAR point cloud is available, the first processing step is development of a digital elevation model (DEM) that characterizes the ground topography. In Quick Terrain Modeler or a similar program, all LiDAR returns classified as ASPRS Class 2 (Ground) [31] are filtered and exported to a raster file, in IMAGINE format (.img). The gridding options are Mean Z (LiDAR elevation) and Adaptive Triangulation with no smoothing. The output cell resolution is set according to the nominal post spacing of the LiDAR point cloud, which has been at or near 1 m for most recent projects.

DSM

A similar LiDAR-derived layer is a digital surface model (DSM) that shows the topmost elements of aboveground features, including trees and buildings. It is created by filtering all non-spurious first LiDAR returns, regardless of classification. The gridding options are Max Z and Adaptive Triangulation with a smoothing filter (Radius, 1.00 Bins; Z Tolerance, 3 m). The output cell size is again the nominal post spacing. An example 0.457-m GSD DSM for Virginia Beach, Virginia, USA is shown in Figure 1b; for comparison, a true-color 0.152-m GSD orthoimage is shown in Figure 1a.

DTM

A third surface is a digital terrain model (DTM), which helps characterize the lower sections of three-dimensional objects. In contrast to the DSM, the DTM is created by filtering last LiDAR returns from the original point cloud, regardless of classification. The gridding options are the same as those used for the DSM and the final cell resolution is based on the nominal LiDAR post spacing. This surface capitalizes on LiDAR’s penetration of partially pervious features such as trees by one or more of the pulse returns; some of the returns capture topmost leaves and branches, some penetrate to interior branches and the tree bole, and others reach the ground. An example DTM is shown in Figure 1c.

nDSM

After creating the initial LiDAR derivatives, a normalized DSM (nDSM) is created by subtracting the DEM from the DSM. This raster-processing step can be performed in many different programs; the SAL usually uses the Two Image Function option with the subtraction operator in ERDAS Imagine. The final product indicates the height of features aboveground, including the upper leaves and branches of tree canopy. An example nDSM is shown in Figure 1d.

nDTM

A normalized DTM (nDTM) is also created using the Two Image Function option, subtracting the DEM from the DSM. The resulting surface superficially resembles the nDSM, but for trees the layers capture fundamentally distinct elements in the vertical plane: the topmost features of the canopy (nDSM) versus a range of features from the canopy to mid-level branches to the tree bole (nDTM). For buildings, the nDSM and nDTM are essentially the same. An example nDTM is shown in Figure 1e.

nDSM/nDTM Difference

The difference between the nDSM and nDTM is a key modeling parameter in the preferred approach to tree canopy mapping, distinguishing highly variable, partially pervious features such as trees from impervious features with smooth surfaces, such as buildings. When the difference is high, indicating high variability in the vertical complexity of LiDAR returns, a feature is likely a tree; when the difference is low, suggesting smooth features, it is probably a building. The SAL usually calculates this layer in eCognition using Layer Arithmetics after importing the nDSM and nDTM layers, but it could be created in any raster-processing software. An example nDSM/nDTM Difference layer is shown in Figure 1f.

Intensity

LiDAR intensity values indicate the strength of individual returns and are extracted from point clouds by importing them along with all first returns, regardless of classification. The gridding options are the same as for the DSM. Intensity is generally not useful for discrimination of trees from other vegetation, but it does help distinguish impervious surfaces and thus can be helpful for removing developed features from consideration during classification.

Mean Number of Returns

The SAL’s current protocol focuses on derivation of nDSM and nDTMs, but for past projects other LiDAR derivatives have been useful for characterizing the texture of aboveground features. To create a layer showing the mean number of returns per unit area, all returns are imported regardless of classification and then summarized according to specified statistical criteria. In Quick Terrain Modeler, this operation is performed using Grid Statistics by specifying the variable as Number of Returns and the statistic as Mean. To ensure adequate statistical precision, a coarser resolution is necessary; as a rule of thumb, the output cell size should be three times the nominal post spacing.
Figure 1. Common source datasets for GEOBIA-based tree-canopy mapping. (a) A 0.152-m GSD true-color orthoimage for Virginia Beach, Virginia, USA. (b) A LiDAR-derived 0.457-m GSD digital surface model (DSM). (c) Digital terrain model (DTM). (d) Normalized digital surface model (nDSM). (e) Normalized digital terrain model (nDTM). (f) Difference between nDSM and nDTM. Trees are shown as high positive values (white); other aboveground features (e.g., buildings) have values equal to or near zero (black).
Figure 1. Common source datasets for GEOBIA-based tree-canopy mapping. (a) A 0.152-m GSD true-color orthoimage for Virginia Beach, Virginia, USA. (b) A LiDAR-derived 0.457-m GSD digital surface model (DSM). (c) Digital terrain model (DTM). (d) Normalized digital surface model (nDSM). (e) Normalized digital terrain model (nDTM). (f) Difference between nDSM and nDTM. Trees are shown as high positive values (white); other aboveground features (e.g., buildings) have values equal to or near zero (black).
Remotesensing 06 12837 g001

Z Deviation

A similar LiDAR-derived texture layer focuses on variability in height values (Z). It can be created in Quick Terrain Modeler by specifying the variable as Z and the statistic as Standard Deviation, and the output resolution should also be three times the nominal post spacing. Both Z Deviation and Mean Number of Returns show trees as highly variable surfaces while buildings and other anthropogenic features are smooth. The SAL has found that the higher-resolution nDTM is a better indicator of object texture in GEOBIA modeling, but Z Deviation and Mean Number of Returns can be useful derivatives from lower-quality LiDAR datasets containing few returns per pulse.

2.1.3. Thematic GIS Layers

Many municipalities have developed detailed vector GIS datasets representing buildings, road surfaces, road centerlines, parking lots, sidewalks, hydrology, and other elements of the natural and built environments. These datasets can substantially aid development of tree canopy maps by identifying landscape features where trees are unlikely to grow or are heavily circumscribed by anthropogenic elements. For all projects, the SAL requests all pertinent datasets from the individual municipalities or downloads them from public GIS clearinghouses. It then reviews the metadata for these layers and layers themselves, assessing whether they are sufficiently accurate and informative for discriminating trees from other features.
Occasionally, the SAL also develops its own thematic GIS datasets, usually to help distinguish easily confused landscape features. Examples include bare soil and utility rights-of-way. This work is performed in ArcGIS using heads-up digitizing techniques. In certain instances, existing vector datasets may be manually updated to reflect the conditions present in the source remotely sensed data. GEOBIA routines are often used to screen the vector data to highlight areas of possible change. The overall goal in these manipulations is not to produce stand-alone products for subsequent use by municipalities but rather to directly inform automated feature extraction. For example, the utility rights-of-way polygons that are developed are rough approximations of areas in which utility lines are present. The decision whether to create or modify thematic GIS layers is contingent on the quality of the available multispectral and LiDAR datasets: Are the remote sensing datasets sufficiently detailed to ensure high-accuracy GEOBIA output? If the answer is no or uncertain, it may be more efficient to first identify and map confounding landscape features for subsequent use in automated feature extraction.

2.2. Automated Feature Extraction

The SAL performs GEOBIA by building eCognition (Trimble Navigation Limited, Westminster, USA) rule sets that segment and classify input LiDAR and multispectral datasets, refining initial classifications with available thematic GIS datasets and context-based evaluation criteria. Rule-set development is an iterative process requiring many cycles of algorithm experimentation, testing, and refinement, proceeding to the highest possible level of accuracy and visual coherence until additional improvement is unnecessary or unprofitable (i.e., when much effort is needed to effect small, incremental refinements). As an expert system, it also requires a thorough understanding of the strengths and weaknesses of the input datasets; efficient rule sets maximize the value of high-quality layers while extracting as much information as possible from lesser-quality inputs.
For tree canopy mapping, LiDAR derivatives are generally the most informative source inputs, driving the initial segmentation steps that create preliminary objects. As the SAL’s mapping protocol has evolved, many combinations of LiDAR derivatives have been used in object creation, along with various segmentation algorithms: Multiresolution, Multi-threshold, Quadtree Based, Spectral Difference, and others. Segmentation parameters (e.g., scale, shape compactness) have similarly varied by project, requiring experimentation to identify values that produce objects of the appropriate size and shape. In all cases, the goal has been to create objects that accurately and realistically depict tree canopy to the scale of individual, isolated trees (i.e., individual trees in contiguous forest patches are not mapped separately). When no LiDAR exists or available datasets have low quality, segmentation necessarily focuses on multispectral inputs.
Classification of preliminary objects is optimally driven by a combination of LiDAR and multispectral criteria, with NDVI and Visual Brightness (i.e., sum of Red, Green, and Blue bands) among the most important spectral variables. Landscape context is also central to classification. The SAL’s UTC assessment protocol specifies development of a 7-class land-cover map: Tree Canopy, Grass/Shrubs, Bare Soil, Water, Buildings, Roads/Railroads, and Other Paved Surfaces. These classes are essential to end users of the data, to not only show the location of existing trees but also where new trees could theoretically be planted (e.g., lawns and shrubby areas, modifiable impervious surfaces such as parking lots, sidewalks, and road medians, etc.). However, these classes also help refine preliminary classifications by highlighting initial tree objects that are likely erroneous: linear objects rimming buildings; objects on elevated roadways; small, tall objects immediately adjacent to roads; etc. Accordingly, even in projects that need only the Tree Canopy class, one or more other classes are usually mapped on a preliminary basis to facilitate contextual evaluation of the initial objects. Where available, thematic GIS datasets are included in development of 7-class maps and context-based object evaluation, particularly building footprints and road surfaces.
After preliminary objects have been created, classified, and refined, a final set of processing algorithms improves representation of individual tree canopy objects and their scene-wide logical consistency. Morphological routines (Morphology and Pixel-based Object Resizing) ensure that tree canopy objects are realistically depicted, smoothing angularities and filling small gaps, and clean-up routines eliminate very small objects. These steps help produce objects that best resemble trees as observed from above, a prerequisite for high end-user confidence regardless of the actual statistical accuracy.
The data fusion capabilities of eCognition make tree canopy mapping possible; multiple inputs are needed to segment, classify, and refine objects, and to compensate for limitations in one or more datasets. They also present important logistical challenges, requiring importation and efficient processing of potentially billions of imagery pixels. Mapping for large study areas is greatly facilitated in eCognition by enterprise processing, which permits simultaneous rule-set operation on multiple tiles. Tiles are distributed across multiple processing cores using the eCognition Server.
In early projects, enterprise processing required creation of a data stack containing all LiDAR and multispectral inputs, divided into reasonably sized tiles with exactly the same dimensions prior to importation with a customized import routine. Our more recent approach, in which LiDAR point clouds serve as the primary data source, relies on the tiling structure of the LiDAR data. This process avoids extensive pre-importation tiling. Additional raster and vector layers are imported as part of the rule set at run time using the Create/Modify Project algorithm.
Running time for enterprise processing varies by tile size, computing cores, and the number of software licenses; large geographic extents (e.g., metropolitan areas and counties) may require a day or more of continuous processing. This is no small accomplishment, however, considering that the volume of pixels handled by these mapping efforts is typically in the tens to hundreds of billions. Classified maps for each tile are exported to output raster files, usually in Imagine format. The tiles are then mosaicked into a single tree-canopy map (or 7-class land cover) covering the entire study area.

2.3. Manual Corrections

No map produced by automated feature extraction will achieve perfect accuracy; tree canopy is usually too variable, especially in densely developed areas with small street trees and other plantings that are easily confused with anthropogenic features (e.g., utility poles). To gain a final measure of accuracy and to eliminate obvious, non-systematic errors, the SAL reviews all draft tree canopy maps in ArcGIS and manually draws polygons around errors of commission and omission over the raster land cover. The scale at which maps are reviewed depends on the initial quality of the map, the time and labor costs needed, and the desired level of accuracy, but often it is about 1:3000. All reviews are performed by geospatial technicians who are well versed in the fundamentals of image interpretation and land cover mapping. A separate eCognition rule set is used to incorporate manual corrections into a final land cover map. This routine segments and re-classifies tree canopy objects according to the thematic correction polygons: erroneous objects are removed and omitted objects are added.

2.4. Accuracy Assessment

For projects requiring an accuracy assessment, an adequate number of random points are selected for a study area, following standard recommendations for per-pixel assessments [32]. Although various studies have attempted assessments more appropriate for GEOBIA-produced maps, no standards or rules-of-thumb have yet been developed that can be efficiently and consistently applied to tree canopy maps. The actual land cover at each random point is determined by reviewing high-resolution reference data, usually in ArcGIS. An error matrix is produced containing the user’s, producer’s, and overall accuracies.

2.5. Retrospective Change-Detection Analysis

After production of a high-resolution tree canopy map, it is possible to gauge recent changes in forest cover by performing retrospective change analyses. A retrospective analysis assumes that the best possible GEOBIA methods and high-resolution input datasets have been used to produce a map for a near-current time period. This map is then compared to data sources from an earlier time period and modified to accommodate changes in tree canopy: gain or loss. One method involves manual comparison of a high-resolution tree canopy map to earlier multispectral imagery, drawing polygons around gained or lost tree canopy in a GIS program. The polygons are then incorporated into a new raster layer that indicates change with three thematic classes: No Change, Gain, and Loss. As with manual corrections, the SAL has found that the change-detection layer is most easily produced in eCognition by segmenting and classifying tree canopy for the latter time period using the thematic polygons. The driving factor behind this approach is the fact that the latter time period data is typically higher in quality (e.g., better horizontal accuracy and higher LiDAR point cloud density).
An alternative automated method of change detection is possible when two time periods of interest have similar high-resolution input datasets, preferably LiDAR and its derivatives. In this method, tree canopy for the latter period is mapped first using GEOBIA and then subjected to manual review. After incorporating any necessary corrections, a separate eCognition rule set is created that maps tree canopy for the earlier time period and then compares it to the corrected map for the later period, creating the three thematic classes of interest (No Change, Gain, and Loss). Last, the draft change-detection layer is reviewed relative to the high-resolution inputs for the earlier time period and corrected as necessary.

3. Results and Discussion

3.1. Program Summary

The SAL has conducted UTC assessments for more than 70 cities and counties in the U.S. and Canada, from New York City, USA to Honolulu, Hawaii, USA (Table 1), and including ecosystems as diverse as eastern temperate forests, chaparral, oak woodlands, and tropical dry forests. The input datasets for individual assessments have included billions of pixels, and some study areas have been as large as 2500 km2. The enterprise processing capabilities of eCognition made analysis of this huge volume practical, greatly expediting modeling for large cities such as New York City and Philadelphia and also for counties in Maryland, West Virginia, and elsewhere. Cumulatively, the 24 counties mapped for the state of Maryland included a land area of 25,640 km2 [26]. These projects encompassed multiple permutations of data availability, from municipalities with excellent LiDAR, multispectral imagery, and thematic GIS layers to others with multispectral data only. The available datasets also varied widely in data quality; some LiDAR datasets had holes or substantial locational offsets relative to other datasets, while certain multispectral datasets had relatively low spatial resolution, noticeable shading, and cloud cover. Similarly, thematic datasets were occasionally outdated, incomplete, or error-filled. GEOBIA’s data-fusion capabilities handled this variability, extracting usable information from individual layers and combining it in a coherent mapping protocol.
Table 1. Tree Canopy (UTC) assessments conducted for selected cities and counties in the United States and Canada, with total area of each study area (km2), area of mapped tree canopy (km2), and tree canopy as a percentage of total area (%).
Table 1. Tree Canopy (UTC) assessments conducted for selected cities and counties in the United States and Canada, with total area of each study area (km2), area of mapped tree canopy (km2), and tree canopy as a percentage of total area (%).
City/CountyState/ProvinceCountryTotal Area (km2)Tree Canopy Area (km2)% Tree Canopy
Allegheny CountyPennsylvaniaUSA1929106155
BaltimoreMarylandUSA2385724
Berkeley CountyWest VirginiaUSA51121041
BramptonOntarioCanada2672911
CambridgeMassachusettsUSA18527
Caroline CountyMarylandUSA84028934
Cuyahoga CountyOhioUSA119044437
Des MoinesIowaUSA2145626
Fairfax County (Region)VirginiaUSA114659652
HartfordConnecticutUSA471225
Honolulu (Region)HawaiiUSA66122534
Howard CountyMarylandUSA65732950
Jefferson CountyWest VirginiaUSA54920337
KitchenerOntarioCanada1373022
LancasterPennsylvaniaUSA254971428
MarkhamOntarioCanada2123818
Mecklenburg CountyNorth CarolinaUSA135666449
MississaugaOntarioCanada2924415
Montgomery CountyMarylandUSA131363048
New HavenConnecticutUSA521936
New YorkNew YorkUSA78815920
NewarkNew JerseyUSA121318215
Peel (Region)OntarioCanada124616213
PhiladelphiaPennsylvaniaUSA3697019
PittsburghPennsylvaniaUSA1516040
Prince Georges CountyMarylandUSA129164650
Richmond HillOntarioCanada1012525
RockvilleMarylandUSA351544
San JoseCaliforniaUSA3916015
SyracuseNew YorkUSA661828
TorontoOntarioCanada63017628
WashingtonDistrict of ColumbiaUSA1775531
Virginia BeachVirginiaUSA79524831
Wicomico CountyMarylandUSA103532131
Worcester (Region)MassachusettsUSA33719257

3.2. Case Studies

3.2.1. Segmentation Workflows

nDSM/nDTM with Building Footprints

The overall workflow for tree canopy mapping is shown in Figure 2. To date, the most effective segmentation steps used by the SAL have relied on a combination of the nDSM and nDTM layers, first focusing on the height of features aboveground and then on their textural differences. The key to this method is an understanding of LiDAR returns and the different information they contain: the first-return nDSM indicates the height of all features, while the last-return nDTM captures objects lower in the vertical LiDAR profile, assuming they are partially pervious. Figure 3 shows an example processing sequence for a recent project in Virginia Beach, Virginia, USA. High-quality LiDAR acquired in 2012 was available for this municipality, permitting development of 0.457-m GSD LiDAR derivatives (Figure 1). High-resolution orthoimagery (0.152-m GSD) acquired in 2013 was also available, although it contained only the visible bands (Figure 3a). The first step used a Multi-threshold Segmentation with the available nDSM to distinguish tall objects from short ones, specifying a threshold of 0.61 m (Minimum Object Size, 1) (Figure 3b). This low height threshold was necessary to capture the short, scrubby trees common in the Virginia Beach coastal environment. All tall objects, initially classed in the temporary category _Tall, were in turn segmented into meaningful sub-objects using a combination of a Quadtree Based Segmentation (Model, color; Scale, 40; weighted by multispectral bands, nDSM, nDTM, nDSM/nDTM Difference, and nDSM Slope, a slope layer created in eCognition from the nDSM) and a Multiresolution Segmentation Region Grow (Scale, 10; Shape, 0.3; Compactness, 0.8; weight values of 1 for nDSM, nDTM, nDSM/nDTM Difference, and nDSM Slope) (Figure 3c). Good thematic building footprints were available for Virginia Beach, so these were incorporated into the classification with a Chessboard Segmentation (Figure 3d). Tree canopy that overhangs buildings and other impervious surfaces plays an important role in runoff mitigation and temperature regulation, so the next step evaluated objects within thematic building footprints using nDSM/nDTM Difference; objects <0.61 m were assumed to be buildings rather than tree canopy (Figure 3e). To create a Tree Canopy class, the preliminary objects were error checked with routines evaluating building edges, roads, hydrology, utility line rights-of-way, and low LiDAR Intensity values, and their overall visual coherence and realism were improved with morphological smoothing routines (Figure 3f). The available orthoimagery did not contain a Near Infrared band, so Visible Brightness was used in lieu of NDVI whenever spectral criteria were needed in error checking and refinement.
Figure 2. Processing workflow for tree canopy mapping, from initial segmentation through manual corrections. The initial segmentation in this workflow is based primarily on LiDAR derivatives (e.g., nDSM, nDTM). The objects generated as part of the initial segmentation are then subjected to an iterative process in which they are classified, segmented, and refined using morphology algorithms. As this iterative process progresses, the classification moves from using basic imagery and LiDAR metrics to contextual information defined by spatial relationships. After automated processing, manual corrections are performed to eliminate non-systematic errors. The workflow is modified as necessary to accommodate data inputs of varying content and quality (e.g., a usable nDSM exists but no nDTM is available).
Figure 2. Processing workflow for tree canopy mapping, from initial segmentation through manual corrections. The initial segmentation in this workflow is based primarily on LiDAR derivatives (e.g., nDSM, nDTM). The objects generated as part of the initial segmentation are then subjected to an iterative process in which they are classified, segmented, and refined using morphology algorithms. As this iterative process progresses, the classification moves from using basic imagery and LiDAR metrics to contextual information defined by spatial relationships. After automated processing, manual corrections are performed to eliminate non-systematic errors. The workflow is modified as necessary to accommodate data inputs of varying content and quality (e.g., a usable nDSM exists but no nDTM is available).
Remotesensing 06 12837 g002

nDSM/nDTM without Building Footprints

The availability of thematic building footprints simplifies and expedites evaluation of overhanging tree canopy. However, these GIS datasets do not exist for every municipality, especially semi-rural or rural areas. If no building footprints are available for a study area or are poorly mapped or incomplete, additional routines focusing on nDSM/nDTM Difference are necessary, along with more extensive routines evaluating multispectral imagery and ancillary data (e.g., LiDAR Intensity). Figure 4 shows this process for Caroline County, Maryland, USA, an area for which 1-meter GSD LiDAR derivatives dating from 2005 were available. The available multispectral imagery included 1-meter GSD, 4-band NAIP acquired in 2011 (Figure 4a), an imagery dataset whose NIR band permitted calculation of NDVI. As with Virginia Beach, the initial segmentation step was a Multi-threshold Segmentation that differentiated tall from short objects; in this case, an nDSM threshold of 2 m (Minimum Object Size, 1) was necessary to adequately differentiate the study area’s taller trees from shrubs, creating the temporary category _Tall (Figure 4b). The secondary segmentation process was again a Quadtree Based Segmentation (Model, color; Scale, 40; weighted by spectral bands, nDSM, nDTM, nDSM/nDTM Difference, and nDSM Slope) followed by a Multiresolution Segmentation Region Grow (Scale, 10: Shape, 0.3; Compactness, 0.8; weight values of 1 for nDSM, nDTM, nDSM/nDTM Difference, and nDSM Slope). The lower resolutions of the input datasets produced larger objects relative to Virginia Beach, but these objects still provided meaningful divisions between trees and impervious surfaces (Figure 4c). The next step used the nDSM/nDTM Difference and spectral criteria to identify objects very likely to constitute trees; objects with a difference >2 m and NDVI > 0.1 were classified as _Candidate Tree Canopy (Figure 4d). Similarly, objects likely to be buildings were classified as _Candidate Buildings when the nDSM-nDSM Difference values were <0.2 m and NDVI < 0.1. Additional classification steps using nDSM-nDSM Difference, spectral criteria (NDVI, Visual Brightness), and object characteristics (Slope, Area, Adjacency) then mapped the remaining _Tall objects into either Tree Canopy or Buildings (Figure 4e). A set of error-checking and morphological smoothing routines produced a final tree canopy map (Figure 4f). The six-year temporal discrepancy between the LiDAR derivatives and the imagery was not ideal, producing higher rates of omission for recently-planted trees (e.g., street trees). Nonetheless, the data fusion approach maximized the value of existing datasets in capturing most of the county’s trees, and missing trees were subsequently added during manual corrections.

Other LiDAR Derivatives

An nDSM combined with other LiDAR derivatives can also be effective, including the Mean Number of Returns or Z Deviation. Like the nDSM/nDTM Difference layer, these derivatives help distinguish trees from other aboveground features by highlighting textural differences between surfaces; buildings usually have a low number of returns per unit area and the height values are precise (i.e., low Z Deviation), while trees have higher values of both. Among SAL projects, the most notable example of these statistics-based LiDAR derivatives was use of Z Deviation in high-resolution tree canopy mapping for New York City [24]. In general, however, use of these layers has been superseded by the nDSM/nDTM Difference layer; the original point cloud is necessary to produce all of these textural layers, so there is usually no need to produce the Mean Number of Returns or Z Deviation layers unless the number of returns per pulse is limited (i.e., when last returns are absent or sparse). The necessity of using a larger cell size when calculating means and standard deviations also reduces the value of these layers relative to nDSM/nDTM Difference, which usually has a cell size comparable to the point spacing of the original LiDAR point cloud.
Figure 3. Segmentation workflow for tree canopy mapping using a combination of an nDSM, nDTM, and thematic building footprints. (a) A 0.152-m GSD true-color orthoimagery for Virginia Beach, Virginia, USA. (b) Initial segmentation of tall objects (temporary class _Temp, in blue) using an nDSM. (c) Re-segmentation of tall objects. (d) Segmentation and classification of thematic building footprints (red objects). (e) Evaluation of thematic buildings for overhanging tree canopy using nDSM/nDTM Difference (overhanging tree canopy in red, uncovered building rooftops in orange). (f) Final Tree Canopy class after error checking and morphological smoothing.
Figure 3. Segmentation workflow for tree canopy mapping using a combination of an nDSM, nDTM, and thematic building footprints. (a) A 0.152-m GSD true-color orthoimagery for Virginia Beach, Virginia, USA. (b) Initial segmentation of tall objects (temporary class _Temp, in blue) using an nDSM. (c) Re-segmentation of tall objects. (d) Segmentation and classification of thematic building footprints (red objects). (e) Evaluation of thematic buildings for overhanging tree canopy using nDSM/nDTM Difference (overhanging tree canopy in red, uncovered building rooftops in orange). (f) Final Tree Canopy class after error checking and morphological smoothing.
Remotesensing 06 12837 g003
Figure 4. Segmentation workflow for tree canopy mapping using a combination of an nDSM and nDTM. (a) 1-m GSD 4-band NAIP Imagery for Caroline County, Maryland, USA. (b) Initial segmentation of tall objects (temporary class _Temp, in blue) using nDSM. (c) Re-segmentation of tall objects. (d) Identification of objects very likely to be trees based on nDSM/nDTM Difference and NDVI (light green objects). (e) Differentiation of trees from buildings using nDSM/nDTM Difference, spectral criteria, object size, slope, and adjacency (tree canopy in dark green, buildings in red). (f) Final Tree Canopy class after error checking and morphological smoothing.
Figure 4. Segmentation workflow for tree canopy mapping using a combination of an nDSM and nDTM. (a) 1-m GSD 4-band NAIP Imagery for Caroline County, Maryland, USA. (b) Initial segmentation of tall objects (temporary class _Temp, in blue) using nDSM. (c) Re-segmentation of tall objects. (d) Identification of objects very likely to be trees based on nDSM/nDTM Difference and NDVI (light green objects). (e) Differentiation of trees from buildings using nDSM/nDTM Difference, spectral criteria, object size, slope, and adjacency (tree canopy in dark green, buildings in red). (f) Final Tree Canopy class after error checking and morphological smoothing.
Remotesensing 06 12837 g004

nDSM Only

For some projects, the SAL did not have access to original LiDAR point clouds, necessitating use of derivatives created by other entities. When an nDTM or other layers that capture height-based textural differences are unavailable, spectral data, object size, and shape characteristics become more important. For example, the SAL had access only to a 0.305-m GSD nDSM for Fairfax County, Virginia, USA and adjacent municipalities, acquired in 2008, but 2-m GSD WorldView-2 multispectral imagery and an accompanying 0.5-m GSD panchromatic dataset were also available, acquired in 2011 (Figure 5a). The 8-band multispectral dataset contained two near infrared bands, permitting the use of NDVI in error checking routines evaluating initial tree canopy objects. Similar to other projects, tall features were initially identified using a Multi-threshold Segmentation with the available nDSM, based on a threshold of 6 m (Minimum Object Size, 1) (Figure 5b). Note that this threshold was 3 times the one used for Caroline County, Maryland; no single height criterion is appropriate for all study areas, and this parameter generally must be established through experimentation. Although available thematic GIS layers for building footprints were outdated for parts of the study area, they were incorporated into the initial _Tall class because at minimum they highlight the approximate location of most buildings (Figure 5c). Missing buildings isolated from other _Tall features were identified using a combination of NDVI (<0.3), Size (>112.4 m2), and Rectangular Fit (>0.75) (Figure 5d). Edges of buildings missed by the available thematic layers were similarly identified with NDVI (≤4), Size (>112.4 m2), and shape (Rectangular Fit > 0.5) criteria. Objects within all thematic or estimated buildings were then evaluated for the presence of overhanging tree canopy by performing a 2-step segmentation process: a Multiresolution Segmentation (Scale, 10; Shape, 0.2; Compactness, 0.7; weight value of 1 for the WorldView-2 panchromatic band) followed by a Spectral Difference Segmentation (Maximum Spectral Difference, 5; weight values of 1 for the Red, Green, and Blue bands and a value of 2 for the Near Infrared band). The objects created by this sequence were evaluated with multiple NDVI, Visible Brightness, Size, and Adjacency criteria (i.e., overhanging trees must be adjacent to already-identified tree canopy objects) (Figure 5e). Multiple spectral criteria were also used to eliminate erroneous draft tree canopy objects in other parts of the landscape (e.g., unmapped structures and utility poles), and the Tree Canopy class was finalized with morphological smoothing (Figure 5f).

Multispectral Imagery Only

LiDAR greatly improves the speed and accuracy of tree canopy mapping, but this data type is not always available for every project area. When necessary, the SAL has based tree canopy segmentation and classification entirely on multispectral data. A Multiresolution Segmentation weighted by the available multispectral bands produces meaningful objects, but it can be relatively slow computationally with very high resolution imagery (e.g., 0.15-m GSD imagery). In such cases, a Quadtree Based Segmentation followed by a Multiresolution Segmentation Region Grow may be a better option. Subsequent classification of draft objects usually depends on NDVI combined with texture variables such as GLCM Homogeneity (quick 8/11) or edge rasters derived using the Lee Sigma algorithm, permitting discrimination of vegetation from non-photosynthetic features and then highly-textured landscape elements (i.e., tree canopy) from more homogenous features (i.e., lawns and other grassy areas). However, the lack of height information generally introduces high rates of commission, meaning that many false tree canopy objects are captured by classification routines based exclusively on multispectral imagery. These errors ultimately require more effort in manual review and editing, underscoring the value of LiDAR to tree canopy mapping. Several Maryland counties, for example, had significant gaps in the available LiDAR derivatives, necessitating tree canopy mapping routines focused on NAIP multispectral imagery.
Figure 5. Segmentation workflow for tree canopy mapping using a combination of an nDSM and multispectral data. (a) 2-m GSD 8-band WorldView-2 imagery for the Fairfax County region, Virginia, USA. (b) Initial segmentation of tall objects (temporary class _Temp, in blue) using an nDSM. (c) Incorporation of thematic building footprints (white objects). (d) Classification of missing buildings using spectral, shape, and size criteria (yellow objects). (e) Evaluation of tree canopy overhanging buildings (draft tree canopy in light green, buildings in white, potential overhanging tree canopy in magenta). (f) Final Tree Canopy class after error checking and morphological smoothing.
Figure 5. Segmentation workflow for tree canopy mapping using a combination of an nDSM and multispectral data. (a) 2-m GSD 8-band WorldView-2 imagery for the Fairfax County region, Virginia, USA. (b) Initial segmentation of tall objects (temporary class _Temp, in blue) using an nDSM. (c) Incorporation of thematic building footprints (white objects). (d) Classification of missing buildings using spectral, shape, and size criteria (yellow objects). (e) Evaluation of tree canopy overhanging buildings (draft tree canopy in light green, buildings in white, potential overhanging tree canopy in magenta). (f) Final Tree Canopy class after error checking and morphological smoothing.
Remotesensing 06 12837 g005

3.2.2. Classification

Object Features

Initial segmentation of tree canopy objects is most effectively accomplished with LiDAR; this data type provides an unmatched capability to discriminate aboveground features from other landscape elements. As already demonstrated, however, classification usually requires a data fusion approach; a combination of LiDAR, spectral criteria, and thematic datasets provides the best method for discriminating actual tree canopy from the multitude of other tall features in urban and suburban landscapes. It is difficult, if not impossible, to highlight any single method or classification protocol; every project is different, necessitating a case-by-case evaluation of the available datasets and their ability to discern the floristic and physiognomic diversity of trees in a given study area. It is equally challenging to condense classification workflows into a brief list of illustrative steps. Some projects require long sequences of classification and re-segmentation steps that methodically work through the objects that constitute complex landscapes. Several trends are nonetheless apparent (Figure 2). NDVI remains one of the simplest and most effective spectral parameters for discriminating tree canopy objects from other features, along with Visible Brightness and the Near Infrared band. Basic geometry measurements are also informative, including Area, Length/Width, Compactness, Rectangular Fit, and Roundness. These measurements are especially helpful when error-checking draft tree canopy objects, identifying irregular or small objects that are unlikely to constitute actual trees. Texture measurements similarly help, especially when mapping trees in agricultural areas with strongly contrasting vegetation, but they can be computationally intensive. The textural indices GLCM Homogeneity (quick 8/11) and Edge Extraction Lee Sigma have proven most informative in distinguishing highly textured tree canopy from other vegetation types.
Thresholds are generally preferred with direct object classification, providing the simplest, most intuitive method for evaluating features. In some cases, however, object variability cannot be adequately addressed by thresholds that assume a linear distribution. Fuzzy classifiers with multiple variables have proven effective in these instances, although they often require substantial experimentation to identify meaningful variable combinations and classification probabilities. In several Maryland counties, for example, some building types were differentiated from tree canopy using a class definition containing Area, Mean nDSM, Mean nDSM Slope, Mean nDSM/nDTM Difference, mean nDTM, and NDVI. The buildings captured by this fuzzy classifier possessed specific roofing materials and geometries that were not captured by other classification criteria.

Contextual Refinement

Although often neglected in published GEOBIA workflows [33], landscape context is essential to refinement of initial tree canopy objects, exploiting GEOBIA’s capacity to evaluate features relative to their neighbors and other objects in a scene (Figure 2). The goal is to use an object’s spatial relationship with nearby objects to gauge the likelihood that it represents actual tree canopy. For example, it is often difficult to error check small trees using spectral data alone because these features may have sparse canopies with low NDVI values. To address this issue in a draft tree canopy map for Wicomico County, Maryland, USA, the SAL used thematic centerline data in combination with an nDSM, NDVI, and adjacency criteria to identify utility poles, light standards, and other false canopy objects along roadways. The initial map was created using 0.5-m GSD LiDAR derivatives (nDSM, nDTM, nDSM/nDTM Difference), acquired in 2011, and 1-m GSD NAIP multispectral imagery (4-band), also acquired in 2011 (Figure 6a). A raster layer showing the distance to roads, as represented by an available road centerlines dataset, was created in eCognition using the Distance Map algorithm (Figure 6b). All aboveground features were mapped with the nDSM/nDTM segmentation workflow, and building footprints were used to distinguish structures from tree canopy (Figure 6c). Draft tree canopy near roads was then identified by highlighting _Candidate Tree Canopy within a mean distance of 25 m (Figure 6d), and these features were further filtered with two sequential steps: (1) selecting objects <50 m2 in area and surrounded by Unclassified; and (2) selecting from the reduced set objects with Maximum Pixel Value >13 m and NDVI < 0.3 (Figure 6e). These criteria identified small, tall objects with relatively low NDVI and wholly encompassed by ground-level land cover classes, a combination that is unlikely for recently planted street trees but entirely possible for structural features along roads. In the final 7-class land cover map, erroneous tree canopy surrounded by short vegetation was assigned to the Grass/Shrubs class (Figure 6f).
Many additional routines are possible for error checking isolated tree canopy objects, and contextual analysis also can be applied to other landscape features when developing a comprehensive land cover map. In the SAL’s UTC assessments requiring 7-class land cover, context is often used to evaluate ground surfaces that are easily confused with multispectral data, including devegetated, compacted soils vs. concrete sidewalks (i.e., Bare Soil vs. Other Paved Surfaces) and recently turned agricultural soils vs. gravel farm roads (i.e., Grass/Shrubs vs. Roads/Railroads; agricultural fields are assigned to Grass/Shrubs in the UTC protocol). These land-cover distinctions are important to end users who need to know not only where trees are located but where additional trees could theoretically be planted, given favorable political, economic, and social circumstances. Context is where GEOBIA truly becomes an expert system; the human analyst must impose rules that translate observable landscape patterns into pertinent and usable information.

3.2.3. Manual Corrections

Automated feature extraction is performed iteratively to improve representation and accuracy of tree canopy objects, but ultimately the law of diminishing returns applies: how much time and money will be required to make small, incremental improvements? When all systematic errors have been removed from an automated classification, the SAL’s experience has shown that it is more efficient to conduct manual corrections (Figure 2). In this effort, the focus is shifted to non-systemic errors such as buildings and utility lines misclassified as trees or omitted street trees. The misclassified objects typically have non-standard attributes (e.g., a building with high NDVI and complex roof structure) that experienced image analysts can easily identify but automated algorithms miss. Interestingly, these errors have very little or no effect on statistical accuracy because automated feature extraction usually captures the overwhelming majority of trees in a given study area [25]. Qualitatively, however, the effect can be profound; non-technical end users tend to focus on realism rather than statistical accuracy, especially in areas that are well known to individual reviewers. An obvious omission or glitch in these areas will erode end-user confidence in the final product, even when the overall statistical accuracy is high.
Figure 6. Contextual error-checking of draft tree canopy along roads. (a) 1-m GSD 4-band NAIP imagery for Wicomico County, Maryland, USA. (b) Thematic road centerlines (red) were used to create a distance-to-roads raster. (c) Draft tree canopy (green) was mapped with the nDSM/nDTM segmentation workflow. (d) Tree canopy objects were identified by selecting objects within 25 m of roads (blue). (e) Erroneous tree canopy was identified by selecting objects with Area <50 m2, Maximum Pixel Value >13 m, NDVI < 0.3, and surrounded by ground-level features (Unclassified). (f) Final 7-class land cover map with erroneous tree canopy assigned to surrounding features.
Figure 6. Contextual error-checking of draft tree canopy along roads. (a) 1-m GSD 4-band NAIP imagery for Wicomico County, Maryland, USA. (b) Thematic road centerlines (red) were used to create a distance-to-roads raster. (c) Draft tree canopy (green) was mapped with the nDSM/nDTM segmentation workflow. (d) Tree canopy objects were identified by selecting objects within 25 m of roads (blue). (e) Erroneous tree canopy was identified by selecting objects with Area <50 m2, Maximum Pixel Value >13 m, NDVI < 0.3, and surrounded by ground-level features (Unclassified). (f) Final 7-class land cover map with erroneous tree canopy assigned to surrounding features.
Remotesensing 06 12837 g006
The labor devoted to manual corrections in SAL projects has varied widely, but where possible 25% or more of the total budget is reserved for manual corrections. More effort may be needed for study areas with no or low-quality LiDAR or for projects with specific accuracy requirements. For example, excellent LiDAR exists for New York City, but 65% of the total effort in its UTC assessment was allocated to manual review and editing, requiring 2025 person-hours of effort [24]. This substantial investment was justified by the need for high-resolution tree canopy maps appropriate for an epidemiological study of lung disease and air quality [34], and also by the challenge of mapping newly-planted street trees in a dense urban environment. Other projects have not had the same compelling need for extensive manual review, and still others simply did not have a budget large enough to accommodate it. Nonetheless, the value of adding a final increment of quality via manual review cannot be overstated; in a practical setting with end users who have little or no knowledge of GEOBIA, it is critical to eliminate non-systematic or obvious errors that have no bearing on total accuracy but detract from visual quality.

3.2.4. Accuracy Assessment

Accuracy assessments for selected SAL tree canopy mapping projects are shown in Table 2. In all cases, accuracies exceeded 90% for Tree Canopy and were usually higher. Because of these high accuracies, combined with the need to limit costs for many projects, accuracy assessments are not performed for every SAL tree canopy mapping project; the time and money that could be devoted to a statistical analysis are instead devoted to manual corrections that improve the overall representation of tree canopy and avoid obvious errors, as described above. This emphasis is also partly a reflection of the difficulty in properly assessing GEOBIA maps with traditional per-pixel methods; a high overall accuracy surely provides useful information, but it does not gauge the degree of overlap or visual realism. Without standard and well-tested protocols for examining objects rather than pixels, and given the documented success of GEOBIA-produced tree canopy maps, accuracy assessments become a secondary consideration in projects with tight schedules and budgets.
Table 2. Per-pixel accuracy assessments for selected Tree Canopy (UTC) assessments in the United States and Canada; user’s and producer’s accuracies refer to Tree Canopy class only.
Table 2. Per-pixel accuracy assessments for selected Tree Canopy (UTC) assessments in the United States and Canada; user’s and producer’s accuracies refer to Tree Canopy class only.
City/CountyState/ProvinceCountryUser’s AccuracyProducer’s AccuracyTotal Accuracy# Pixels
All 24 CountiesMarylandUSA99%98%-- a25,639,031,536
BaltimoreMarylandUSA97%96%94% b284,742,776
Cuyahoga CountyOhioUSA98%96%94% c12,818,033,654
New YorkNew YorkUSA98%97%98% b33,917,474,301
PhiladelphiaPennsylvaniaUSA97%96%95% b3,984,373,375
San JoseCaliforniaUSA97%92%93% d4,160,591,793
a Tree Canopy only class in final land-cover map. b Total accuracy calculated for 7-class land-cover map: Tree Canopy, Grass/Shrubs, Bare Soil, Water, Buildings, Roads/Railroads, and Other Paved Surfaces. c Total accuracy calculated for 5-class land-cover map: Tree Canopy, Soil, Water, Buildings, and Roads/Paved. d Total accuracy calculated for 8-class land-cover map: Tree Canopy, Irrigated Grass/Shrubs, Non-irrigated Grass/Shrubs, Bare Soil, Water, Buildings, Roads/Railroads, and Other Paved Surfaces.

3.2.5. Retrospective Change-Detection Analysis

Manual Review of Earlier Time Period

Multiple temporally distinct LiDAR datasets have not yet been collected for many municipalities in the United States and Canada, making automated change detection difficult. Theoretically, an automated approach could be used with two time periods captured by high-quality multispectral data, but in practice this method is too coarse; the detail provided by LiDAR is necessary to confidently map change as small as individual street trees while minimizing the volume of manual corrections. However, manual review is a feasible option for change detection in moderately sized study areas where the later time period can be mapped with automated methods.
This approach was necessary for an analysis of tree canopy loss precipitated by the Asian long-horned beetle (ALB) in Worcester, Massachusetts, USA and four adjoining towns, where thousands of trees were removed during the period 2008–2010 [11]. No LiDAR existed for either the pre- or post-removal periods, but 1-m GSD, 4-band NAIP datasets were available for both intervals. Preliminary tree canopy was mapped for the 2010 period with the Multispectral Imagery Only segmentation workflow, relying on a Multiresolution Segmentation (Scale, 25; Shape, 0.3; Compactness, 0.5; weight values of 3 for NDVI, 2 for Near Infrared and Red, and 1 for the Green band) and a sequence of classification criteria featuring textural (GLCM Homogeneity (quick 8/11)) and spectral (NDVI) thresholds. Exclusive reliance on multispectral imagery predictably resulted in overestimation of tree canopy, especially in shrubby areas with textures similar to trees, necessitating thorough manual review and editing. Tree canopy change was assessed manually by comparing the 2010 map (Figure 7a) to the 2008 NAIP (Figure 7b) and delineating losses and gains (Figure 7c). Tree loss across the entire study area was relatively small (−2%) but the spatial distribution of change was uneven; some neighborhoods suffered substantial losses, including nearly 70% of tree cover in north-central Worcester. Gains in tree canopy were negligible during the 2-year period. This study demonstrated that GEOBIA can assist retrospective change detection even when only sub-optimal input datasets are available, revealing patterns of loss that will inform subsequent damage control and remediation efforts.

Automated Comparison of Two Time Periods

Some communities have funded or received multiple LiDAR collections, making automated comparison more efficient than manual review in change detection. For example, an 8-year analysis of tree canopy change in Virginia Beach, Virginia, USA was based on 0.457-m GSD LiDAR acquired in 2012 and a 0.61-m GSD dataset acquired in 2004. As described above, the 2012 tree canopy map was produced using the nDSM/nDTM segmentation workflow aided by 0.152-m GSD true-color orthoimagery dating from 2013. Initial tree canopy objects for the 2004 map were similarly based on the nDSM/nDTM workflow, with error checking dependent on the best available multispectral imagery (1-m GSD, 4-band NAIP acquired in 2008). In areas with recently exposed bedrock or soil visible in the 2008 NAIP (i.e., very low NDVI), this temporal gap prompted error checking routines to erroneously remove some draft tree canopy objects, but otherwise the 2004 LiDAR produced a reasonable starting point for change detection. An additional eCognition routine directly compared the draft 2004 tree canopy to the reviewed and corrected 2012 layer, creating the No Change, Gain, and Loss categories. Small Gain objects (0–34.8 m2) wholly or partly enclosed by No Change objects were merged into No Change to reduce overestimation and to improve visual continuity. The draft layer was then thoroughly reviewed and modified as necessary, producing a final change detection map (Figure 8a).
An estimated increase in tree canopy area of 7.6 km2 was observed across the study area during the period 2004–2012, with a gain of 33.1 km2 offsetting a loss of 25.6 km2. Overall, the relative gain was 3.1%. Land cover changes produced the most noticeable losses, as contiguous blocks of trees were removed for residential or commercial development (Figure 8b), but established neighborhoods showed a patchwork of gains and losses; some trees were lost to attrition while new ones were planted. However, the largest apparent source of gain was the growth of existing tree canopies, particularly along the margins of forest patches and around single, isolated trees (Figure 8c). Different acquisition parameters for the two LiDAR datasets may have contributed to the magnitude of the observed increase, but the 8-year period was likely long enough to capture actual canopy growth. Change detection based on comparable LiDAR datasets was predictably more efficient for Virginia Beach, permitting analysis of a land area 2.4 times that of Worcester. Extensive quality control was still needed to fix errors in the change detection classes attributable to the 2004 tree canopy map, but subsequent change analyses will require much less effort if later time periods are mapped with similar GEOBIA techniques and LiDAR inputs.
Figure 7. Retrospective change detection analysis for Worcester, Massachusetts, USA based on manual comparison of GEOBIA-produced high-resolution tree canopy map to earlier multispectral imagery. (a) 2010 Tree Canopy (green) superimposed on 1-m GSD 4-band NAIP imagery from which it was modeled. (b) Manual comparison of 2010 Tree Canopy to 1-m GSD 4-band NAIP imagery acquired in 2008, with drawn lines representing tree canopy lost during 2-year period (red) and tree canopy gained (yellow). (c) Final change-detection layer superimposed on 2010 NAIP, with No Change (green), Loss (red), and Gain (yellow) classes.
Figure 7. Retrospective change detection analysis for Worcester, Massachusetts, USA based on manual comparison of GEOBIA-produced high-resolution tree canopy map to earlier multispectral imagery. (a) 2010 Tree Canopy (green) superimposed on 1-m GSD 4-band NAIP imagery from which it was modeled. (b) Manual comparison of 2010 Tree Canopy to 1-m GSD 4-band NAIP imagery acquired in 2008, with drawn lines representing tree canopy lost during 2-year period (red) and tree canopy gained (yellow). (c) Final change-detection layer superimposed on 2010 NAIP, with No Change (green), Loss (red), and Gain (yellow) classes.
Remotesensing 06 12837 g007
Figure 8. Retrospective change detection analysis for Virginia Beach, Virginia, USA based on automated comparison of GEOBIA-produced high-resolution tree canopy maps for two time periods, 2004 and 2012. (a) 2004–2012 change-detection layer with No Change (green), Loss (red), and Gain (yellow), superimposed on 0.152-m GSD 3-band orthoimagery. (b) Loss attributable to land cover conversion from forest to developed land uses. (c) Gain along margins of tree canopy.
Figure 8. Retrospective change detection analysis for Virginia Beach, Virginia, USA based on automated comparison of GEOBIA-produced high-resolution tree canopy maps for two time periods, 2004 and 2012. (a) 2004–2012 change-detection layer with No Change (green), Loss (red), and Gain (yellow), superimposed on 0.152-m GSD 3-band orthoimagery. (b) Loss attributable to land cover conversion from forest to developed land uses. (c) Gain along margins of tree canopy.
Remotesensing 06 12837 g008

4. Conclusions

It is possible to move GEOBIA from a purely academic setting to a production environment, shifting focus from parameter experimentation to timely conversion of huge volumes of data into usable information. Administrators and planners at all levels of government need current and accurate tree canopy maps, but for many parts of the United States and Canada the available data are too coarse or limited in extent. GEOBIA techniques offer solutions to this information gap, facilitating simultaneous analysis of multiple pertinent inputs, including many high-resolution datasets already acquired at substantial cost but underutilized or overlooked. Many possible modeling approaches exist, and indeed the specific segmentation and classification procedures for individual projects must be tailored to the available datasets, time, funds, and processing capabilities; no single solution is optimal for all cases.
LiDAR derivatives such as Normalized Digital Surface Models (nDSMs) and Normalized Digital Terrain Models (nDTMs) have consistently proven to be the most important remote sensing datasets for mapping tree canopy, permitting efficient discrimination of tall, highly textured objects from other aboveground features. However, GEOBIA can accommodate many different types of data and varying levels of quality, providing a robust method for differentiating trees from surrounding landscape features. Context is also key to tree canopy mapping; procedures that identify contrast with other features best approximate the way humans perceive variability among landscape features. Combined with GEOBIA’s versatile data fusion capabilities, this contextual, expert system based approach provides both a quantitative accuracy and high visual realism that per-pixel, probabilistic methods cannot match. As more communities acquire second rounds of temporally distinct LiDAR, it will be possible to perform more retrospective change detection analyses that help document the effects of anthropogenic and natural environmental change.
Production environment GEOBIA is feasible for many other landscape mapping efforts covering large geographic extents, including comprehensive land use/land cover (LULC) mapping, vegetation classification, wetlands delineation, agricultural monitoring, forest characterization, and river channel modeling. As always, overall quantitative accuracy will depend on the quality and resolution of the input imagery, LiDAR, and thematic datasets, along with a fundamental understanding of the features of interest and their relationship to other landscape objects. Where output products are expected to be highly scrutinized, manual editing helps remove non-systematic inconsistencies that affect end-user confidence, even if its effect on quantitative accuracy is negligible. GEOBIA thus integrates the old with the new, joining sophisticated automated feature extraction with time-tested elements of image interpretation into a single, efficient expert system.

Acknowledgments

Funding for the work presented in this article has come from collaborative research with the USDA Forest Service’s Northern Research Station, the National Science Foundation (awards 1053559, 1065785, and DEB 1027188), the National Aeronautics and Space Administration (grant NNX12AN07G), and from the numerous communities with which we have partnered. We would like to acknowledge the contributions made by our collaborators, particularly J. Morgan Grove from the USDA Forest Service, Austin Troy from the University of Colorado, Dexter Locke from Clarke University, Michael Galvin from the SavATree Consulting Group, and Ralph Dubayah from the University of Maryland.

Author Contributions

Jarlath O’Neil-Dunne conceived the overall approach to urban tree canopy assessments and developed the initial eCognition rule sets for mapping tree canopy from LiDAR derivatives and multispectral data. He has also supervised many of the UTC assessments for individual cities and municipalities. Sean MacFaden has similarly developed rule sets for many individual UTC projects and guided delivery of maps and summary data to end users. Anna Royar coordinated review and manual correction of the draft tree canopy maps for many projects. Jarlath O’Neil-Dunne and Sean MacFaden drafted and refined the manuscript, while Anna Royar provided review and editorial assistance.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Brack, C.L. Pollution mitigation and carbon sequestration by an urban forest. Environ. Pollut. 2002, 116, S195–S200. [Google Scholar] [CrossRef] [PubMed]
  2. Nowak, D.J.; Dwyer, J.F. Understanding the benefits and costs of urban forest ecosystems. In Urban and Community Forestry in the Northeast, 2nd ed.; Kuser, J.E., Ed.; Springer Netherlands: Heidelberg, Germany, 2007; pp. 25–46. [Google Scholar]
  3. Nowak, D.J.; Crane, D.E.; Stevens, J.C. Air pollution removal by urban trees and shrubs in the United States. Urban. For. Urban. Green. 2006, 4, 115–123. [Google Scholar] [CrossRef]
  4. Nowak, D.J.; Crane, D.E. Carbon storage and sequestration by urban trees in the USA. Environ. Pollut. 2002, 116, 381–389. [Google Scholar] [CrossRef] [PubMed]
  5. Mansfield, C.; Pattanayak, S.K.; McDow, W.; McDonald, R.; Halpin, P. Shades of green: Measuring the value of urban forests in the housing market. J. For. Econ. 2005, 11, 177–199. [Google Scholar]
  6. Troy, A.; Grove, J.M.; O’Neil-Dunne, J. The relationship between tree canopy and crime rates across an urban-rural gradient in the greater Baltimore region. Landsc. Urban. Plan. 2012, 106, 262–270. [Google Scholar] [CrossRef]
  7. Kuo, F.E.; Bacaicoa, M.; Sullivan, W.C. Transforming inner-city landscapes: Trees, sense of safety, and preference. Environ. Behav. 1998, 30, 28–59. [Google Scholar] [CrossRef]
  8. Nowak, D.J.; Greenfield, E.J. Tree and impervious cover change in U.S. cities. Urban. For. Urban. Green. 2012, 11, 21–30. [Google Scholar] [CrossRef]
  9. Haack, R.A.; Hérard, F.; Sun, J.; Turgeon, J.J. Managing invasive populations of Asian longhorned beetle and citrus longhorned beetle: A worldwide perspective. Annu. Rev. Entomol. 2010, 55, 521–546. [Google Scholar] [CrossRef] [PubMed]
  10. Kovacs, K.F.; Haight, R.G.; McCullough, D.G.; Mercader, R.J.; Siegert, N.W.; Liebhold, A.M. Cost of potential emerald ash borer damage in U.S. communities, 2009–2019. Ecol. Econ. 2010, 69, 569–578. [Google Scholar] [CrossRef]
  11. Dodds, K.J.; Orwig, D.A. An invasive urban forest pest invades natural environments—Asian longhorned beetle in northeastern US hardwood forests. Can. J. For. Res. 2011, 41, 1729–1742. [Google Scholar] [CrossRef]
  12. Desprez-Loustau, M.; Marçais, B.; Nageleisen, L.; Piou, D.; Vannini, A. Interactive effects of droughts and pathogens in forest trees. Ann. For. Sci. 2006, 63, 597–612. [Google Scholar] [CrossRef]
  13. Irland, L.C. Ice storms and forest impacts. Sci. Total Environ. 2000, 262, 231–242. [Google Scholar] [CrossRef] [PubMed]
  14. Sisinni, S.M.; Zipperer, W.C.; Pleninger, A.G. Impacts from a major ice storm: Street-tree damage in Rochester, New York. J. Arboric. 1995, 21, 156–167. [Google Scholar]
  15. Allen, C.D.; Macalady, A.K.; Chenchouni, H.; Bachelet, D.; McDowell, N.; Macalady, M.; Kitzberger, T.; Rigling, A.; Breshears, D.D.; Hogg, E.H.; et al. A global overview of drought and heat-induced tree mortality reveals emerging climate change risks for forests. For. Ecol. Manag. 2010, 259, 660–684. [Google Scholar] [CrossRef]
  16. Intergovernmental Panel on Climate Change (IPCC). Climate Change 2014: Impacts, Adaptation, and Vulnerability. Part. A: Global and Sectoral. Aspects. Contribution of Working Group II to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change; Field, C.B., Barros, V.R., Dokken, D.J., Mach, K.J., Mastrandrea, M.D., Bilir, T.E., Chatterjee, M., Ebi, K.L, Estrada, Y.O., Gonova, R.C., Eds.; Cambridge University Press: Cambridge, UK and New York, NY, USA, 2014; Chapter 4; pp. 33–40. [Google Scholar]
  17. MillionTreesNYC. A PlaNYC Initiative with NYC Parks and New York Restoration Project. New York City Parks & Recreation and New York Restoration Project (2007). Available online: http://www.milliontreesnyc.org/html/home/home.shtml (accessed on 20 June 2014).
  18. DC Urban Forestry Administration. District of Columbia Assessment of Urban Forest Resources and Strategy. Urban Forestry Administration, District Department of Transportation, Government of the District of Columbia (2010). Available online: http://ddot.dc.gov/sites/default/files/dc/sites/ddot/publication/attachments/dc_assessment_urban_forest_resources_strategy_2010–06.pdf (accessed on 20 June 2014).
  19. City of Philadelphia. Greenworks Philadelphia. The City of Philadelphia, Mayor’s Office of Sustainability (2009). Available online: http://www.phila.gov/green/greenworks/pdf/Greenworks_OnlinePDF_FINAL.pdf (accessed on 20 June 2014).
  20. Nowak, D.J.; Crane, D.E.; Stevens, J.C.; Hoehn, R.E.; Walton, J.T.; Bond, J. A ground-based method of assessing urban forest structure and ecosystem services. Arboric. Urban. For. 2008, 34, 347–358. [Google Scholar]
  21. i-Tree—Tools for Assessing and Managing Community Forests. USDA Forest Service. Available online: http://www.itreetools.org (accessed on 20 June 2014).
  22. Walton, J.T.; Nowak, D.J.; Greenfield, E.J. Assessing urban forest canopy cover using airborne or satellite imagery. Arboric. Urban. For. 2008, 34, 334–340. [Google Scholar]
  23. Benz, U.C.; Hofmann, P.; Willhauck, G.; Lingenfelder, I.; Heynen, M. Multi-resolution, object-oriented fuzzy analysis of remote sensing data for GIS-ready information. ISPRS J. Photogramm. Remote Sens. 2004, 58, 239–258. [Google Scholar] [CrossRef]
  24. MacFaden, S.W.; O’Neil-Dunne, J.P.M.; Royar, A.R.; Lu, J.W.T.; Rundle, A.G. High-resolution tree canopy mapping for New York City using LiDAR and object-based image analysis. J. Appl. Remote Sens. 2012, 6. [Google Scholar] [CrossRef]
  25. O’Neil-Dunne, J.P.M.; MacFaden, S.W.; Royar, A.R.; Pelletier, K.C. An object-based system for LiDAR data fusion and feature extraction. Geocarto. Int. 2013, 28, 227–242. [Google Scholar] [CrossRef]
  26. O’Neil-Dunne, J.; MacFaden, S.; Royar, A.; Reis, M.; Dubayah, R.; Swatantran, A. An object-based approach to statewide land cover mapping. In Proceedings of ASPRS 2014 Annual Conference, Louisville, KY, USA, 23–28 March 2014.
  27. Blaschke, T. Object based image analysis for remote sensing. ISPRS J. Photogramm. Remote Sens. 2010, 65, 2–16. [Google Scholar] [CrossRef]
  28. USDA Farm Service Agency, Aerial Photography Field Office. National Agriculture Imagery Program (NAIP). Available online: http://www.fsa.usda.gov/FSA/apfoapp?area=home&subject=prog&topic=nai (accessed on 24 June 2014).
  29. Lillesand, T.M.; Kiefer, R.W.; Chipman, J.W. Remote Sensing and Image Interpretation, 6th ed.; John Wiley & Sons: Hoboken, NJ, USA, 2008; pp. 464–466. [Google Scholar]
  30. Digital Globe, Inc. Digital Globe Content Collection: WorldView-2 Details. Available online: http://digitalglobe.com/about-us/content-collection#worldview-2 (accessed on 24 June 2014).
  31. American Society for Photogrammetry & Remote Sensing (ASPRS). LAS Specification Version 1.4-R11. Available online: http://info.asprs.org/society/committees/standards/LAS_1_4_r11.pdf (accessed on 19 December 2014).
  32. Congalton, R.G. A review of assessing the accuracy of classifications of remotely sensed data. Remote Sens. Environ. 1991, 37, 35–46. [Google Scholar] [CrossRef]
  33. O’Neil-Dunne, J.P.M.; MacFaden, S.W.; Pelletier, K.C. Incorporating contextual information into object-based image analysis workflows. In Proceedings of ASPRS 2011 Annual Conference, Milwaukee, WI, USA, 1–5 May 2011.
  34. Lovasi, G.S.; O’Neil-Dunne, J.P.M.; Lu, J.W.T.; Sheehan, D.; Perzanowski, M.S.; MacFaden, S.W.; King, K.L.; Matte, T.; Miller, R.L.; Hoepner, L.A.; et al. Urban tree canopy and asthma, wheeze, rhinitis, and allergic sensitization to tree pollen in a New York City birth cohort. Environ. Health Perspect. 2013, 121, 494–500. [Google Scholar] [PubMed]

Share and Cite

MDPI and ACS Style

O'Neil-Dunne, J.; MacFaden, S.; Royar, A. A Versatile, Production-Oriented Approach to High-Resolution Tree-Canopy Mapping in Urban and Suburban Landscapes Using GEOBIA and Data Fusion. Remote Sens. 2014, 6, 12837-12865. https://doi.org/10.3390/rs61212837

AMA Style

O'Neil-Dunne J, MacFaden S, Royar A. A Versatile, Production-Oriented Approach to High-Resolution Tree-Canopy Mapping in Urban and Suburban Landscapes Using GEOBIA and Data Fusion. Remote Sensing. 2014; 6(12):12837-12865. https://doi.org/10.3390/rs61212837

Chicago/Turabian Style

O'Neil-Dunne, Jarlath, Sean MacFaden, and Anna Royar. 2014. "A Versatile, Production-Oriented Approach to High-Resolution Tree-Canopy Mapping in Urban and Suburban Landscapes Using GEOBIA and Data Fusion" Remote Sensing 6, no. 12: 12837-12865. https://doi.org/10.3390/rs61212837

APA Style

O'Neil-Dunne, J., MacFaden, S., & Royar, A. (2014). A Versatile, Production-Oriented Approach to High-Resolution Tree-Canopy Mapping in Urban and Suburban Landscapes Using GEOBIA and Data Fusion. Remote Sensing, 6(12), 12837-12865. https://doi.org/10.3390/rs61212837

Article Metrics

Back to TopTop