Loss of Relict Oak Forests along Coastal Louisiana: A Multiyear Analysis Using Google Earth Engine

Abstract

Coastal forests along the southeastern Gulf of Mexico are known to be diminishing at an alarming rate. The live-oak dominant chenier forests of southeast Louisiana are amongst those exhibiting the steepest declines. The remnant stands have experienced numerous hurricanes and intense storm events in recent years, calling into question the current status and immediate future of this imperiled natural resource. Despite their noted ecological and physiographic importance, there is a lack within national geographic data repositories of accurate representations of forest loss and wetland extent for this region. Supervised machine learning algorithms in the Google Earth Engine were used to classify and process high-resolution National Agricultural Image Product (NAIP) datasets to create accurate (>90%) tree cover maps of the Louisiana Chenier Plains in Cameron and Vermilion Parishes. Data from three different years (2003, 2007, and 2019) were used to map 2302 km² along the southwestern coast of Louisiana. According to the analyses, there was a 35.73% loss of forest cover in this region between 2003 and 2019. A majority of the land-use change was from tree cover to saltmarsh, with losses in pastoral land also documented. We found variable rates of loss with respect to elevation. Forest cover losses corresponded strongly to rises in mean sea level. These findings deliver a baseline understanding of the rate of forest loss in this region, highlighting the reduction and potentially the eventual extirpation of this imperiled ecosystem.

1. Introduction

Coastal forests are exposed to a unique suite of disturbances. Rising sea levels, human encroachment, salinification of soils, and exposure to extreme weather, all of which are intricately tied to global change, are leading to irreversible changes in these environments [1,2,3]. Cheniers and their associated plains are products of sediment outwash and fluvial sedimentation from major river systems along coastal regions. They are found globally across South America, Europe, North America, Africa, Australia, and Asia [4,5,6]. Our study focuses on cheniers along the Texas–Louisiana Coastal Marsh Ecoregion formed from the Mississippi River within the last 3000 years [7,8]. These cheniers once contained a rich diversity of closed-canopy hardwood species, the dominant being the coastal live oak (Quercus virginiana Mill.) [9]. Colloquially, the term ‘chêne’ translates to “place of oaks” in Acadian French, which correctly recognizes the stature and prevalence of this wide-spreading tree as a keystone member of the ecosystem. Agrarian settlement and infrastructure development has driven land-use change to non-forest means from the early 1900s to the present day [10]. In this study, we used remotely sensed imagery to census forest area and assess rates of loss in tree cover from degradative natural and anthropogenic influences.

The National Land Cover, National Wetland Inventory, Coastal Change Analysis Program, and Global Forest Cover datasets (PALSAR-2) have failed to identify these forests, due to the coarse nature of classifying imagery [11]. As we continue to see rapid declines in forest areas along coastlines, it is appropriate to recognize and inventory areas not previously considered in the analysis of forest loss in the southeastern United States (US). A lack of inclusion of this land-cover type has led to an underrepresentation of over 2790 ha. of forested land area. Existing data for southeastern Louisiana, the Floridian Gulf Coastal forests, and coastal forests in North and South Carolina all provide estimations of areal forest and marsh extent from the mid-1980s to the present day [12,13,14]. In comparison, this area is data-poor regarding land use and change, limiting the ability of governmental agencies or non-governmental organizations to draft policies that include this forest ecosystem. Forest change modelling can be done using a variety of methods, from spectral to spatial, dependent on the imaging sensor characteristics (multispectral, hyperspectral, SAR, etc.). Multi-band imagery can be utilized to create novel vegetation indices that can map forest change successfully at variety of scales (30 m–1 km pixel size) [15]. Sub-pixel mapping allows investigators to create geospatial models that are more highly granular than the original pixels, by using auxiliary datasets. Sub-pixel methodologies, such as spectral unmixing or sub-pixel swaps, present opportunities to create more refined maps. Valjarević et al. [16] used historic topographic maps to create socio-economic associations between forest cover, population abundance, and average rates of deforestation. Spatially based methodologies can be done using clustering, grouping, or organizing of pixels into fixed arrays, with regard to the similarity or lack thereof to neighboring pixels. Cluster mapping such as Fuzzy K-Means, fusion algorithms, or various machine-learning algorithms (RF, CNN, SVM) can perform very accurate (>90%) classification of a variety of remote sensing products [17,18].

The preceding studies show the importance of “big data” and remote sensing in functional ecological studies. Big datasets can present difficulties for analyses, as their volume prevents traditional tools from accessing underlying information of value [19]. Large time-series datasets (i.e., satellite imagery) are a powerful way to exemplify changes at the global or local scale for any vegetated land use [20]. Very high resolution (VHR) imagery (<2 m) extends the benefit by providing the ability to clearly demarcate change and assess fine patterns of variation in land use and land-cover [21]. The VHR National Image Agricultural Product (NAIP) used in this study is voluminous; a coterminous orthoimage of the US would be over 65TB [22]. For researchers without access to super-computers, cloud-based computing increases accessibility to these data. One can access a wide variety of data sources or iteratively test algorithms in the creation of a processing workflow, something that could take days or weeks if processed locally [23,24]. In this study, we used Google Earth Engine (GEE) to tackle this large-scale mapping and analysis, intuitively stitching and classifying custom extents of desired study areas.

