A Comprehensive Approach for Mapping Ecological System Types Using Sentinel‐2 Satellite Imagery

1. Introduction

Spatial depictions of land-cover are integral to effective ecological and environmental management. Additionally, the thematic resolution of land-cover and vegetation maps contributes significantly to the effective development of policies and land management decisions based on spatial data. The importance of classifying and mapping vegetation has been widely acknowledged given the essential role it plays in wildlife habitat, biodiversity, ecosystem services and climate regulation [1,2]. Other considerations, such as the implementation of ecosystem restoration or preservation programs, also require that vegetation be mapped throughout time to serve as reference points for future comparisons.

Spatial ecosystems classifications can be broadly lumped into one of two categories: 1) geophysical-based classifications derived mainly from geophysical settings and abiotic factors, and 2) integrated classifications that combine remotely sensed land cover data with abiotic and geophysical factors. While ecological mapping efforts based on stratifications of various geophysical data layers provide a good foundation for ecological assessments, they are conceptual in nature and depict reference states rather than the actual location of different vegetation or ecological types. Rather, they serve more as a guide to potential ecological states in various locations [3,4].

Some past efforts have been made to link remote sensing and ancillary modeling approaches to produce detailed ecosystem and vegetation maps. Jung et al. (2020) used global 100-meter resolution land-cover along with various rulesets to generate a 47-class global map of terrestrial habitat types [5]. In another study, a 250-meter resolution global dataset depicting 431 world terrestrial ecosystems was created using a combination of climate regions, landform models and land-cover, and biogeographic realm [6]. However, while such global datasets serve as useful tools for national or multi-state planning efforts, they do not offer the precision that is often required for practitioners at county, state, or provincial scales to meet more localized objectives, either in terms of spatial or thematic resolution [7]. There have been some notable attempts to incorporate remotely sensed data to produce improved resolution ecosystem and vegetation maps. For instance, the Gap Analysis Project (GAP) and LANDFIRE Existing Vegetation Type (EVT) datasets provide terrestrial ecological systems classifications for the entire United States, driven by classifications of 30-meter Landsat satellite imagery combined with bio- and geophysical setting data [8,9]. Yet these efforts are conduct nationally, and their accuracy is generally low when compared with field-based data and other land-cover datasets [10,11,12]. In addition, their ability to realistically circumscribe vegetative features using a 900 m² grid cell resolution is limited [13,14]. However, the detection and realistic spatial depiction of landscape features, particularly vegetation, is of great importance for a wide range of environmental and conservation management concerns. Specifically, the spatial and thematic resolution of vegetation data has implications for water estimations in hydrological models, fire intensity and extent in predictive wildfire models, sustainability planning, and wildlife habitat management [15,16,17]. Furthermore, the satellite spatial resolution underlying a given ecological systems classification is critical since large amounts of tree and shrub cover can go completely undetected at coarser scales [18]. This is of particular importance in landscapes such as mixed woodland-grassland ecosystems, such as savannas, or transitional zones where forest and herbaceous types intermix [19,20].

Given that widely available land-cover and vegetation datasets are generally lacking in either thematic resolution or spatial resolution, there is ongoing need to integrate modeling techniques and ancillary data with traditional satellite image classifications to produce more detailed ecosystem maps over large areas. Results from such mapping efforts can serve as effective decision tools and model inputs for local and regional conservation planners and environmental modelers. To address this, we present a comprehensive approach using satellite imagery, field training data, machine learning, image segmentation, and ancillary modeling approaches to develop high-resolution ecological system type maps statewide for Arkansas, USA (Figure 1). The overarching goal of this work was to address the shortcomings described for other ecological and vegetation mapping efforts, and to provide a framework for future mapping efforts to facilitate improved natural resource management, wildlife conservation, and environmental modeling.

2. Materials and Methods

2.1. Study Area

The study area included the entire state of Arkansas, which is situated in the Southeastern United States. The state resides at the confluence of several significant ecoregions, including the Ozark Highlands, Arkansas Valley, South Central Plains, Ouachita and Boston Mountains, Mississippi Alluvial Plain, and Mississippi Valley Loess Plains [21]. These regions represent diverse land-cover and ecosystem types, ranging from rugged mountainous areas in the central-western and northwestern portions of the state to low-lying areas along the Mississippi River in the east.

