Sentinel Hub & comparing NDVI over time: an exploration of a ghost forest

The New Jersey Pinelands are an expansive (1.1 million acres) and biodiverse Atlantic coastal pine barrens ecosystem that sits atop an estimated 17 million gallons of drinking water. Like many Atlantic coastal areas in the eastern United States, the Pinelands are becoming increasingly affected by climate change, specifically through sea level rise and the intrusion of saltwater, which has to led to the emergence of "ghost forests." Ghost forests are characterized by dead stands of trees -- usually species that can normally withstand wet soil, such as pines and cedars -- in coastal regions. The emergence of ghost forests is typically gradual (over the course of multiple decades), but it is a clear sign of climate change.

I spend a lot of time exploring the Pinelands, and in recent years, I've started to notice the presence of ghost forests. I know from existing research on ghost forests (Cheriegate et al., 2023 and Martinez et al., 2024 to name only a couple) that remote sensing techniques are valuable in analyzing ghost forests in the eastern United States. NDVI, which measures the health of vegetation, is adept at identifying transformations from forest to marsh and therefore may be particularly well-suited for analyzing the progression of ghost forests in this exploration.

Because the processes driving ghost forests evolve over multiple decades, many researchers use Landsat imagery (which dates back to the 1970s) to analyze ghost forests over long periods of time. Although Landsat imagery may be well-suited for this exploration, I've been meaning to explore Sentinel Hub and Sentinel-2 imagery (due to its high spatial resolution and frequent captures), so I decided to apply it in this quick exploration. My interest in using SentinelHub is rooted in its use of cloud infrastructure, which lessens the burden of this work on my computer, and its easy handling of pre-processing steps.

# for accessing and processing sentinel imagery
import datetime
import os
import getpass
from sentinelhub import (
    SHConfig,
    CRS,
    BBox,
    DataCollection,
    DownloadRequest,
    MimeType,
    MosaickingOrder,
    SentinelHubDownloadClient,
    SentinelHubRequest,
    bbox_to_dimensions,
)
from oauthlib.oauth2 import BackendApplicationClient
from requests_oauthlib import OAuth2Session
from PIL import Image
# for plotting
import matplotlib.pyplot as plt
# for linear regression
import numpy as np
# for masking rasters
from rasterio.transform import from_bounds
from rasterio.features import geometry_mask
# for downloading polygon as geopandas geodataframe
import io
from io import StringIO
import geopandas as gpd
import requests


# The following is not a package. It is a file utils.py which should be in the same folder as this notebook.
from utils import plot_image