The geomorphology of Louisiana’s coastline has been in a state of intense flux over the last century, impacting biotic components which rely on stable stratigraphy and hydrology [25]. Much of the research done on this dilemma has focused on dynamics related to land loss or on the marsh aspects of the ecosystem [26,27]. Neyland and Meyer [28] carried out the first study on the composition, extent, and ecology of the Louisiana Chenier Plain woodlands in 1997 yet did not give estimates of the relative abundance of remnant forests. The forests of the Louisiana Chenier Plain were decimated before being properly inventoried; historic plant communities and abundances are largely unknown [29]. Gosselink [9] estimates a historic extent of 46,500 ha. forested land-cover between cheniers, spoil banks, and bottomland swamp forests. Our results show a significant change in the area and location of forest cover today.

Largely, climactic factors and sea level have determined the fate of coastal forests, and by extension the whole region [30,31]. Extreme weather is also a growing concern. Since the year 2000, the chenier woodlands have suffered the wrath of hurricanes Rita, Gustav, Laura, and Delta. These storms were all category 2 or greater when making landfall on the coast. Each hurricane brought the dual threats of storm surge, which can cause inundation up to 4m above ground level, and structural damage from wind-gusts up to 177 km/h [32]. There is an eminent need to quantify the impacts on the vegetation of coastal forests, as tropical storm events are likely to increase in frequency and intensity in the near future [33]. Our observations help establish baseline data on the current extent of these unique and vulnerable forest ecosystems within the region, making clear their status and the vulnerability this ecosystem.

2. Materials and Methods

2.1. Study Area

The study area encompassed the coastline of Cameron and Vermilion Parishes in southwestern Louisiana (Figure 1). A varying buffer of 6 to 14 km from the coast was maintained to capture pertinent physiographic features. The area is classified as a low-profile, microtidal, storm-dominated coast that experiences both erosive and accumulative sediment loads [8]. It is distinctive in that the land experiences subsidence in addition to above-average rates of mean sea level (MSL) rise [34]. Broadly, this area is divided into two geomorphic regions, the Chenier Plains and the marginal deltaic plains of SE Louisiana. The chenier plain ridges are laterally accreted sections of sand and shells deposited episodically. Accretion heights vary between 0–4 m, with lengths ranging from 10–90 km and widths less than 1 km across, arrayed in a linear-to-concave manner [20]. These ridges form complexes, defined as two or more parallel or sub-parallel sections, that stretch for lengths well beyond 400 km. Chenier complexes can directly about the coast, while the most landward are located up to 15 km away. The two soil complexes of the region are Hackberry and Mermenteau [35]. Hackberry soils form the ridges and are deep, poorly drained, and moderately permeable, composed of sand and shells. The Mermenteau soil complex forms the edaphic base of the saltmarshes. The complex is very deep, with reduced permeability, and is primarily composed of clayey over loamy sediment. Contextually, these ridges form the only topography above MSL within this region, and are thus visually striking when seen at distance.