We identified the prevailing ecological systems for the entire state (Figure 2) using a combination of two potential natural vegetation datasets. The first of these was developed by grouping SSURGO [22] map unit polygons to create soil eco-groups based on shared characteristics and classifying the grouped map units based on NatureServe Terrestrial Ecological Classifications [23]. The grouping process was facilitated by the NRCS ecological site descriptions, which are attached to SSURGO soil components [24,25]. We assigned soil eco-group using expert opinion where the NRCS had not yet assigned ecological sites to soil components in some eastern portions of the state. The second dataset used covered the Mississippi Alluvial Valley region, and depicted potential natural vegetation based on basin-level functional wetland assessments, which were drafted from a merging of surface geology, soils maps, and flooding maps [26,27].

2.2. Training Data Generation & Classification of Sentinel-2 Imagery

Training data for the landcover classification were generated primarily from photointerpretation of high-resolution aerial imagery using a random grid of points within each of three mapping zones covering the state (Figure 3). A total of 12 land cover class targets were identified (Table 1). Field collected plots were not directly used as samples for the landcover classification but were used to aid in photo survey sample interpretation, particularly for distinguishing between eastern redcedar versus pine landcover classes. The total number of samples collected within each of the primary mapping zones ranged from 2500 to 3000.

Remotely sensed data from the Copernicus Sentinel-2 earth observation mission [28] was surveyed in order to identify cloud-free imagery for the entire state of Arkansas. The Sentinel-2 mission consists of a constellation of two satellites (Sentinel-2A and Sentinel-2B), each carrying a Multispectral Instrument (MSI) designed to measure the reflected radiance of Earth in 13 spectral bands. Of interest for this project were: visible and near-infrared bands at 10 meters, and red edge and SWIR bands at 20 meters [29]. Imagery for the study region was selected based on a combination of recency and clarity (i.e., the absence of clouds). Given the extent of the study area, it was necessary to assess and compile imagery for multiple zones (Zone 1, Zone 2, and Zone 3), each corresponding to Sentinel-2 orbits (Figure 3). Multiple imagery dates from 2020—2021 (early growing season, late growing season, leaf-off) were selected for each of the three primary Sentinel-2A/2B orbits/zones passing over the state of Arkansas. After selecting the imagery, Sentinel-2 Level-2A (atmospherically corrected surface reflectance) image tiles were loaded into the Google Earth Engine (GEE) cloud computing platform for analysis. A subset of the available spectral bands for each image date was used. The selected bands included: band 2 (blue), band 3 (green), band 4 (red), bands 5, 6, and 7 (red-edge 1, 2, and 3), band 8 (near-infrared 1), band 8a (near-infrared 2), band 11 (shortwave infrared 1), and band 12 (shortwave infrared 2). The spatial resolution of bands 2 through 8 was 10 meters. Bands 5, 6, 7, 8A, 11 and 12, whose native resolutions are 20 meters, were resampled to 10 meters. In addition to the spectral bands, the Normalized Difference Vegetation Index (NDVI) [30], Enhanced Vegetation Index-2 (EVI2) [31], and the Normalized Difference Water Index (NDWI) [31], were calculated. Since it was designed as a complementary index to capture changes in canopy water content, the NDWI was calculated only for the late leaf-on dates, as they reflected peak vegetation/canopy growth.

The photo interpreted training samples were then used to train a RandomForest [32] classifier. Parameter values for all model runs were set as follows: number of trees = 1500, and number of variables randomly sampled at each split = 6. In order to optimize the landcover classification results, multiple model runs were executed for each zone based on iteratively modified training data. The final RandomForest model for each of the zones was then applied to their respective image stack in order to generate a prediction raster, one for each zone. Finally, additional samples for pine woodland and forest and eastern redcedar woodland and shrubland were obtained for zones 1 and 2. These samples were later used to modify the landcover assignment for all pixels that were initially classified as pine or eastern redcedar in earlier iterations of the classification. The final raster zones were merged to produce a single statewide land-cover dataset to be used in the subsequent modeling process.

2.3. Image Classification Refinements