/Users/evelyngorey/Documents/projects/ghostforests/ghostforest-rs/env/lib/python3.9/site-packages/urllib3/__init__.py:35: NotOpenSSLWarning: urllib3 v2 only supports OpenSSL 1.1.1+, currently the 'ssl' module is compiled with 'LibreSSL 2.8.3'. See: https://github.com/urllib3/urllib3/issues/3020
  warnings.warn(

1. Get Sentinel-2 imagery for 2016 and 2024

1.1 Set up configuration

The code below initializes the Sentinel Hub configuration and checks if OAuth credentials are set. For this to work, I created a new OAuth client from my Sentinel Hub dashboard and copied my unique client ID and secret into my config.toml file.

# get config w/ Sentinel Hub client ID and secret
config = SHConfig()

if not config.sh_client_id or not config.sh_client_secret:
    print("Warning! To use Process API, please provide the credentials (OAuth client ID and client secret).")

1.2 Set study area

The code below sets the specific area of interest, which I found using bbox finder, as well as the resolution (10m) and CRS (WGS84) of my image.

This area, which is only about 6.8 square miles, surrounds Batsto Village (coincidentally a ghost town) in Wharton State Forest. The area is relatively small, as I'm just testing the methods below.

pb_coords = (-74.682857, 39.635805, -74.617969, 39.6641237)

# initialize bbox of interest 
resolution = 10
pb_bbox = BBox(bbox=pb_coords, crs=CRS.WGS84)
pb_size = bbox_to_dimensions(pb_bbox, resolution=resolution)

print(f"Image shape at {resolution} m resolution: {pb_size} pixels")

Image shape at 10 m resolution: (555, 316) pixels

1.3 Explore true color mosaic for one year

1.3.1 Evalscript

The evalscript below (taken from the Sentinel Hub Process API example notebook) specifies that when we query our mosaic, we want just the true color bands (4, 3, and 2). This is helpful for quickly visualizing the study area.

evalscript_true_color = """
    //VERSION=3

    function setup() {
        return {
            input: [{
                bands: ["B02", "B03", "B04"]
            }],
            output: {
                bands: 3
            }
        };
    }

    function evaluatePixel(sample) {
        return [sample.B04, sample.B03, sample.B02];
    }
"""

1.3.2 Get and plot image

In the code block below, we use the evalscript from above to query a mosaic comprising the least cloudy (leastCC) acquisitions within our timeframe of interest (January 1st, 2025 through March 1st, 2025). The most recent (mostRecent) image(s) from all available images are used in the mosaicking order by default.

request_true_color = SentinelHubRequest(
    evalscript=evalscript_true_color,
    input_data=[
        SentinelHubRequest.input_data(
            data_collection=DataCollection.SENTINEL2_L1C,
            time_interval=("2025-01-01", "2025-03-01"),
            mosaicking_order=MosaickingOrder.LEAST_CC,
        )
    ],
    responses=[SentinelHubRequest.output_response("default", MimeType.PNG)],
    bbox=pb_bbox,
    size=pb_size,
    config=config,
)

plot_image(request_true_color.get_data()[0], factor=3.5 / 255, clip_range=(0, 1))

1.4 Get raw band mosaics for multiple years, calculating NDVI

1.4.1 Evalscript

The evalscript below defines our data request from Sentinel Hub. We specify that we want all Sentinel-2 spectral bands with raw values. For a quick download, we want raw numbers (as opposed to reflectances) downloaded in INT16 format, so the input units are set to DN (digital numbers) and the output image type is a TIFF. We also specify the NDVI calculation and use this to get an NDVI band for each year. As mentioned before, NDVI has proven to be especially helpful when distinguishing ghost forests.

evalscript_all_bands = """
    //VERSION=3
    function setup() {
        return {
            input: [{
                bands: ["B01","B02","B03","B04","B05","B06","B07","B08","B8A","B09","B10","B11","B12"],
                units: "DN"
            }],
            output: {
                bands: 14,
                sampleType: "FLOAT32"
            }
        };
    }

    function evaluatePixel(sample) {
        var ndvi = (sample.B08 - sample.B04) / (sample.B08 + sample.B04);  // NDVI calculation

        return [sample.B01,
                sample.B02,
                sample.B03,
                sample.B04,
                sample.B05,
                sample.B06,
                sample.B07,
                sample.B08,
                sample.B8A,
                sample.B09,
                sample.B10,
                sample.B11,
                sample.B12,
                ndvi];
    }
"""

1.4.2 Get imagery for every year from 2016-2025

We use the evalscript defined above and the SentinelHubRequest class to return all Sentinel-2 spectral bands with raw values in our study area and time period, which is January 1st, 2016 through March 1st, 2025. I chose these years because Sentinel-2 imagery becomes available in 2016 and 2025 is the most recent year. I chose January through March to ensure that I'm consistently looking at the same season for each year, so I know that any drastic changes in NDVI likely aren't attributed to typical seasonal vegetation variation. Also, this can help focus more on evergreen forests, as deciduous foluage is dormant during this time...

# dictionary to store results for each year
all_band_data = {} 

# define loop - getting imagery for every year from 2016-2025
for year in range(2016, 2026):  
    request = SentinelHubRequest(

        # save imagery in respective folders
        data_folder=f"s2_{year}_dir",
        # use the evalscript we defined above
        evalscript=evalscript_all_bands,

        input_data=[
            SentinelHubRequest.input_data(
                # define sensor (Sentinel)
                data_collection=DataCollection.SENTINEL2_L1C,
                # same time frame every year (January 1st through April 1st for consistency)
                time_interval=(f"{year}-01-01", f"{year}-04-01"),
                # choose images with least amout of cloud cover
                mosaicking_order=MosaickingOrder.LEAST_CC,
            )
        ],
        # imagery in TIFF format
        responses=[SentinelHubRequest.output_response("default", MimeType.TIFF)],
        # set the area of interest
        bbox=pb_bbox,
        size=pb_size,
        config=config,
    )

    # download and store data
    data = request.get_data(save_data=True)
    all_band_data[year] = data
    print(f"Downloaded data for {year}")

Downloaded data for 2016
Downloaded data for 2017
Downloaded data for 2018
Downloaded data for 2019
Downloaded data for 2020
Downloaded data for 2021
Downloaded data for 2022
Downloaded data for 2023
Downloaded data for 2024
Downloaded data for 2025

1.4.3 Plot NDVI for each year

Because NDVI is such a valuable predictor of ghost forests, we can plot this band for each year to see if there are any noticeable changes that we can spot. NDVI is always supposed to be a value between 1 (healthy vegetation) and 0. We can easily distinguish areas with denser, healthier vegetation (darker green in the images, or a value closer to 1) versus water (white, or a value of 0) and less dense vegetation like grass (lighter green, or a value somewhere in between 0 and 1).

# set subplots according to number of years
n_years = len(all_band_data)
ncols = 5
nrows = (n_years + ncols - 1) // ncols

fig, axes = plt.subplots(nrows=nrows, ncols=ncols, figsize=(20, 4 * nrows))
axes = axes.flatten()

# plot each year's NDVI band
for i, (year, data) in enumerate(sorted(all_band_data.items())):
    img = data[0]  # get the single image for each year
    # handle different shape cases
    if img.shape[0] == 14:  # shape: (bands, H, W)
        ndvi = img[13, :, :]
    elif img.shape[-1] == 14:  # shape: (H, W, bands)
        ndvi = img[:, :, 13]

    ax = axes[i]
    im = ax.imshow(ndvi, cmap='YlGn', vmin=0, vmax=1)
    ax.set_title(f"NDVI: {year}")
    ax.axis('off')

# hide unused axes
for j in range(i + 1, len(axes)):
    axes[j].axis('off')

# add colorbar and legend to the right
fig.subplots_adjust(right=0.92)
cbar_ax = fig.add_axes([0.94, 0.15, 0.015, 0.7])
fig.colorbar(im, cax=cbar_ax, label='NDVI')

plt.tight_layout(rect=[0, 0, 0.93, 1])
plt.show()

/var/folders/13/3rdcyfxx2v72vql4x3d2v9j40000gn/T/ipykernel_86669/4025519471.py:32: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.
  plt.tight_layout(rect=[0, 0, 0.93, 1])

2. Mask to only include Freshwater Forested/Shrub Wetland land cover

2.1 Import National Wetlands Inventory wetlands data

Since we'll be comparing NDVI by year, it would be helpful to focus only on the type of land cover we know ghost forests will likely be in. I know that low-lying forests near marshes or open water make up the type of ecosystem that could turn into ghost forests. Elevation data, such as a DEM or LiDAR, might be helpful here, but to keep it simple and use what's easily accessible, I'll use the National Wetlands Inventory to help me focus on ghost forest-prone areas in particular. The National Wetlands Inventory identifies multiple different wetlands classes in the Pinelands, including areas that have been classified as "freshwater forested/shrub wetland." I'll use this class to narrow my focus.

In the code below, I download the National Wetlands Inventory data (filtered to only include "freshwater forested/shrub wetland" polygons) for my study area via the Wetlands Map Service API and read them in as a geopandas GeoDataFrame. I visualize the polygons in a dynamic Leaflet map using the geopandas explore() method.

You can see that the polygons do not include agricultural areas, open fields, or urban features such as streets and buildings. Water bodies and dense forests that are not considered freshwater forested or shrub wetland are also excluded.

# define the layer URL
service_url = "https://fwspublicservices.wim.usgs.gov/wetlandsmapservice/rest/services/Wetlands/MapServer/0/query"

# build the query
params = {
    "where": "1=1",
    "geometry": ",".join(map(str, pb_bbox)),
    "geometryType": "esriGeometryEnvelope",
    "inSR": "4326",
    "spatialRel": "esriSpatialRelIntersects",
    "outFields": "*",
    "outSR": "4326",
    "f": "geojson"
}

# make the request
response = requests.get(service_url, params=params)

# read in as geodataframe with geopandas
geojson_data = response.text
gdf = gpd.read_file(StringIO(geojson_data))
# filter to only include freshwater forested/shrub wetland
gdf = gdf[gdf["Wetlands.WETLAND_TYPE"] == "Freshwater Forested/Shrub Wetland"]

# interactive map to explore layer
gdf.explore(tiles="Esri.WorldImagery")

Make this Notebook Trusted to load map: File -> Trust Notebook

2.2 Mask yearly imagery with wetlands data

The function mask_wetlands() below uses the rasterio package to mask (or clip) my imagery that contains mosaics for all years by the GeoDataFrame that I loaded above. The result is a raster that only contains the areas that are considered to be Freshwater Forested/Shrub Wetland -- open water, dense non-freshwater forests, urban features like buildings and roads, agriculture, and all other land covers will not be included in the output. This ensures that we're only comparing NDVI across the land covers that are most likely to contain ghost forests, and we're minimizing any noise that would come from comparing the NDVI of other land covers.

def mask_wetlands(all_band_data, gdf, pb_bbox, pb_size, crs="EPSG:4326"):
    # map raster pixel coordinates to real-world coordinates
    width, height = pb_size
    minx, miny, maxx, maxy = pb_bbox
    transform = from_bounds(minx, miny, maxx, maxy, width, height)

    # reproject geodataframe to the same CRS as the raster
    gdf = gdf.to_crs(crs)

    # rasterizes polygons and masks -- if pixels are inside polygon -> true, else -> false
    mask = geometry_mask(
        gdf.geometry, out_shape=(height, width), transform=transform, invert=True
    )

    # empty dictionary to store the masked arrays for each year
    masked_data = {}
    # loops over each year's raster
    for year, data_list in all_band_data.items():
        # get the first item and convert to numpy array to store in dictionary
        arr = np.asarray(data_list[0])

        # standardize array shape (bands, height, width)
        if arr.ndim == 2:  # single band (height, width)
            arr = arr[np.newaxis, :, :]
        elif arr.ndim == 3 and arr.shape[-1] <= 20:  
            # assume height, width, bands
            arr = arr.transpose(2, 0, 1)

        # apply the mask -- where true, keep pixel value; where false, replace with NaN
        masked_arr = np.where(mask[np.newaxis, :, :], arr, np.nan)
        # store results
        masked_data[year] = masked_arr

    return masked_data

# run function
masked_data = mask_wetlands(all_band_data, gdf, pb_bbox, pb_size)

2.3 Plot masked yearly NDVI

We can plot the results to see what of our original yearly images was kept. Again, we see that open water, urban features such as roads and buildings, agriculture/fields, and other dense forest are excluded now.

def plot_ndvi_years(masked_data, years=None):
    # sort years chronologically
    if years is None:
        years = sorted(masked_data.keys())
    
    # adjust columns/subplots according to the number of years
    n = len(years)
    ncols = 3
    nrows = int(np.ceil(n / ncols))
    
    # creates grid of subplots (all 5x5)
    fig, axes = plt.subplots(nrows, ncols, figsize=(5*ncols, 5*nrows))
    axes = axes.flatten() # easier to loop through
    
    # loop through each year
    for ax, year in zip(axes, years):
        # gets the last band of each year's raster (NDVI)
        ndvi = masked_data[year][-1] 
        # assigns color scheme and fixed sclae for NDVI
        im = ax.imshow(ndvi, cmap="RdYlGn", vmin=-1, vmax=1)
        ax.set_title(str(year))
        ax.axis("off")
    
    # remove extra subplots if any
    for ax in axes[len(years):]:
        ax.axis("off")
    
    # adjust layout and colorbar/legend
    plt.tight_layout()
    plt.colorbar(im, ax=axes, orientation="vertical", fraction=0.02, pad=0.01)
    plt.show()

plot_ndvi_years(masked_data)

3. Compare NDVI by year

3.1 Calculate mean and median NDVI per year

One very simple way to attempt to understand how ghost forests are emerging/changing over time is by seeing how the mean and median NDVI changes across years. In the code block below, the mean and median NDVI for each year is calculated.

# empty dictionary to store the NDVI values per year
ndvi_summary = {}

# loop through each year
for year, arr in masked_data.items():
    # get the last band only, which is NDVI
    ndvi = arr[-1] 
    # only include non-NaN pixels
    valid_ndvi = ndvi[~np.isnan(ndvi)]
    
    # specify which summary statistics we want for NDVI
    ndvi_summary[year] = {
        "mean": np.mean(valid_ndvi),
        "median": np.median(valid_ndvi),
        "std": np.std(valid_ndvi),
        "min": np.min(valid_ndvi),
        "max": np.max(valid_ndvi)
    }

# print results
for year, stats in ndvi_summary.items():
    print(year, stats)

2016 {'mean': np.float32(0.48649815), 'median': np.float32(0.50368327), 'std': np.float32(0.07916312), 'min': np.float32(-0.16950673), 'max': np.float32(0.66997874)}
2017 {'mean': np.float32(0.5100661), 'median': np.float32(0.52590615), 'std': np.float32(0.07966204), 'min': np.float32(-0.0738255), 'max': np.float32(0.71281296)}
2018 {'mean': np.float32(0.47524548), 'median': np.float32(0.5002192), 'std': np.float32(0.108640105), 'min': np.float32(-0.114922814), 'max': np.float32(0.6862031)}
2019 {'mean': np.float32(0.5058639), 'median': np.float32(0.5300601), 'std': np.float32(0.09281689), 'min': np.float32(-0.10233393), 'max': np.float32(0.70484245)}
2020 {'mean': np.float32(0.49970126), 'median': np.float32(0.5205171), 'std': np.float32(0.085757256), 'min': np.float32(-0.08712121), 'max': np.float32(0.6887097)}
2021 {'mean': np.float32(0.50415915), 'median': np.float32(0.5241972), 'std': np.float32(0.08876265), 'min': np.float32(-0.13504274), 'max': np.float32(0.6918667)}
2022 {'mean': np.float32(0.5028465), 'median': np.float32(0.52116585), 'std': np.float32(0.08588857), 'min': np.float32(-0.024822695), 'max': np.float32(0.69151056)}
2023 {'mean': np.float32(0.46895048), 'median': np.float32(0.48846745), 'std': np.float32(0.10058494), 'min': np.float32(-0.07722772), 'max': np.float32(0.6834593)}
2024 {'mean': np.float32(0.4842454), 'median': np.float32(0.49968925), 'std': np.float32(0.09759679), 'min': np.float32(-0.08988764), 'max': np.float32(0.69106346)}
2025 {'mean': np.float32(0.4841035), 'median': np.float32(0.49714267), 'std': np.float32(0.08633356), 'min': np.float32(-0.20008285), 'max': np.float32(0.6790378)}

3.2 Plot NDVI mean and median by year

The mean and median NDVI for each year are plotted below. The results look somewhat meaningful -- we can see that the mean and median have very similar patterns throughout the years, with the median always being a bit higher than the mean; this indicates that there is some variation that the median doesn't necessarily pick up. 2019 and 2023 have the lowest values for both mean and median of all the years.

From our NDVI plots, we see lower NDVI values among vegetated areas in 2019 than in other years. Could this be due to the March 2019 Spring Hill Wildfire and/or drought conditions leading up to it? I couldn't find a particular wildfire that could have caused 2023's dip in NDVI, but it is very likely that a fire could be contributed for that year as well.

Regardless of these two dips, we are seeing that the mean and median NDVI are trending downward over time. Whether this is due to the emergence of ghost forests is a different question.

# sort years in chronological order
years = sorted(ndvi_summary.keys())
# define mean and median only
means = [ndvi_summary[y]["mean"] for y in years]
medians = [ndvi_summary[y]["median"] for y in years]

plt.figure(figsize=(8,5))
# plot mean and median only
plt.plot(years, means, marker="o", label="Mean NDVI")
plt.plot(years, medians, marker="x", label="Median NDVI")
plt.xlabel("Year")
plt.ylabel("NDVI")
plt.title("NDVI trend across years")
plt.legend()
plt.grid(True)
plt.show()

3.3 Linear regression to identify changed pixels

We see that the NDVI among the freshwater forests/shrubs in our study area is decreasing over time, but it would be helpful to know where (which pixels) this decline is specifically occurring. A linear regression allows us to determine whether the NDVI is increasing or decreasing by year for each pixel and at what rate.

# order the years chronologically
years = sorted(masked_data.keys())
# build a stack of NDVI rasters
ndvi_stack = np.array([masked_data[y][-1] for y in years])  # (time, row, col)
# each number represents one year in the sequence
x = np.arange(len(years))

# calculate the difference between 2016 (first year) and 2025 (most recent year) -- produces a map of change
ndvi_diff = ndvi_stack[-1] - ndvi_stack[0] # positive = NDVI increased, negative = NDVI decreased

# flatten raster stack (makes regression easier)
n_time, rows, cols = ndvi_stack.shape
pixels = rows * cols
ndvi_2d = ndvi_stack.reshape(n_time, pixels)

# apply regression using polyfit
slopes = np.full(pixels, np.nan)
for p in range(pixels):
    y = ndvi_2d[:, p]
    # skip pixels w/ NaN values for all years
    if np.all(np.isnan(y)):
        continue
    # get slope -> rate of NDVI change per year
    coeffs = np.polyfit(x, y, 1) 
    slopes[p] = coeffs[0]

# reshape slopes back to a raster for mapping pixel change
ndvi_trend = slopes.reshape(rows, cols)

Now we can plot the different types of change we calculated. The first is simply looking at the difference between the NDVI values of the first year of data (2016) and the most recent year (2025). This isn't really as helpful as seeing the change over every year, but it can give a quick understanding of what's changing and where. Note that we see very little change overall; we see some areas experiencing a decline in NDVI, but very few areas experienced an increase. Visualizing the rate of change (below) across all years may be more meaningful.

plt.figure(figsize=(8,6))
plt.imshow(ndvi_diff, cmap="RdYlGn", vmin=-1, vmax=1)
plt.colorbar(label="NDVI change")
plt.title(f"NDVI change from {years[0]} to {years[-1]}")
plt.axis("off")
plt.show()

We can also visualize pixel-wise NDVI direction and rate of change per year, across all years (2016-2025). NDVI change values are between -1 in the red (meaning NDVI is steadily decreasing across years), 0 in the yellow (no change), and 1 in the green (NDVI is steadily increasing).

It's difficult to tell without a basemap and context of the area, but the areas with the steadily decreasing NDVI (the red pixels) are along the Mullica River, which other sources (e.g, here and here) have confirmed is a hostpot of ghost forest emergence. We also see a bit of red next to Batsto Lake (the main water body in these images) and in the Sleeper Branch area to the north of the agricultural activity. Whether these are statistically significant changes is another question?

Though we still cannot confirm that the decrease in NDVI is due to the emergence of ghost forests, this exploration of the forests' health over time is a start.

plt.figure(figsize=(8,6))
plt.imshow(ndvi_trend, cmap="RdYlGn")
plt.colorbar(label="NDVI change")
plt.title("Pixel-wise NDVI trend across 2016-2025")
plt.axis("off")
plt.show()

Closing thoughts

This was just a quick, simplistic exploration of the applications of Sentinel Hub and Sentinel-2 imagery. Although I may not pay for a personal Sentinel Hub subscription, I did enjoy the ease of the pre-processing and mosaicking steps, and the high spatial and spectral resolution of the imagery was useful for such a small study area.

If I were to expand on this analysis, I would use a larger study area and study the change over at least 15 years, which would entail switching to Landsat or other multispectral imagery that is available for a longer timeframe. I would also go beyond simply comparing NDVI over time by training a supervised classification model to identify perhaps multiple stages of ghost forests over time (e.g., emerging, longstanding). Because ghost forests usually occur among species that are tolerant of wet soils, such as pines and cedars, the ability to distinguish between tree species and focus just on the species of greatest concern could also be helpful. There are also, of course, dozens of other data sources, processes, and considerations that could be added to this analysis to make it far more meaningful.