The extant vegetation communities form a closed-canopy hardwood forest composed of live oak (Quercus virginiana Mill.), southern hackberry (Celtis laevigata (Kunth) Spreng., American elm (Ulmus americana L.), black willow (Salix nigra Marshall.), honey locust (Gleditsia triacanthos L.), and persimmon (Diospyros virginiana L.), among others, as well as now the increasingly common invasive Chinese tallow (Triadica sebifera L.). Shrubs such as American holly (Ilex opaca Aiton.), coastal mallow (Kosteletzkya pentacarpos L.), common hawthorn (Crataegus viridis L.), toothache-tree (Zanthoxylum americanum Mill.), and chinaberry (Melia azedarach L.), among many more, are also frequently encountered [22]. The intact community forms a dense stand that prohibits the establishment of shade-intolerant species. These communities have been described as “cathedral-like” due to the tall and wide branching stature of the live oaks. The leaves of live oaks are dark green in hue and maintain foliage year-round, making their foliage a strong visual indicator in aerial imagery.

Recent vegetation surveys show a shift to a community more akin to an open saltmarsh swamp with thick stands of southern amaranth (Amaranthus australis (A.Gray) J.D.Sauer), roseau cane (Phragmites australis (Cav.) Steud.), narrow-leaved cattail (Typha angustifolia L.), and other wetland obligate species interspersed with dead or dying trees (unpublished data, Bhattacharjee, J). This newly emerging ecotype represents a breakdown of the traditional segregation of upland forest and lowland haline swamp. This change is broadly represented by the encroachment of hydrophytic vegetation becoming dominant in the understory. This change is identifiable when seen in VHR imagery, although the spectral characteristics of the scene can appear differently. Much of the structural complexity of the forest has been lost. The physiognomy appears as a brushy single layer with limited overstorey, mirroring the surrounding saltmarsh to some degree. Some areas have become entirely dominated by dense herbaceous vegetation, preventing the establishment of native upper-story species. The site has always experienced anthropogenic influence in community composition; the native Attakapas people contributed to the original diversity of vegetation by introducing and maintaining an abundance of food-bearing plants within the landscape [36,37]. The study site does not have any mangrove swamps, which can be found further east and south along the Gulf Coast.

2.2. Data Sources

The National Agricultural Image Program of the United States Department of Agriculture (USDA) houses aerial orthoimage sets for the contiguous United States, in 3.75° × 3.75° quadrangles. The 2003 and 2007 imagery for Louisiana were collected at 2 m and 1 m Ground Sample Distance (GSD), respectively, in red, green, and blue bands. The 2019 imagery has 0.5 m GSD with the addition of a near-infrared band. In total, 234 NAIP image tiles were used for analysis. Finer scale imagery was not downsampled to coarser resolution. Additionally, all bands were used in the classification scheme to provide the greatest discriminatory accuracy. An effort was made to classify imagery from around a similar date to capture equivalent phenophases, although image quality was given precedence over temporal constraints. The aerial images used in the data analysis pipeline for the years 2003, 2007, and 2019 were collected on 8/24/03, 8/10/07, and 9/5/19, respectively. These particular years were chosen as they represented the highest quality imagery within the study period. It should be noted the NAIP is not an annual dataset, collections vary according to agreements with state governments and are often performed on a biennial or triennial cycle. Limiting factors in many of the interior year image sets were cloud cover or high local variance in scene reflectance, which were prohibitive for the manner of land cover classification used in the study. Intense glare from surrounding water bodies caused significant reduction in quality of less optimal image correction after the fact. It became difficult to distinguishing between many similar vegetation types, as even moderate local variance can introduce noise, degrading algorithm effectiveness. Imagery from 2013 had begun to be processed but was abandoned after initial development due to scene variability and poor classification performance. Only image sets with less than 5% cloud cover were used.

The Coastal National Elevation Database (CoNED) Digital Elevation Model (DEM) was used for elevation analysis. The CoNED is generated by the United States Geological Survey (USGS) using LiDAR and has a 3 m GSD. It is housed in the National Map Data repository. The raster was downsampled to 1 m to match the resolution of the GEE classification output. Mean sea level (MSL), salinity, and temperature values were extracted from two USGS stream gages (#08017118, #07387050) located at opposite ends of the study area, providing data for the entirety of our study period. The data were extracted from the USGS StreamStats interactive web map. The National Hydrography Dataset (NHD), a vector repository of waterbodies (NHDWaterbody, NHDArea) and wetland features (NHDWaterbody FCode: 46600, 46601, 46602), was used for feature input in the creation of post-processing masks. The NHD is also a product of the USGS.

2.3. Classification and Processing Algorithms

Google Earth Engine offers a suite of supervised and unsupervised Machine Learning Classifiers for image classification. Through iterative testing of all available options, the Classification and Regression Tree (CART) algorithm was chosen for this study. CART is a decision-tree-based supervised learning model able to handle nonlinear relationships among features [38]. Because class balancing is recommended for this algorithm, a stable ratio of tree-to-pasture and tree-to-marsh training sets were created. Parametrization was not modified from default properties. The primary consideration in choosing this algorithm was robust classification for forest cover. Visually, we observed that many other contenders had higher incidence rates of non-forest-based classification of tree cover when conducting initial testing, leading to a general under-representation of forest cover. We use a reductionist approach for accuracy in this study, where pixel values for forest cover were reduced rather than added, which makes underrepresentation problematic. The post-processing vectors allowed extremely fine detail in the modification of the original CART output.

The study area was divided into four sections (A, B, C, and D) of equal approximate length, as seen in Figure 1. Sections C and D were further divided into three subsections due to image variability in the 2003 dataset, and were merged post-classification. This facilitated a more targeted training dataset given local variation, while also expediting processing within the workflow. We used six total land-cover classes in the classification scheme—tree, pasture, pond, river, saltmarsh, and impervious area. The classes were comprised of manually delineated Point FeatureCollections within GEE and were given specific unique IDs. Tree cover was interpreted as the visual presence of dark foliage, branching, or scrubbiness. As such, these selections include trees located on spoil bank levees (Figure 2B), swamp forests (Figure 2C), or residential areas, in addition to the live-oak communities found on the main chenier ridges (Figure 2A,D). The rationale was to include any type of vegetation that promotes beneficial land use and habitat. Pasture was interpreted as any area in which there are grazing activities, or managed lands on which regular mowing occurs, such as yards or parks. The Point FeatureCollections were created using 2019 NAIP imagery and were manually validated before use in each preceding year to ensure classes remained accurate. Modifications were made when class values changed between years. The CART classification was subject to a 3-pixel square image reducer, thereby homogenizing values to reduce the “salt and pepper effect” common in pixelwise classifications [39]. The data were then exported to Google Drive for further processing in Quantum GIS (QGIS). A detailed workflow for the procedure is given in Figure 3.

2.4. Post-Processing of GEE Outputs

Two-layer masks created from the NHD were used to refine the classified output. The masks were originally ascertained and modified in a vector format, and were converted to a raster-type format when used for modifying the GEE outputs. The two layers, which were in the same projection and scale, were employed to perform mathematical operations on a per-pixel basis. Raster calculations can be applied to overlapping areas according to their underlying geospatial properties, affecting how data is represented [40]. As the NHD is incomplete regarding wetlands and pond features, care was taken to digitize areas of the 2003 NAIP imagery that were wholly saltmarsh or water [41]. An assumption was made that any saltmarsh or water feature present in 2003 would maintain its non-forested land use during the study period. This allowed a set standard for the modification of GEE outputs. The saltmarsh mask was set so that any pixel registering as a tree would be converted into saltmarsh, without the modification of any other land-cover class. The water mask was absolute in terms of pixel modification. All pixels overlapping the mask were converted, regardless of their class. Tree cover vectors were then extracted from each year’s classification. All cover area calculations, including different elevation products, were performed on this vectorized form, due to the superior computational ability for calculating absolute area of georeferenced polygons in a vector format. Polygons were also manually scanned over the respective years of NAIP imagery and were either deleted or adjusted following a set standard for Quality Assurance (QA) to ensure best representation.

2.5. Data Analysis and Statistics

Because post-processing and QA significantly modified the original CART classification, baseline statistics on algorithm confidence would not be representative. Instead, producer’s accuracy (PA) was used as the primary determinant of classification success (Equation (1), where reference sites denote created validation points).

$$\text{Producer’s Accuracy}=\frac{\#\mathrm{of}\mathrm{correctly}\mathrm{classified}\mathrm{reference}\mathrm{sites}}{\#\mathrm{of}\mathrm{total}\mathrm{reference}\mathrm{sites}}\text{×}100$$

(1)

The creation of points for sections A, B, C, and D, to assess accuracy, was carried out using three tools in the QGIS Processing Toolbox. First, a ‘Difference’ was run to extract tree cover in an original year’s dataset but not in the succeeding year. A ‘Dissolve’ was run on this product to create one large MultiPolygon object. Finally, ‘Random Points’ inside the Polygon function was used to create 200 point features that were manually validated. Three columns indicating original year accuracy (Y/N), whether land-use change was correctly attributed (Y/N), and what it had changed to (six classes), were populated for each point. The last metric had absolute accuracy irrespective of the year’s classification, as it was assigned visually. Validation was carried out against the respective year’s NAIP dataset. Because we did not assess land-use change in 2019, only accuracy was computed for that. A Cohen’s Kappa value was used to assess interrater reliability using Equation (2), where Pr(a) represents the actual observed agreement and Pr(e) represents chance agreement. The equations were derived from and expanded on by McHugh [42].

$$\mathsf{\kappa }=\frac{\mathrm{Pr}\left(a\right)-\mathrm{Pr}\left(e\right)}{1-\mathrm{Pr}\left(e\right)}$$

(2)

The CoNED DEM was subdivided into 0–1, 1–2, 2–3, and 3+ meter classes using the Raster Calculator tool in QGIS. The class values were then vectorized to run a geospatial difference iteratively with each year’s tree cover vectors. The difference was validated and then used for area calculations regarding differential rates of loss for elevation gradients. Gage height and temperature values were extracted from both stream gages and were combined into a site average. We used gage height as a proxy for MSL due to the spatial proximity to the coast and study sites, as well as for considering the response to freshwater influence, which is an important hydrological component of the Chenier Plains. Average values for the three wettest months (July, August, and September) were used in the regressions and analysis.

3. Results

3.1. Forest Loss, Transition, and Classification Accuracy

From 2003 to 2019, an overall decline of 35.73% was seen for forest cover and 22.38% for pasture. The annualized rate of loss for forest cover was approximately 2% higher within the period 2003–2007, losing 4.2% of tree cover per year assuming a constant change (

Table 1. The amount and percentage loss of three land-cover types—tree, marsh, and pasture—over each section and the total area spanning two parishes, Calcasieu and Vermillion, on the western Louisiana Gulf Coast. Sections A, B, C, and D represent areas delineated for this study within the two parishes (for details see Figure 1).

Landcover	2003 (ha.)	2007 (ha.)	2019 (ha.)	Net Change (%) (2003–2019)
Tree
A	462.74	376.39	179.88	−61.13
B	643.89	386.34	274.86	−57.31
C	996.23	827.91	650.50	−34.70
D	2317.75	2151.92	1735.80	−25.11
Total	4420.61	3742.56	2841.04	−35.73
Marsh
A	11,822.44	11,480.60	12,386.36	4.77
B	11,130.45	10,847.38	13,954.53	25.37
C	20,517.30	22,777.85	22,044.12	7.44
D	65,986.32	68,393.19	63,599.52	−3.62
Total	109,456.51	113,499.03	111,984.53	2.31
Pasture
A	5101.46	4603.23	3243.48	−36.42
B	1363.95	890.19	770.21	−43.53
C	6962.09	7545.00	6530.27	−6.20
D	11,520.47	10,360.39	8819.60	−23.44
Total	24,947.96	23,398.81	19,363.55	−22.38

). Losses were greatest in the western two reaches of the Chenier Plains, averaging over 21.3% more than sections C or D. The largest amount of tree cover was located within section D, with significant contributions from spoil banks. Chenier Au Tigre (location: −92.20665, 29.56885) had the largest total tree cover within the study area, albeit spread over a large expanse. Annual loss rates were more constant in sections containing primarily chenier ridge forests, in contrast to swamp or spoil bank forests. Section D has a large saltmarsh (1–10 km) bisecting the chenier from the coast, which contributed to a larger land and forest area. There was a cumulative gain of 2.31% for the saltmarshes, adding 25.28 km². In the 2003–2007 period, the only gains in land cover type were in saltmarsh which was maintained into 2019 (Figure 4).

We surveyed a total of 2123 points across the study area. Producer’s accuracy was 91.6%, 81.6%, and 91.8% for the years 2003, 2007, and 2019, respectively. The Cohen’s Kappa indexes of 0.91, 0.82, and 0.93, for 2003, 2007, and 2019, respectively, showed strong agreement between raters across all years. We considered any value above 0.85 to be a valid assessment. The year 2007 scored lower in both PA and Kappa values, likely due to the nature of the base imagery, however, it still fell into an acceptable range. Extreme contrast and saturation in the 2007 NAIP imagery caused larger aggregations of tree cover pixels, as grasses and trees have a similar hue, accentuated by the pixel reduction post-classification. Validation points that were very close to tree cover but not categorized as such were not considered accurate, which contributed to this source of error. The 2019 classifications were generalized as being more correct from the onset, with little need for QA.

Changes in land-cover were not evenly distributed across different sections, as can be seen in Table 1. Typically, cheniers are more proximally located to the coast within the eastern reaches, with depositional bars spreading the length. One noted outlier in the assessment of change was observed from 2003 to 2007 in section B, where the loss of pasture accounted for over 65.6%. The quantitative measurements for changes in land cover were also ascertained from these samples. We observed changes of 60.6%, 28.2%, 3.4%, and 7.6% for saltmarsh, pasture, water, and impervious area, respectively.

Observed results deviated from hypothesized values, specifically in the large proportion of conversion to pasture. A lone tree existing one year after transitioning to an open field was common across all areas. It was assumed the percentage of conversion from tree to saltmarsh would be even greater. However, impervious area transitions were also high, given that considerable forest cover directly adjoins roads and houses in this area. The underlying structure can then be seen following the mortality of vegetation. Change to water was also observed for several reasons; in areas along waterbodies, a loss of tree cover exposed the water, while severe areas of deforestation prompted a change to a more wetland hydrology, which also required water replacement pixels.

3.2. Rates of Loss with Respect to Elevation

A net decrease in forest cover was seen across all elevation categories (

Table 2. Differences in the amount of forest area within each elevation band and the associated percentage loss of forest cover between years. Sections A, B, C, and D represent areas delineated for this study within the two parishes (for details see Figure 1).

	Sites and Forest Cover (ha.)					Percentage Loss
	Year—2003
Elevation (m)	A	B	C	D	Total:
0–1	130.83	255.25	663.34	1872.36	2921.78
1–2	271.68	363.74	255.83	390.48	1281.72
2–3	19.34	21.88	64.81	46.40	152.42
>3	3.61	1.70	0.05	3.43	8.79
	Year—2007
	A	B	C	D	Total:	2003–2007
0–1	183.71	165.02	563.49	1667.76	2579.97	−11.70
1–2	178.34	202.10	212.89	42.05	635.38	−50.43
2–3	12.28	17.83	41.18	54.46	125.75	−17.50
>3	0.48	2.08	0.02	4.04	6.62	−24.74
	Year—2019
	A	B	C	D	Total:	2003–2019
0–1	71.61	82.27	377.62	1421.02	1952.53	−33.17
1–2	102.83	129.13	246.56	44.01	522.54	−59.23
2–3	15.06	9.32	57.89	42.53	124.80	−18.12
>3	0.45	1.41	0.02	4.02	5.90	−32.88

, Figure 5). The largest losses were in the 1–2 m category, with a net reduction of 59.3% from 2003 to 2019. There was a temporal shift in the rate of change, with the greatest amount of loss between 2003–2007 contributed by section D. Values within the 2–3 m category remained steady from 2007 to 2019, with a reduction of only 0.76%. The large majority (74.9%) of the remaining forest was seen to be located within the 0–1 m band. The losses were not evenly distributed within this band, buffered to some degree by section D. A comparison of sections A vs. C shows a decrease of 49.3%, with losses centered around the period 2007–2019. One issue to note is that the GSD for the DEM was 3 m, which is wider than some spoil banks. While the local elevation may in fact rise to 2–3 m, the value was reduced due to the effect of averaging. The losses in the 3 m+ band are discussed below; a note should be made of the limited representation. Overall, the losses varied both spatially and temporally across the region, largely separated into eastern and western reaches.

3.3. Interaction with Water Level

Water plays a critical role in the natural function of this ecosystem. As such, we considered the interplay of temperature and MSL over a decadal time scale, and the local rates of temperature and gage height from both USGS stations (Figure 6). In Figure 6, forest extent is superimposed for each year’s classification to show visual correspondence between the two drivers. Temperature and rise in sea-level were negatively correlated until an inflection point around 2010, after which there was a relatively stronger correlation, increasing to a maximum of 0.25. The correlative increase lends evidence to the correspondence between climate and water pressure in coastal areas. An increase in the average MSL has been observed for the site. The two stream gages showed an annualized MSL rise of 13.6 mm/year for Vermilion Bay and 5.9 mm/year for the Calcasieu River. There was a weak but positive trend for temperature to increase slightly from the onset of the study period. Time proved to be a good indicator for increase in water level and for decrease in forest cover (Figure 6). The lack of tree-cover samples prohibited many statistical comparisons, but visually the correspondence seems evident.

4. Discussion

4.1. Classification Reliability and Constraints

NAIP imagery proved efficacious for fine-scale analysis in this region, shown by the strong average producer’s accuracy and Kappa coefficient. This is a positive indicator that variability between years or areas did not influence classification accuracies to a large extent, and that the final classification is a reliable interpretation [42]. Many other studies conducting land-cover classifications using NAIP imagery have also shown robust accuracies [43,44,45]. However, there were noted issues seen in certain parts of the workflow. A problem preeminent in the initial classification was a limitation of spectral resolution. This problem became pronounced in areas with similar shades of green, or between dark blues and black hues. Without the use of ancillary datasets, the classification accuracy was poor between dark water and shadows of trees, with considerable overlap (misclassifications of approximately 50%). Mixed land use can also prove problematic in NAIP land-cover classifications [46]. This became especially pronounced in transitional forest/marsh areas in which there was no dominant spectral signature for a particular land use. Local variation between adjoining image sets highlighted another shortcoming of NAIP imagery. The viewing geometry, illumination, and time between the capture of neighboring image tiles created edge effects, including artificially high or low pixel values, and refraction problems on water [47]. Each year’s dataset was collected by an external contractor, causing differences in sensor instrumentation and parametrization between years. By 2019 many of these problems had been mitigated, and images were of much higher quality.

Compared to a large swath sensor such as Landsat-8 or Sentinel-2, an innate non-uniformity poses problems for pixel-wise classifications of VHR imagery. Typically, data from satellite imagery have been the primary choice for analyzing coastal forests, marshes, and forests at large [48,49,50]. The spaceborne sensors obtain information over eight bands, allowing much finer spectral discrimination of certain classes of objects [51]. Even amongst diverse forest types, space-borne Sentinel-2 and ASTER imagery has been used to map distinct species clusters of broadleaf and deciduous plant species. The large spatial scale enables the use of these created geospatial products to analyze the efficacy of governmental reforestation programs in remote areas of the world [52]. The usage of normalized difference indices, including the Normalized Difference Vegetation Index (NDVI), Enhanced Vegetation Index (EVI), and Enhanced Water Index (EWI), can also provide greater confidence in the separation of water and marsh features from deciduous broadleaf forests [53]. However, a spatial resolution of 10m is still too coarse for the detection of remnant and landscape trees, especially those alongside homogenous flat features such as maintained grass or water bodies. The most practical solution would be to assess this land with a higher resolution dataset, such as unmanned aerial systems (UAS) imagery, which can capture a much finer GSD in a plethora of bands if multispectral or hyperspectral sensors are used [54]. In fact, scientists are now pushing the boundaries of the scale of analysis conducted through remote sensing. Broad groupings (e.g., clades or forest types) have been the historic norm in the classification of imagery for use in land management. A greater mixture of classes allows more accurate measurements for those interested in landscape-scale ecological processes. However, VHR can be used to map stands at the species level, which allows forest managers an unparalleled insight into the dynamics of recruitment, disease, or resource availability [55].

Google Earth Engine allowed us to utilize NAIP imagery at scale, which had for some time been a limiting factor in processing this data source locally. Even so, some limitations hindered the functional utility of GEE for detailed land-cover classification. In our relatively uniform environment, we discriminated between many shades of green. When viewing geometries cause variations in hue, different classes begin to overlap within the same spectral regions. Solutions using Geographic Object-Based Image Analysis (GEOBIA) can circumvent this issue, yet this segmentation algorithm is currently not supported by GEE. Many consider GEOBIA the optimal tool for segmenting landscapes and for the production of land-cover classifications from this source of imagery [38,56]. Its successful integration could substantially elevate the use and utility of mapping within GEE. Additionally, the allocation of increased processing bandwidth on GEE servers may help overcome some limits in classification. Greater bandwidth allows larger and more diverse training sets, in addition to more complex parameterizations of machine learning classifiers. If one is proceeding on local workstations, the analysis of NAIP imagery requires an intimate knowledge of data structures and process optimization. Custom architectures, processing pipelines, and simpler binary classifications can allow an expansion of the study area to cover coterminous areas on a larger scale [16].

4.2. Forest Loss Contextually

We mapped forest cover change dynamics between 2003–2019 on the southwest coast of Louisiana. 1580 hectares of tree cover, 1.3% of the terrestrial area, was converted to different land use. The total loss of forest cover in the study period was 35.73% over sixteen years, which is indicative of a larger rate of forest loss than at sites in a similar physio-geographic region, yet is certainly not a global outlier [18]. These findings validate the sentiments of people in the region and underscore a major irreversible trend in the area. Our estimate should be considered conservative, due to the way we approached the definition of “forest cover” within the landscape. Abandoned pastoral land exhibiting dense scrubbiness and shrub cover was demarcated as forest cover, even though the target species (live oak) was probably not present. Spoil banks, which are the byproduct of historic oil exploration in the region, indicating possible mortality, were also often considered a component of forest cover if the imagery could not indicate the absence of high vegetation. The primary method of identifying whether or not forest cover was present was the observation of a shadow adjoining the vegetation, indicating a vertical structure.

The Big Bend region of Florida, one of the longest undisturbed coastlines in the continental US, offers some insight into the vegetation dynamics this region has been experiencing. It has been the focus of numerous ecological studies, has a similar climate, and is exposed to storms with similar regularity. Raabe and Stumpf [57] found that from the early to late 1900s over 82 km² of forest land was converted to marsh, and 66 km² was in a transitional state to marsh. McCarthy et al. [58] document an annual forest cover loss of 7.44% from 2010–2017, an acceleration in the trend within that area. A new study by the same team found rates of loss totaling approximately 10 km² a year, 800% larger than previously stated in the literature [59]. Rates of annual mortality ranged from 4%–16% in Southeastern Louisiana from 1987 to 1997 [60]. In comparison, the annualized rate of loss found in our study never exceeded 5%. Plots located in South Carolina that have been less significantly impacted by MSL rise show rates of loss closer to 1%–2%, in line with our observed rate of loss from 2007–2019. However, we cannot assume that percentages of forest loss across the country have been constant within the last century, therefore it is critical to carry out intensive geospatial studies at regular intervals.

In the adjoining southeastern part of the state of Louisiana, in the Atchafalaya and Terrebonne basins, forested wetlands occupy 121,871 hectares or 41% of all vegetated land use [12]. In comparison, the terrestrial forested land use in our study area was 2.1% in 2019. The wide variation can be interpreted in several ways, each with important implications. Historic rates of deforestation differ in the two regions, with a greater emphasis on the preservation of wilderness areas in the Atchafalaya basin [61]. The Atchafalaya and Terrebonne basins also have productive land used in the growth of timber, which promotes long-term forest cover. The regeneration of forests in our study region also differs substantially. Chenier ecosystems and traditional alluvial bottomland hardwood forests are compositionally different in vegetation and soil, even when accounting for their unique hydrologic regimes. Communities of live oak may be less resilient to the types of forcings that cause changes in community composition within these ecosystems, furthered by proximity to the sea. If stocking density rates are similar to the adjoining basin, as they were in the past, we would see an areal extent well over 50,000 hectares.

The analysis vis-à-vis elevation shows weak trends in the manners of loss on elevational gradients. Logically, in the paradigm of coastal retreat, there should be an aggregation of forest resources towards the highest bands of topography in a region. Our analysis does not follow this chain of thought, with the largest losses observed in the 1–2 m category rather than the lowest 0–1 m (Table 2: 2007 C, 2019 D). While we do see a large loss by percentage in the highest elevation band, it is important to note that a small sample causes inherent variance. Forest cover on the highest reaches of the chenier ridges is typically less than 1% of total forest cover for that region. The reduction can be attributed to the loss of 23 features in section A accounting for 3.3 hectares bordering the Sabine River (Figure 1: −93.86286, 29.73355). The loss of forest cover was driven by intense weather. Vegetation in the area was completely destroyed, and substantial amounts of sedimentation could be observed, as evident in the 2007 NAIP imagery. The only gain in forest cover was seen in the 2–3 m elevation band in the period 2007, in an area that had a high density of pastoral land that was abandoned. Forest retreat is often discussed as one method of buffering forest loss along coasts [62]. Unfortunately, the geomorphology of this area does not allow great spatial movement as the current extent is already located within the highest bands of elevation.

4.3. Land Use Conversion and Driving Factors

The leading cause of land-use conversion from tree cover was to a saltmarsh ecotype (60.6%), often following mortality and replacement of vegetative biomass. The loss of forest cover can largely be typified into two categories. The ingression of saltwater causing salinification within the soil matrix paired with increased flooding caused mortality of overstorey vegetation over a multi-year timescale [63,64]. Hurricanes or other storm events caused rapid mortality through wind-throws, severe defoliation, saltwater inundation during a storm surge, or canopy damage. Following intense weather events between 2003 and 2007, the vast majority of cheniers located within a 1km buffer of the coast were unable to support tree cover, as seen in Figure 7A,B. Clearing for industry or development also caused immediate shifts in forest composition, as can be seen in Figure 8A,B. This change was more pronounced on reaches A and B, especially on Hackberry Beach. This area was the site of landfall for Hurricane Rita, which caused massive devastation to other components of the ecosystem [65,66]. Singular storm events have a transformative impact on vegetation communities, and the weight of this driver on annual rates of forest loss is a topic that should be explored further. Wayward shifts in species composition are a trend that has been observed globally, especially in forests that border some source of water. Trees like Juniperus virginiana L. that have distinctive water-use patterns in comparison to other coastal tree species are susceptible to salinification at much faster rates, as evidenced through isotope analyses [67]. As with many ecological systems, the entangled effects of climate and vegetation require diverse means of looking at change. Isotope analyses, dendrochronology, field sensors, and remote sensing data are needed to characterize land-use change at the stand and community levels.

Large forest remnants such as Cow Island (−92.83581, 29.73842), Tiger Island (−92.78439, 29.72832), and Grand Chenier Ridge (−92.81175, 29.72986), among others, have diminished more slowly. By 2019, these stands had lost an average of 15% of their maximum extent. The land use changed to both marsh and managed pasture, largely dependent on the surrounding land use. The conversion to managed pasture (28.2%) accounted for a larger than expected amount of land-cover change. Many of these changes were in areas where tree cover may have been more “ornamental” in nature. Trees beside houses or alongside roads are some of the largest specimens of live-oak to be found. However, these trees are often solitary and perhaps more vulnerable to wind-borne disturbance. The concern is that their loss due to abiotic drivers may not be easily mitigated, leading to a permanent decline. Our data support evidence on the whittling of forest remnants and the eventual fate of these trees becoming ‘relict forests’ unable to recruit new generations in the short term; and in the long-term, total disappearance of these cheniers along the Louisiana Gulf Coast [68].

Beyond loss in forest cover, there is a noticeable trend towards an altered hydrology where pastoral areas and intact forest remnants are associated with pockets of standing water. An increase in the correlative factor between temperature and water levels could be a worrying trend; a direct conversion of tree cover to water at 3.4% could be an unforeseen outcome. Chenier ridges in their geomorphic context should not support non-ephemeral ponding. It is not documented how increased tidal influence, which is one factor of rising MSL, influences the hydrology of this ecosystem and its associated communities. Albeit in a more inland context, Magolan and Halls [69] documented a transition of over 40% of land to emergent marsh from closed-canopy cover from 1949 to 2018 due to hydrologic influence. Recent rates of relative sea rise have been documented as high as 11.9 mm/year compared to background Holocene rates of 1.7 mm/year [70]. This has resulted in a rate of lateral coastal forest retreat between 2–14 times greater than pre-industrial levels [26]. The rate of shoreline retreat in this region increased from ~6.1 to ~9.4 m/year, based on the 100- and 30-year moving averages, amounting to over 10,000 hectares of land loss every year from 1800 to 2005 [27].

4.4. Future of Chenier Forests

The chenier complexes of Louisiana do not have the natural benefaction of maintaining large original tracts untouched by axe or dozer. The extent map through the years (Figure 9) exemplifies the scattered nature of forest remnants. The largest single continuous chenier forest (Grand Chenier Ridge—18.5 ha.) encompasses less than 1% of the total forest area. Exposed forest edges pose a problem for the continued integrity of the stand as adapted species fail to compete in these new microclimates [71]. Slow-growing saplings of desirable tree species (live oak, hackberry) are subject to unfavorable growing conditions from weather and light availability. Vegetation surveys following aspect and inundation along the Florida coast have all but confirmed the eventual extirpation of overstory vegetation through poor recruitment of native species [31]. Chinese tallow (Triadica sebifera), a non-native tree, is likely to fill empty niche space. This aggressive invader is problematic throughout the Gulf Coastal Plain from Texas to Florida [72]. The tree has a wide means of dispersal, being primarily distributed through frugivory via avian populations. A high rate of germination compounded with quick growth and a sexual maturity of less than three years make it a difficult species to control. Its invasive capability is seen in nearly all forest types in the region, although it is more pervasive in coastal forests than inland. Modeling studies show a high likelihood of Chinese tallow becoming the dominant species within this ecosystem and in surrounding areas [73]. Future publications from this lab will provide more in-depth data and analysis on the state of existing stands and recruitment trends within remnants of the Chenier Plains via field sampling.

Louisiana has enacted a historic response to the problem of land-loss following a public outcry after devastating storms in 2005. The Coastal Protection and Restoration Authority (CPRA) have been tasked with enacting a management plan to combat coastal erosion. The Comprehensive Master Plan details proposals to create marshes, build sediment diversions, form new hydrologic structures, and promote barrier island restoration through built and natural means. As of 2017, over 135 projects had been completed impacting over 14,568 ha [74]. An update to this document is due to be released in 2023. The project has an overall operating budget of over 50 billion US dollars with equal shares dedicated to restoration and risk reduction. The primary funding source is a settlement associated with the Deepwater Horizon Oil Spill, in addition to various initiatives at the state and federal levels. However, there are no current plans to supplement vegetation or perform any geoengineering work on the cheniers for afforestation.

Coastal wetlands rank among the most important ecosystems in terms of the valuation of ecosystem services [75]. Provisioning services such as water filtration or storm-surge protection are commonly considered when discussing coastal areas. Beyond these commonalities, chenier forests provide vital wildlife habitats being the first encountered source of vertical structure for migratory neo-tropical birds and migrating monarch butterflies in the Mississippi alluvial valley flyway [76]. This en-route habitat is of critical importance due to its closeness to the Gulf of Mexico. During favorable weather conditions, 10% of birds migrating through the Gulf use these forests as a stopover spot, rising dramatically to 80% during adverse weather conditions [29]. With losses mounting, many species may not be able to seek refuge during storms. The area also has a unique suite of amphibians and reptiles, being home to 13 obligate chenier forest species of herpetofauna [77].

Lu et al. [78] used the eddy covariance method (ECM) to calculate strong carbon sink effects for a chenier ecosystem driven by forest vegetation. Net energy production values were significantly higher than surrounding marshland, signifying the importance of dense vegetative cover on a per-unit basis. Counter to this, there is also strong evidence that degraded forest ecosystems act as net carbon sources. A study in North Carolina also using ECM found that a seasonally inundated tide-influenced forest was a net carbon source annually, and only functioned as a sink for 3 months of the growing year [79]. Similarly, a study in the northeastern part of the state showed net negative rates of carbon sequestration for an ephemerally flooded forest [80]. Unique environments demand individual attention, a characterization of this ecosystem for nitrogen or carbon flux with respect to a similar system may not be functionally sound [81].

5. Conclusions

The cheniers of southwestern Louisiana were found to be declining in size at a concerning rate. Many of the driving factors are ecosystem-wide trends that have been increasing in intensity within the last century. The results of this study suggest that coastal forests are an ecosystem in peril and the ecology of such systems should be a closely monitored subject. The 36.5% loss of forest cover is an alarming but not wholly unexpected result of many decades of intense perturbation and human exploitation. A dynamic of continued forest loss is evident. Unless significant investment is made in engineering solutions to mitigate storm effects and land subsidence, we will continue to see the decline of chenier forests in Louisiana. Through targeted mitigation, governmental programs such as the Wetland or Conservation Reserve Program may be able to supplement existing forests, as well as the incomes of local people who have lost their livelihoods. In this study, we have begun the process of accurate areal estimates of forest extent, and have shed light on ecosystem processes in the region. As this community continues to diminish, we see the forest, not by the trees, but rather the lack of them.