Several raster-based mask layers were developed in order to further refine the RandomForest landcover predictions. These included: (1) an urban area mask, outside of which pixels classified as urban landcover were recoded to barren, (2) a crop area mask, outside of which pixels classified as crop were coded to herbaceous vegetation, (3) a water mask, within which pixels classified as urban and other types of land cover were converted to water to alleviate the issue of errant pixels situated within water bodies, and (4) a pine landscape mask, which was used to model pine plantation types.

2.3.1. Urban Masks

An urban area correction mask was created from the Microsoft computer generated building footprints for the United States database [33], as well as the National Land Cover Database (NLCD) [34] urban landcover classes. Building footprints converted into a 10-meter raster layer and were expanded by five pixels. The NLCD urban classes (classes 21 through 24) were re-sampled from 30-meter to 10-meter pixels and also expanded by five pixels (i.e., 50 meters). This mask was used to reclassify all pixels outside the mask that were originally classified as Developed Impervious Cover or Developed Mixed Intensity Urban to Barren in order to minimize the amount of spurious urban pixels across the landscape.

2.3.2. Crop and Water Masks

A crop mask was produced from the annual Cropland Data Layer (CDL) [35], where a given pixel was included if it was classified as crop in at least one of the five years from 2017 through 2021. Row Crop landcover outside of the mask was reclassed to Herbaceous landcover to minimize the amount of errantly classified Row Crop cover. A water mask was produced using the Copernicus Programme Global Surface Water dataset (GSW) [36]. Multiple 30-meter surface water occurrence tiles corresponding to the statewide extent of Arkansas were mosaicked and resampled to 10-m pixels. Areas where water occurred for at least 25 years during the period from 1984 to 2021 were then selected, and the result was shrunk by two pixels. This mask was then used to reclassify errant landcover pixels residing within areas of water.

2.4. Modeling Ecological Systems

Modeling of the ecological mapping targets was a complicated, multiple-step process. This process is outlined in detail below.

2.4.1. Generation of Image Objects

Trimble eCognition software was used to generate image objects. Two dates of Sentinel-2 data were used for image object generation: late season leaf-on and leaf-off. Data stacks included the red, green, blue, and near infrared (RGBN) bands. The red and blue bands were assigned a weight of 1 and the near infrared and green bands were weighted 2. The scale parameter was 40, shape 0.1, and compactness 0.7. Objects were developed by 100 km square tile. This process produced a vector feature class containing polygons that represented relatively homogeneous units (relative to the Sentinel-2 datasets). The average size of the generated polygons was approximately 0.3 hectares. Since this represents the same area as 30, 10 m Sentinel-2 pixels, any ‘salt and pepper’ effects produced during landcover classification from landcover assigned to 10-meter pixels was substantially smoothed.

2.4.2. Attribution of Landcover and Geophysical Setting to Image Objects

