Spatiotemporal patterns of rising annual plant abundance in grasslands of the Willamette Valley, Oregon (USA)

Context: Plant communities are undergoing compositional changes that affect ecosystem function. These changes are not always uniform across the landscape due to heterogenous topographic and edaphic conditions. To predict areas most at risk of change, it is necessary to identify the landscape drivers affecting plant abundance.

Objectives: Annual plants are increasing across the Wwestern USA, largely driven by non-native annual invasions. Here, we quantified change in annual plant abundance and identified landscape factors contributing to that change over the past 35 years.

Methods: We focused on Willamette Valley (Oregon) grasslands because they represent a new example in this phenomenon. To understand the spatiotemporal patterns of annual plant abundances between 1986 and 2020, we combined a remote-sensing vegetation cover dataset from the rangeland analysis platform with gridded soils data and topographic variables. We determined the rate of change in percent cover for each 30 × 30 m pixel and regressed cover against heat load, soil depth and sand content for > 5975 hectares to determine areas most sensitive to rising annual cover.

Results: We found a tendency toward increasing annual cover, with a median gain of + 15% cover among significantly increasing pixels. However, change was uneven across the landscape, with annual cover increasing markedly in areas with high heat load and shallower soils.

Conclusions: We identified steep, south-facing slopes as being particularly sensitive to rising annual cover. Annual plant invasions may be lagging in this region compared to elsewhere in the Wwestern USA, but trends here suggest it may just be a matter of time.

To perform this study, this dataset provides all necessary variables (as .csv tables), as well as a shapefile for the study area:

"df_RAP.csv": Annual forb and grass (AFG) and perennial forb and grass (PFG) cover data collected directly from the Rangeland Analysis Platform public datasets via Google Earth Engine. Total (herbaceous) cover was calculated by adding AFG and PFG cover.
"df_ratechange.csv": Rate of change data calculated from linear regressions of AFG cover against year from the "df_RAP.csv" dataset.
"df_xvars.csv": Soil and topographic variables collected/calculated directly from public platforms.
"sites16shape": Zipped folder containing shapefiles for the 16 study sites of interest.

For a detailed description of the datasets, please refer to the README file.

Methods
df_RAP.csv: Yearly (1986-2020) vegetation cover data of annual and perennial forbs and grasses for the Willamette Valley were downloaded as TIFF files (30-m resolution) from the Rangeland Analysis Platform's Google Earth Engine catalog. To eliminate forested and other non-grassland areas, we masked RAP data to grassland/herbaceous, pasture/hay, and shrub/scrub classifications in the National Land Cover Database 2019 Land Cover dataset using the raster package in R. Prior to this, we had previously resampled the land cover dataset from 10-m to 30-m resolution and reprojected to WGS 1984 using ArcMap 10.5. For each year in the 1986-2020 RAP dataset, we calculated new rasters for total herbaceous cover by summing the annual and perennial layers. We then filtered to cells with ≥20% average total herbaceous cover across the 35 years to avoid areas which may have skewed estimates due to low herbaceous cover in general (e.g., grassland borders near forests). Finally, we filtered to the 16 sites of interest using shapefiles of site polygon boundaries.

df_ratechange.csv: We used our processed RAP data in R to conduct pixel-by-pixel linear regressions of annual percent cover against year. We then generated a new raster layer for the rate of change using the coefficients from these models. 

df_xvars.csv: For soil data (depth, sand, silt, and clay), we obtained 2019 gSSURGO data (10-m resolution) for the state of Oregon through the USDA Natural Resources Conservation Service web portal (https://gdg.sc.egov.usda.gov/). Within this geodatabase, each soil map unit corresponds to one or more soil polygons represented in raster form. Each map unit contains several soil component records, which each have a percent attribute indicating their relative proportional composition within the map unit. A component record may have several associated ‘chorizon’ records, which each indicate their top and bottom depths. For soil depth, we used the ‘brockdepmin’ variable in the ‘muaggatt’ (map unit-aggregated) table when available. For map unit records with no data for ‘brockdepmin’, we calculated a minimum depth by choosing the deepest non-rock horizon of each component (i.e., the max ‘hzdept_r’ from the ‘chorizon’ table), and then averaging these numbers weighted by the component percent. For percent sand, silt, and clay, we used the attribute for the topmost horizon in the dominant components of the mapunit (e.g., ‘sandtotal_r’ from the ‘chorizon’ table where ‘hzdept_r’ = 0). We then converted these layers to WGS 1984 and resampled from 10-m to 30-m resolution using ArcMap 10.5. For topographic variables (elevation, slope, aspect, and heatload), we downloaded a 30-m resolution digital elevation model file from OpenTopography (https://portal.opentopography.org/raster?opentopoID=OTSRTM.082015.4326.1). To calculate slope and aspect, we used the Slope and Aspect tools from the Spatial Analyst toolkit in ArcMap 10.5. Heat load is a unitless index of the amount of heat a ground surface receives through solar radiation based on its slope, aspect, and latitude. We calculated heat load in R using the McCune and Keon (2002) method with south-facing slopes (180°) at a maximum. Once we had calculated raster layers for all these soils and topographic variables, we masked them to the RAP raster data layers.