Landcover class and geophysical setting information were assigned to each image object. Percent slope and solar radiation were calculated from a 10 m digital elevation model (DEM) that was derived by resampling a 1 m resolution DEM (https://gis.arkansas.gov/product/dem-1m-2018/). The means of continuous variables were calculated for each object using eCognition and these variables included percent slope, solar radiation, and soil moisture. A unique identifier for the soil polygons within the SSURGO dataset was also attributed to each object based on the majority soil polygon occupying the object. Soil eco-group was defined by grouping Ecological Sites (referred to as ecoclassid and ecoclassname within the SSURGO database) that were assigned to soil map units within the NRCS gSSURGO database [37]. We assigned soil eco-group using expert opinion where the NRCS had not yet assigned Ecological Sites to soil polygons, mainly in the eastern Ozarks and Ouachitas. Heads-up inspection and re-grouping of soil map units was accomplished based on field plot data information, as well as other information extracted from gSSURGO data including percent slope, soil texture, flooding, impervious layer, and parent material.

2.4.3. Ecological Type Modeling and Mapping (Master Model)

Modeling and mapping for most ecological types was accomplished by assigning each combination of landcover and prevailing ecological system to a unique ecological type (see Figure 1). Hence, a list of eco-groups or potential natural vegetation types were in one column and landcover class in other columns, and each cell in the table was assigned an ecological type. Python scripts were used to implement the modeling rules. This provided a flexible, repeatable, and transparent method for model implementation.

2.4.4. Ecological Type Modeling and Mapping (Special Cases)

A number of ecological types were modeled using methods not outlined in the master model. These types are described below.

Developed types were classified using a mask based on the NLCD 2019 Developed Impervious Descriptor dataset [38] to identify vegetation likely to be influenced by urbanization impacts. This NLCD raster product was resampled to 10-meter resolution. Resampled pixels belonging to NLCD classes 20 (primary road) or 29 (energy production) were expanded by one pixel and pixels classed as 24 (non-road/non-energy impervious) or 25 (Microsoft buildings) were expanded by 3 pixels. All mask patches < 10 hectares were then eliminated. Within the mask, pixels were reclassified to either Developed Herbaceous Vegetation or Developed Wooded Vegetation depending on their respective types in the original landcover classification.

Floodplain ecological type models were applied to objects within the Mississippi Alluvial Plain, Red River floodplain, South-Central Interior, and West Gulf Coastal Plain regions separately. These regions were defined by selection of HUC-12 watersheds to create watershed-enforced boundaries of EPA Level 4 Ecoregions [39,40]. The floodplain modeled types were extracted from the master model and applied to the objects based on their basin location. Flooding frequency modifications were then applied to the high- and low-floodplain and high- and low-flatwoods types in the Mississippi Alluvial Valley. This was accomplished using inundation frequency from a dataset developed from LANDSAT satellite imagery for the years 1983—2011 [41]. Mean percentage inundation frequency was assigned to the objects. High floodplain types and typic or high flatwoods types were modeled in areas that had inundation frequency <= 20% during the period analyzed, and low floodplain and low flatwoods types were modeled in areas that had an inundation frequency > 20%.

Riparian types were classified based on polyline features depicting named streams in the NHDPlus database [42]. Image objects within 20 meters of a named stream were modeled as riparian types and the same watershed-enforced ecoregion boundaries used for the floodplain type modeling were used to inform riparian ecological type assignment.

Pine plantation types were mapped only within a pine landscape mask. This was done by first identifying areas of vegetation loss throughout the state. Vegetation loss information was obtained from the USDA Forest Service Landscape Change Monitoring System database (LCMS) [43] which identifies year-to-year vegetation change. A set of 500 random points were generated within the year-over-year LCMS vegetation loss footprint from 2011 to 2021. These were classified as either pine plantation or not pine plantation based on photointerpretation of high-resolution aerial imagery. A 1 km² moving window neighborhood analysis was then applied around each pixel in the state that: a) counted the number of pine pixels within the window, b) calculated the proportion of positive pine plantation points within the window versus the total number of points in the window, and c) calculated the proportion of negative points outside the window versus the total number of points outside the window. Lines representing the relationships between pine pixel counts and the two calculated proportions were plotted to determine their intersection point (630 pine pixels/1 km² neighborhood), which represented where the ability of the model to positively identify pine plantations and its ability to correctly identify areas that have pine cover and are not plantations was optimized. Consequently, the pine plantation mask was defined as LCMS vegetation change areas that corresponded to the threshold defined by the moving window analysis (1 km² neighborhoods with greater than 630 pine pixels). Within the pine plantation mask, we coded three different pine plantation types: 1) recently harvested, 2) young evergreen, and 3) mature. To do this, all changed pixels for the most recent 10 years were reclassified unless they were mapped as water, developed, wetland, or row crop types. Recently harvested plantations were comprised of reclassed pixels previously mapped as barren, deciduous woodland and forest, deciduous shrubland and young woodland, and herbaceous. Young evergreen pine plantations were comprised of pine woodland and forest, eastern redcedar woodland, and shrubland. Mature plantations were comprised of pixels mapped as pine or eastern redcedar woodland and shrubland that were coincident with sudden LCMS vegetation losses occurring between 1985—2011.

Ozark-Ouachita and Crowley’s Ridge Woodland and Forest Type Modeling - Soil eco-groups were used to model Dry, Dry-Mesic, and Mesic Forest types for the Ozark-Ouachita region. In addition, plot data from this project, the National Park Service (Pea Ridge National Military Park and Hot Springs National Park), and The Natural Conservancy (collected on USFS lands) were aggregated and analyzed. Each plot was labeled as Dry Oak Woodland, Dry-Mesic Oak Forest, or Mesic Hardwood Forest based on field calls. Values for topographic wetness index (TWI), soil moisture, and solar radiation were extracted for each of the plot locations. Solar radiation and TWI were calculated using traditional methods, whereas soil moisture was derived from a five-year average given by estimates from the SMAP-HydroBlocks dataset [44], which combines hyper-resolution land surface modeling with satellite derived soil moisture models. These three variables were used to build a decision tree to differentiate the three vegetation types using recursive partitioning. The model using the three variables was insufficient to differentiate mesic from the other two vegetation types, so plots representing mesic vegetation were excluded from further analysis. This resulted in using 452 dry plots and 646 dry-mesic plots to build a decision tree to separate these types. The result was a model that reclassified objects representing Ozark-Ouachita Dry-Mesic Oak Forest to Ozark-Ouachita Dry Oak Woodland if mean solar radiation for the object was > 60393 WH/sq m and mean soil moisture was < 0.211416 m³/m³. TWI did not prove useful in the differentiation of classes in the model and was therefore omitted. A similar classification tree methodology was used to differentiate Mississippi Alluvial Valley Dry-Mesic Loess Deciduous Forest from Crowley's Ridge Mesic Loess Slope Deciduous Forest in the Crowley’s Ridge ecoregion. Slope was used in conjunction with solar radiation such that objects with mean slope > 22% or with mean solar radiation < 59986 WH/sq m were reclassified to the Crowley's Ridge Mesic Loess Slope type.

To model glade types, glades previously mapped by Nelson [45]. Polygons provided by that project were rasterized (10-m pixels) to produce the glade types (limestone, sandstone, shale, etc.). To identify cedar glade types, the mapped glades were buffered by 200 meters and objects classified as cedar landcover within the buffer were selected. The objects were then dissolved to form larger polygons, and polygons that did not intersect the original mapped glade were eliminated. Polygons resulting were then classified according to their corresponding Nelson glade type.

Crosstimbers Oak Woodland and Forest and the River Floodplain Deep Sand Grassland-Woodland Complex ecological system types were classified by applying field expertise and aerial photography to reassign appropriate image objects. The Ouachita Montane Oak Forest type was mapped using a similar technique based on field expertise, aerial imagery, elevation, and proximity to the same type previously mapped in Oklahoma. The Mississippi Alluvial Valley Dry-Mesic Loess Eastern Redcedar Woodland type in the Crowley’s Ridge ecoregion was mapped in similar fashion.

2.4.4. Production of Final Ecological Type Map

From visual inspection, we noted that incorrect Herbaceous Wetland, Wooded Wetland, and Pine Woodland/Forest landcover often occurred in shadows in mountainous regions (e.g., Boston, Ouachita Mountains). To correct this issue, these classes were replaced where inundation frequency from 1983—2011 was zero [46] using the results of an alternate landcover classification that excluded Wet Herbaceous and Wooded Wetland as classification targets and excluded leaf-off imagery from the explanatory data stack.

Buildings and roads were burned in as a final modification to the statewide map of ecological types. Buildings were based on a rasterized (10-meter cell-size) version of the previously described Microsoft Buildings data layer, whereas roads were based on a multi-class rasterization (10-meter cell-size) of the Arkansas Centerline File, which depicts detailed road linework statewide. Buildings and paved roads were burned in as the urban impervious class, whereas minor, unpaved roads were burned in as mixed intensity urban.

2.5. Detailed Ecological Field Data

Field data were collected in order to assess the mapped ecological types. The Arkansas Department of Transportation coordinated field data collection efforts and a previously developed standardized methodology was employed [47]. The general data collection procedure included:

• Ecologists designed sampling efforts to cover major areas of different vegetation cover by geographic region, with a goal of collecting data in every county. Large areas of uniform landcover, especially row crops, were identified and data collection efforts were reduced in these areas.
• Sample virtual plots were collected on both sides of a road at approximately 1-mile intervals from a random starting point. Locations were geographically referenced using a GPS enabled computer in the field vehicle.
• Sample virtual plots consisted of visually homogeneous landcover within a 50-meter circular radius while maintaining adequate distance from roads to avoid disturbances.
• Air photos and digital soils maps loaded on the laptop were used to help facilitate data collection when field survey views were limited.

A standardized suite of data was collected using a digital feature data form that included drop-down lists for the data fields to reduce input errors. Photographs were taken at each plot location to serve as reference. Collected data included cover by cover class for the tree, shrub, and herbaceous vegetation layer, and the top three species in rank order by cover for each layer. Plant species names were obtained from the USDA Plants database (https://plants.sc.egov.usda.gov/java/). Finally, the landcover type and ecological type were recorded for each plot based on the expert opinion of field ecologists. A total of 1704 plots was collected.

3. Results

3.1. Statewide Landcover Results

Since the 10-meter landcover results for the study region were produced zonally, classification accuracy was assessed by evaluating the final RandomForest out-of-bag (OOB) scores [48] for each of the zones. This assessment included types that were mapped from the initial RandomForest landcover classification, rather than those that were subsequently modeled and improved using the approaches described in previous sections. Overall OOB error rates were 6.88% for zone 1, 6.91% for zone 2, and 5.27% for zone 3. Less prevalent classes tended to exhibit higher error rates relative to classes that comprised larger proportions of the respective classification zones, and these types of issues represent absolute (non-correctable) limitations of the remote sensing classification. The highest error rates were 27.81% for deciduous shrubland and young woodland in zone 1, due primarily to confusion with herbaceous vegetation and deciduous forest; 24.5% for eastern redcedar woodland and shrubland for zone 2 due primarily to confusion with pine woodland and forest; and 14.5% for deciduous shrubland and young woodland in zone 3, due primarily to confusion with herbaceous vegetation and deciduous woodland and forest.

3.2. Ecological Mapping Results

Application of the master and special cases models described above resulted in a total of 123 types mapped statewide for Arkansas (Figure 4). Some types such as flatwoods, saline soil woodlands, and floodplain woodlands may be hard to recognize in the field but were mapped using soils information in the ecological systems model. Many ecological system type concepts have similar aspects and species compositions such that vegetation ecologists may have different interpretations of the same location on the ground and identify different types. Despite these limitations, information from 1704 field-collected virtual plots were summarized according to mapped ecological types to assess the accuracy of the final statewide mapping results. Thus, 80.97% of the field plot sample type calls either perfectly agreed with the mapped EMS type or were a closely related type.

• Perfect Agreement between field vs. mapped type: 37.92%
• Close Agreement between field vs. mapped type: 43.05%
• Disagreement between field vs. mapped type due to modeling: 9.76%
• Disagreement between field vs. mapped type due to landcover difference: 9.27%

Of the 123 mapped ecological system types, 23 made up > 1% of the study area (Table 2). Ozark-Ouachita Dry-Mesic Forest and Row Crops each made up more than 17% of the area, with the latter restricted to the Ozark-Ouachita region, and the former most abundant in the Mississippi Alluvial Valley (Figure 4). Six of the 23 most abundant types were pine plantation in different regions and in different stages of maturity, and these covered 14.27% of the state. Large river floodplain types in the West Gulf Coastal Plain and Mississippi River Alluvial Plain covered 5.31% of the state.

4. Discussion

This study demonstrated an application in which Sentinel-2 satellite imagery can serve as an effective foundation for mapping ecosystem types when combined with other modeling approaches. We demonstrated how increasingly available national-scale ancillary datasets from sources such as the USGS, NRCS, and other researchers can be integrated into an ecological mapping framework to enhance ecological type classifications. The 10-meter resolution offered by the Sentinel-2 MSI platform also represents an ideal tradeoff between granularity and processing demands, allowing vegetation to be depicted more realistically across large areas without onerous computational requirements.

Future research should focus on leveraging cloud-hosted satellite imagery, such as the Sentinel-2 archive, and ancillary modeling frameworks such as those presented here, to develop regularly updated ecological type maps. Such maps can provide a solid basis for effective wildlife conservation, natural resource management, fire modeling, hydrologic modeling, and environmental preservation efforts, as well as serve as a useful tool for evaluating the outcomes of past management practices. There is also further opportunity to explore how other existing and emerging data, including improved soil properties datasets (e.g., POLARIS) [49] and spaceborne LiDAR products such as those provided by ICESat-2 and GEDI [50,51] can enhance ecological type modeling going forward.