Skip to main content

Extracting Time Series from Forest Carbon Diligence

This Jupyter Notebook demonstrates how to derive time-series of statistics from the Forest Carbon Diligence product using a combination of the Statistical API and Processing API.

Statistics (more details below) will be produced for the different variables contained in the Forest Carbon Diligence product:

  • Canopy Height: quantifies the average height of all vegetation per pixel;
  • Canopy Cover: quantifies the horizontal area occupied by tree canopies that are > 5 m tall;
  • Aboveground carbon density: quantifies the density of carbon stored in woody vegetation, primarily trees and shrubs;
  • Biomass: converted from aboveground carbon density.

The examples produced in this Jupyter Notebook are divided into four steps.

  1. After plotting the parcel of interest (imported from a GeoJSON file) on a map, we will extract the following statistics using the Statistical API:
  • Minimum
  • Maximum
  • Average
  • Standard deviation
  • Data/No Data count

  1. The Processing API extracts the sum of aboveground carbon density and biomass over the parcel of interest.

  2. Query the same statistics as in step 1, but masking non-forested areas out. The forest mask is derived from a combination of height and cover thresholds: for the purpose of this demonstration, forested areas are considered as pixel with Canopy Height > 5 m and Canopy Cover > 15% (note that these thresholds can be changed in the script).

  3. The Processing API is used to plot the Canopy Height over the parcel of interest at two different dates. The dates will be chosen based on the time series previously plotted, highlighting how a combination of statistics and visualisation of products can support the analysis of various phenomena.

import getpass
import json
import geopandas as gpd
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Patch
import matplotlib.colors as mcolors
from shapely.geometry import Polygon
from shapely.ops import transform

from sentinelhub import (
    SHConfig,
    geometry,
    CRS,
    SentinelHubStatistical,
    DataCollection,
    MimeType,
    SentinelHubRequest,
    BBox
    )

Credentials​

To obtain your client_id & client_secret, you need to navigate to your Dashboard. In the User Settings, you can create a new OAuth client to generate these credentials. More detailed instructions can be found on the corresponding documentation page.

Now that you have your client_id & client_secret, it is recommended to configure a new profile in your Sentinel Hub Python package. Instructions on how to configure your Sentinel Hub Python package can be found here. This is useful as changes to the config class in your notebook are usually only temporary and by saving the configuration to your profile, you don't have to generate new credentials or overwrite/change the default profile every time you run or write a new Jupyter Notebook.

Currently, the following cell is set up in a way that you need to specify your client_id & client_secret when running it.

config = SHConfig()

if not config.sh_client_id or not config.sh_client_secret:
    print("No credentials found, please provide the OAuth client ID and secret.")
    config.sh_client_id = getpass.getpass("sh_client_id: ")
    config.sh_client_secret = getpass.getpass("sh_client_secret: ")
    config.save()
    print(f"Credentials saved to {SHConfig.get_config_location()}")
else:
    print(f"Using credentials stored here: {SHConfig.get_config_location()}")
Output
    Using credentials stored here: C:\Users\matt.ballard\.config\sentinelhub\config.toml

Set collection IDs​

The Forest Carbon Diligence layers are available as subscriptions through Sentinel Hub's TPDI service. Once the area of interest and variables are subscribed to, the data are then automatically imported into Sentinel Hub and available by specifying the Collection ID.

In this example, we will use the collection IDs available from Planet's Sandbox Data page for Forest Carbon Diligence.

A collection ID will allow you to retrieve the data from the collection, just as you would with a standard dataset (e.g. Sentinel-2). For more information on how to call a collection ID in a request with Python, you can refer to the sentinelhub-py documentation.

# Replace with your own collection IDs
height_collection = DataCollection.define_byoc("f3312c82-edea-42a1-8c9d-ada86ddcc857")

cover_collection = DataCollection.define_byoc("e3d2a21c-cb75-4311-86ac-024385c85b9c")

carbon_collection = DataCollection.define_byoc("cc31cada-80d8-46fe-a746-43ac2f87b5da")

Area of Interest​

First, we define an Area of Interest (AOI), located close to São Félix do Xingu, in Brazil.

In this example, the AOI is defined in a GeoJSON file in EPSG:3857. We will read the geoJSON and plot the area of interest on an interactive map.

You can also explore the area of interest in Sentinel Hub's EO Browser.

# Read a geojson containing a polygon representing a parcel of interest
parcel_data = json.loads(open("./aoi.geojson").read())
parcel_polygon = Polygon(parcel_data["features"][0]["geometry"]["coordinates"][0])

# Convert to a Sentinel Hub geometry
parcel_geo = geometry.Geometry(parcel_polygon, crs=CRS(3857))

del parcel_data
evalscript = """
//VERSION=3
const defaultVis = true;
const max = 100;
const min = 0;

function setup() {
  return {
    input: ["CC", "dataMask"],
    output: { bands: 4, sampleTYPE: "AUTO" },
  };
}

function updateMap(max, min) {
  const numIntervals = map.length;
  const intervalLength = (max - min) / (numIntervals - 1);
  for (let i = 0; i < numIntervals; i++) {
    map[i][0] = max - intervalLength * i;
  }
}

const map = [
  [100, 0x183d19],
  [90, 0x124f24],
  [80, 0x0e6327],
  [70, 0x246d29],
  [60, 0x498418],
  [50, 0x669516],
  [40, 0x859e25],
  [30, 0xa4ab38],
  [20, 0xd3c058],
  [10, 0xddd17c],
  [0, 0xf0f5d5],
];

const visualizer = new ColorRampVisualizer(map, min, max);

function evaluatePixel(sample) {
  let val = sample.CC;
  let imgVals = visualizer.process(val);

  return [...imgVals, sample.dataMask];
}
"""

# Create a GeoDataFrame for plotting
gdf = gpd.GeoDataFrame({"geometry": [parcel_polygon]}, crs="EPSG:3857")
minx, miny, maxx, maxy = gdf.total_bounds
bbox = BBox((minx, miny, maxx, maxy), crs=CRS.POP_WEB)

request = SentinelHubRequest(
    evalscript=evalscript,
    input_data=[
        SentinelHubRequest.input_data(
            cover_collection,
            time_interval=('2016-10-01', '2017-10-31')
        )
    ],
    responses=[SentinelHubRequest.output_response('default', MimeType.JPG)],
    bbox=bbox,
    resolution=(30, 30),
    config=config
)

image = request.get_data()[0]


# Define the color map
map = [
    [100, 0x183d19],
    [90, 0x124f24],
    [80, 0x0e6327],
    [70, 0x246d29],
    [60, 0x498418],
    [50, 0x669516],
    [40, 0x859e25],
    [30, 0xa4ab38],
    [20, 0xd3c058],
    [10, 0xddd17c],
    [0, 0xf0f5d5],
]

# Convert hex colors to RGB
legend_elements = [
    Patch(
        facecolor=mcolors.to_rgb(f"#{color:06x}"),
        edgecolor="black",
        label=f"{value}%",
    )
    for value, color in map
]

# Plot the image and the GeoDataFrame
fig, ax = plt.subplots(figsize=(8, 6))

ax.imshow(image, extent=(minx, maxx, miny, maxy), origin="upper")

gdf.plot(
    ax=ax,
    color="none",
    edgecolor="black",
    linewidth=2,
)

# Add the legend
ax.legend(
    handles=legend_elements,
    title="Forest Cover (%)",
    loc="upper right",
    fontsize=8,
    title_fontsize=10,
)

ax.set_aspect("equal")
ax.set_xticks([])
ax.set_yticks([])
plt.tight_layout()
ax.set_title("Area of Interest with Forest Cover", fontsize=16)

plt.show()

Extract Statistics​

In the following cells we will extract statistics using the Statistical API and parse the results into a Dataframe, making it easier to analyse and plot the results. (Tip: learn more about the API in our webinar).

Results of Statistical API are aggregated statistical values of satellite data instead of entire images. By using Statistical API we can avoid downloading and processing large amounts of satellite data. All general rules for building Evalscripts apply. However, there are some specifics when using evalscripts with the Statistical API:

  • The evaluatePixel() function must, in addition to other output, always return a dataMask output. This output defines which pixels are excluded from calculations. For more details and an example, see here.
  • The default value of sampleType is FLOAT32.
  • The output.bands parameter in the setup() function can be an array. This makes it possible to specify custom names for the output bands and different output dataMask for different outputs, see this example.
evalscript_statistics = """
//VERSION=3

const BIOMASS_CONV = 0.476;

function setup() {
    return {
        input: [{
            datasource: "canopyHeight",
            bands: ["CH", "dataMask"]
        },
        {
            datasource: "canopyCover",
            bands: ["CC", "dataMask"]
        },
        {
            datasource: "carbon",
            bands: ["ACD", "dataMask"]
        }
        ],
        output: [
            { id: "canopyHeight", bands: 1, sampleType: "UINT8" },
            { id: "canopyCover", bands: 1, sampleType: "UINT8" },
            { id: "carbon", bands: 1, sampleType: "UINT16" },
            { id: "biomass", bands: 1, sampleType: "FLOAT32" },
            { id: "dataMask", bands: 1, sampleType: "UINT8" }
        ],
        mosaicking: "ORBIT"
    };
}

function evaluatePixel(samples) {
    var sampleHeight = samples.canopyHeight[0].CH
    var nodataHeight = samples.canopyHeight[0].dataMask
    var sampleCover = samples.canopyCover[0].CC
    var nodataCover = samples.canopyCover[0].dataMask
    var sampleCarbon = samples.carbon[0].ACD
    var nodataCarbon = samples.carbon[0].dataMask

    // Biomass simple scalar conversion
    var sampleBiomass = sampleCarbon / BIOMASS_CONV


    // Combined noData mask
    var allMask = nodataHeight * nodataCover * nodataCarbon

    return {
        canopyHeight: [sampleHeight],
        canopyCover: [sampleCover],
        carbon: [sampleCarbon],
        biomass: [sampleBiomass],
        dataMask: [allMask]
    };
}
"""
# Start year and end year
start_year = 2013
end_year = 2017

# Payload
request = SentinelHubStatistical(
    aggregation=SentinelHubStatistical.aggregation(
        evalscript=evalscript_statistics,
        time_interval=(f"{start_year}-01-01T00:00:00Z", f"{end_year}-12-31T23:59:59Z"),
        # To fetch the yearly data that is set to 1 date / year we set the aggregation to P1D
        aggregation_interval="P1D",
        # Match the resolution of the PV (30 meters). Beware: resolution should be in the units of the coordinate system used
        resolution=[30, 30],
    ),
    input_data=[
        SentinelHubStatistical.input_data(
            # Note that we set the service_url to the US-West-2 (Oregon) region since the data are stored there
            height_collection,
            identifier="canopyHeight",
        ),
        SentinelHubStatistical.input_data(
            cover_collection,
            identifier="canopyCover",
        ),
        SentinelHubStatistical.input_data(
            carbon_collection,
            identifier="carbon",
        ),
    ],
    geometry=parcel_geo,
    config=config,
)
# Run the request and display a sample of the data returned
statistical_response = request.get_data()
statistical_response
Output
    [{'data': [{'interval': {'from': '2013-01-01T00:00:00Z',
         'to': '2013-01-02T00:00:00Z'},
        'outputs': {'biomass': {'bands': {'B0': {'stats': {'min': 0.0,
             'max': 226.89076232910156,
             'mean': 94.24312644359381,
             'stDev': 54.250074546784866,
             'sampleCount': 51975,
             'noDataCount': 27551}}}},
         'carbon': {'bands': {'B0': {'stats': {'min': 0.0,
             'max': 108.0,
             'mean': 44.85972813625931,
             'stDev': 25.823035472491263,
             'sampleCount': 51975,
             'noDataCount': 27551}}}},
         'canopyHeight': {'bands': {'B0': {'stats': {'min': 0.0,
             'max': 27.0,
             'mean': 10.215075335735367,
             'stDev': 7.967104024080437,
             'sampleCount': 51975,
             'noDataCount': 27551}}}},
         'canopyCover': {'bands': {'B0': {'stats': {'min': 0.0,
             'max': 93.0,
             'mean': 51.64469374385916,
             'stDev': 30.24538432050859,
             'sampleCount': 51975,
             'noDataCount': 27551}}}}}},
       {'interval': {'from': '2014-01-01T00:00:00Z', 'to': '2014-01-02T00:00:00Z'},
        'outputs': {'biomass': {'bands': {'B0': {'stats': {'min': 0.0,
             'max': 224.7899169921875,
             'mean': 85.53647475003687,
             'stDev': 54.07485604925614,
             'sampleCount': 51975,
             'noDataCount': 27551}}}},
         'carbon': {'bands': {'B0': {'stats': {'min': 0.0,
             'max': 107.0,
             'mean': 40.71536193907608,
             'stDev': 25.739631454701577,
             'sampleCount': 51975,
             'noDataCount': 27551}}}},
         'canopyHeight': {'bands': {'B0': {'stats': {'min': 0.0,
             'max': 27.0,
             'mean': 8.683671798231284,
             'stDev': 8.108693034728725,
             'sampleCount': 51975,
             'noDataCount': 27551}}}},
         'canopyCover': {'bands': {'B0': {'stats': {'min': 0.0,
             'max': 92.0,
             'mean': 45.848345889289504,
             'stDev': 30.171668913515024,
             'sampleCount': 51975,
             'noDataCount': 27551}}}}}},
       {'interval': {'from': '2015-01-01T00:00:00Z', 'to': '2015-01-02T00:00:00Z'},
        'outputs': {'biomass': {'bands': {'B0': {'stats': {'min': 0.0,
             'max': 224.7899169921875,
             'mean': 83.46840627038036,
             'stDev': 54.78463705777067,
             'sampleCount': 51975,
             'noDataCount': 27551}}}},
         'carbon': {'bands': {'B0': {'stats': {'min': 0.0,
             'max': 107.0,
             'mean': 39.73096134949235,
             'stDev': 26.07748722782563,
             'sampleCount': 51975,
             'noDataCount': 27551}}}},
         'canopyHeight': {'bands': {'B0': {'stats': {'min': 0.0,
             'max': 27.0,
             'mean': 8.2521290533901,
             'stDev': 8.258077736819507,
             'sampleCount': 51975,
             'noDataCount': 27551}}}},
         'canopyCover': {'bands': {'B0': {'stats': {'min': 0.0,
             'max': 92.0,
             'mean': 44.04966426465771,
             'stDev': 30.509428099625637,
             'sampleCount': 51975,
             'noDataCount': 27551}}}}}},
       {'interval': {'from': '2016-01-01T00:00:00Z', 'to': '2016-01-02T00:00:00Z'},
        'outputs': {'biomass': {'bands': {'B0': {'stats': {'min': 0.0,
             'max': 226.89076232910156,
             'mean': 79.92018469411562,
             'stDev': 55.76700474635435,
             'sampleCount': 51975,
             'noDataCount': 27551}}}},
         'carbon': {'bands': {'B0': {'stats': {'min': 0.0,
             'max': 108.0,
             'mean': 38.04200786111999,
             'stDev': 26.5450942422532,
             'sampleCount': 51975,
             'noDataCount': 27551}}}},
         'canopyHeight': {'bands': {'B0': {'stats': {'min': 0.0,
             'max': 27.0,
             'mean': 7.754995086799883,
             'stDev': 8.338305923157964,
             'sampleCount': 51975,
             'noDataCount': 27551}}}},
         'canopyCover': {'bands': {'B0': {'stats': {'min': 0.0,
             'max': 93.0,
             'mean': 41.00601867016062,
             'stDev': 31.15364418737545,
             'sampleCount': 51975,
             'noDataCount': 27551}}}}}},
       {'interval': {'from': '2017-01-01T00:00:00Z', 'to': '2017-01-02T00:00:00Z'},
        'outputs': {'biomass': {'bands': {'B0': {'stats': {'min': 0.0,
             'max': 207.9832000732422,
             'mean': 61.25681939829456,
             'stDev': 39.49880595576853,
             'sampleCount': 51975,
             'noDataCount': 27551}}}},
         'carbon': {'bands': {'B0': {'stats': {'min': 0.0,
             'max': 99.0,
             'mean': 29.158245987553034,
             'stDev': 18.80143167646003,
             'sampleCount': 51975,
             'noDataCount': 27551}}}},
         'canopyHeight': {'bands': {'B0': {'stats': {'min': 0.0,
             'max': 25.0,
             'mean': 4.993162463151044,
             'stDev': 5.842448668295666,
             'sampleCount': 51975,
             'noDataCount': 27551}}}},
         'canopyCover': {'bands': {'B0': {'stats': {'min': 0.0,
             'max': 90.0,
             'mean': 30.478791352767743,
             'stDev': 27.090565658618637,
             'sampleCount': 51975,
             'noDataCount': 27551}}}}}}],
      'status': 'OK'}]

Sorting the Data​

As you can see in the cell above, the Statistical API returns a json structure with the results of the query. This format isn't easy to plot or filter: to make the data easier to work with, we will convert it into a tabular form using a pandas dataframe.

# Convert to geopandas
forest_df = pd.json_normalize(statistical_response[0]['data'])

# Populate a year column
forest_df['year'] = pd.to_datetime(forest_df['interval.from']).dt.year

# Drop unecessary columns
forest_df = forest_df.drop(columns=['interval.from', 'interval.to'])

# Rename columns
forest_df.columns = forest_df.columns.str.replace('outputs.', '').str.replace('.bands.B0.stats', '')

# Set index
forest_df.set_index('year', inplace=True)

# Show the new Dataframe with the results
forest_df
biomass.minbiomass.maxbiomass.meanbiomass.stDevbiomass.sampleCountbiomass.noDataCountcarbon.mincarbon.maxcarbon.meancarbon.stDev...canopyHeight.meancanopyHeight.stDevcanopyHeight.sampleCountcanopyHeight.noDataCountcanopyCover.mincanopyCover.maxcanopyCover.meancanopyCover.stDevcanopyCover.sampleCountcanopyCover.noDataCount
year
20130.0226.89076294.24312654.25007551975275510.0108.044.85972825.823035...10.2150757.96710451975275510.093.051.64469430.2453845197527551
20140.0224.78991785.53647554.07485651975275510.0107.040.71536225.739631...8.6836728.10869351975275510.092.045.84834630.1716695197527551
20150.0224.78991783.46840654.78463751975275510.0107.039.73096126.077487...8.2521298.25807851975275510.092.044.04966430.5094285197527551
20160.0226.89076279.92018555.76700551975275510.0108.038.04200826.545094...7.7549958.33830651975275510.093.041.00601931.1536445197527551
20170.0207.98320061.25681939.49880651975275510.099.029.15824618.801432...4.9931625.84244951975275510.090.030.47879127.0905665197527551

5 rows × 24 columns

In the following cell, we will plot the average value over the parcel of interest for the 4 different variables returned from the Statistical API query over the time period 2013-2017. The standard deviation is represented by error bars for each year.

Looking at the graph, there appears to be a downward trend in all variables between 2013 and 2017, with a larger drop between 2016 and 2017. For example, the Canopy Cover decreases by around 12% between 2013 and 2016, but then by approximately 14% between 2016 and 2017, suggesting an increase in deforestation.

# Plot average values for the area of interest over time
fig, ax = plt.subplots(figsize=(10, 8), nrows=2, ncols=2, sharex=True)

# Cover
ax[0, 0].errorbar(
    forest_df.index,
    forest_df["canopyCover.mean"],
    forest_df["canopyCover.stDev"],
    line,
    marker=".",
    color="g",
    linewidth=0.5,
)
ax[0, 0].set_title("Canopy Cover")
ax[0, 0].set_ylabel("%")

# Height
ax[0, 1].errorbar(
    forest_df.index,
    forest_df["canopyHeight.mean"],
    forest_df["canopyHeight.stDev"],
    line,
    marker=".",
    color="r",
    linewidth=0.5,
)
ax[0, 1].set_title("Canopy Height")
ax[0, 1].set_ylabel("/ m")

# Carbon
ax[1, 0].errorbar(
    forest_df.index,
    forest_df["carbon.mean"],
    forest_df["carbon.stDev"],
    line,
    marker=".",
    color="b",
    linewidth=0.5,
)
ax[1, 0].set_title("Aboveground carbon density")
ax[1, 0].set_ylabel("/ Mg⋅ha-1")
ax[1, 0].set_xlabel("Year")


# Biomass
ax[1, 1].errorbar(
    forest_df.index,
    forest_df["biomass.mean"],
    forest_df["biomass.stDev"],
    line,
    marker=".",
    color="k",
    linewidth=0.5,
)
ax[1, 1].set_title("Biomass")
ax[1, 1].set_ylabel("/ Mg⋅ha-1")
ax[1, 1].set_xlabel("Year")


for axl in [ax[1, 0], ax[1, 1]]:
    axl.tick_params(axis="x", labelrotation=40)

fig.suptitle("Averaged over AOI")
plt.show()

Compute Sum of Carbon and Biomass​

In a second step, we will use the Processing API to extract the sum of aboveground carbon density and biomass over the parcel of interest. The Statistical API doesn't return the sum of the pixels contained in a geometry, so we will use Process API to write the values to a json file.

We build the request according to the API Reference, using the SentinelHubRequest class from the sentinelhub-py Python package. Each Process API request also needs an evalscript.

The subtlety in this request is that we will pass the cumulative carbon and biomass values to the userdata object that is a property of the outputMetadata, sparing us to have to return an image, and thus only accessing the statistics needed for this analysis.

evalscript_sum = """
//VERSION=3
function setup() {
    return {
        input: ["ACD", "dataMask"],
        output: { bands: 1, sampleTYPE: "UINT8" },
        mosaicking: "ORBIT"
    }
}

// Initialise sum variables
var sumVariables = {};

// Set constants
const BIOMASS_CONV = 0.476

function updateOutputMetadata(scenes, inputMetadata, outputMetadata) {
    // Write the output to a metadata file
    outputMetadata.userData = sumVariables
}

function evaluatePixel(samples, scenes, inputMetadata, customData, outputMetadata) {

    // Loop around samples
    samples.forEach((sample, index) => {
        if (sample.dataMask == 1) {
            // Populate the Object containing sum values by adding the pixel value or creating the entry if it doesn't exist
            if (sumVariables.hasOwnProperty(scenes[index].date.getFullYear())) {
                sumVariables[scenes[index].date.getFullYear()]["carbon"]["sum"] += sample.ACD
                sumVariables[scenes[index].date.getFullYear()]["biomass"]["sum"] += sample.ACD / BIOMASS_CONV
            } else {
                sumVariables[scenes[index].date.getFullYear()] = { "carbon": { "sum": sample.ACD }, "biomass": { "sum": sample.ACD / BIOMASS_CONV } }
            }
        }
    });
}
"""
# Payload
carbonsum = SentinelHubRequest(
    evalscript=evalscript_sum,
    input_data=[
        # Set the collection containing Aboveground Carbon Density, Forest Carbon Diligence
        # Here we fetch data from 2013 to 2017
        SentinelHubRequest.input_data(
            data_collection=carbon_collection,
            time_interval=(f"{start_year}-01-01", f"{end_year}-01-01"),
        )
    ],
    # We return 2 outputs: "default" which will be empty and "userdata" which will contain the metadata requested in the Evalscript
    responses=[
        SentinelHubRequest.output_response("default", MimeType.PNG),
        SentinelHubRequest.output_response("userdata", MimeType.JSON),
    ],
    geometry=parcel_geo,
    # Match the resolution of the PV (30 meters). Beware: resolution should be in the units of the coordinate system used
    resolution=(30, 30),
    data_folder=".",
    config=config,
)
# Run the request and display a sample of the data returned
sum_response = carbonsum.get_data()
sum_response[0]['userdata.json']
Output
    {'2013': {'carbon': {'sum': 1095654}, 'biomass': {'sum': 2301794.1176472562}},
     '2014': {'carbon': {'sum': 994432}, 'biomass': {'sum': 2089142.8571429157}},
     '2015': {'carbon': {'sum': 970389}, 'biomass': {'sum': 2038632.3529412486}},
     '2016': {'carbon': {'sum': 929138}, 'biomass': {'sum': 1951970.5882353636}},
     '2017': {'carbon': {'sum': 712161}, 'biomass': {'sum': 1496136.5546219142}}}

Just as we did for the results coming from the Statistical API, we will convert the results into a tabular form using a pandas dataframe. To insert the data into the main dataframe, we will follow the same structure.

sum_response[0]
Output
    {'default.png': array([[0, 0, 0, ..., 0, 0, 0],
            [0, 0, 0, ..., 0, 0, 0],
            [0, 0, 0, ..., 0, 0, 0],
            ...,
            [0, 0, 0, ..., 0, 0, 0],
            [0, 0, 0, ..., 0, 0, 0],
            [0, 0, 0, ..., 0, 0, 0]], dtype=uint8),
     'userdata.json': {'2013': {'carbon': {'sum': 1095654},
       'biomass': {'sum': 2301794.1176472562}},
      '2014': {'carbon': {'sum': 994432}, 'biomass': {'sum': 2089142.8571429157}},
      '2015': {'carbon': {'sum': 970389}, 'biomass': {'sum': 2038632.3529412486}},
      '2016': {'carbon': {'sum': 929138}, 'biomass': {'sum': 1951970.5882353636}},
      '2017': {'carbon': {'sum': 712161}, 'biomass': {'sum': 1496136.5546219142}}}}
# Convert to geopandas
carbon_df = pd.DataFrame()

for x in sum_response[0]['userdata.json']:
    df = pd.json_normalize(sum_response[0]['userdata.json'][x])
    df["year"] = int(x)
    carbon_df = pd.concat([carbon_df, df])

df = None
carbon_df.set_index('year', inplace=True)

# Display the dataFrame
carbon_df
carbon.sumbiomass.sum
year
201310956542.301794e+06
20149944322.089143e+06
20159703892.038632e+06
20169291381.951971e+06
20177121611.496137e+06

Merge Dataframes from Statistical and Processing API​

In the next cell, we will integrate the newly created Dataframe into the one created earlier from the Statistical API results.

statistics_df = forest_df.join(carbon_df)
statistics_df
biomass.minbiomass.maxbiomass.meanbiomass.stDevbiomass.sampleCountbiomass.noDataCountcarbon.mincarbon.maxcarbon.meancarbon.stDev...canopyHeight.sampleCountcanopyHeight.noDataCountcanopyCover.mincanopyCover.maxcanopyCover.meancanopyCover.stDevcanopyCover.sampleCountcanopyCover.noDataCountcarbon.sumbiomass.sum
year
20130.0226.89076294.24312654.25007551975275510.0108.044.85972825.823035...51975275510.093.051.64469430.245384519752755110956542.301794e+06
20140.0224.78991785.53647554.07485651975275510.0107.040.71536225.739631...51975275510.092.045.84834630.17166951975275519944322.089143e+06
20150.0224.78991783.46840654.78463751975275510.0107.039.73096126.077487...51975275510.092.044.04966430.50942851975275519703892.038632e+06
20160.0226.89076279.92018555.76700551975275510.0108.038.04200826.545094...51975275510.093.041.00601931.15364451975275519291381.951971e+06
20170.0207.98320061.25681939.49880651975275510.099.029.15824618.801432...51975275510.090.030.47879127.09056651975275517121611.496137e+06

5 rows × 26 columns

In the next cell we will plot the sum of aboveground carbon and biomass (linearly derived from carbon) for the parcel of interest.

Similarly to the previous figure, there seems to be a downward trend in values over time, with a larger drop between 206 and 2017.

# Plot average values for the area of interest over time
fig, ax = plt.subplots(figsize=(10, 8), nrows=1, ncols=1, sharex=True)

# Carbon
ax.plot(
    statistics_df["carbon.sum"],
    line,
    marker=".",
    color="g",
    linewidth=0.5,
    label="Carbon"
)
ax.set_ylim(0, 3e6)

# Biomass
ax.plot(
    statistics_df["biomass.sum"],
    line,
    marker=".",
    color="r",
    linewidth=0.5,
    label= "Biomass"
)

ax.set_ylabel("/ Mg⋅ha-1")
ax.set_xlabel("Year")

plt.gca().legend()
plt.suptitle("Total aboveground carbon and biomass")
plt.show()

Download to CSV​

In the following cell, we will download the results stored in the pandas dataframe as a .csv file. If you wish to download the raw data from the requests (Statistical or Processing APIs), simply add the save_data=True flag to when executing them, such as:

  • sum_response = carbonsum.get_data(save_data=True)
  • statistical_response = request.get_data(save_data=True)
# Download pandas as csv
statistics_df.to_csv("./fcd_statistics.csv")

Extract Statistics Filtering for Forest Cover​

In the following cells we will extract statistics using the Statistical API and parse the results into a Dataframe, as we did in Step 1. The difference here is that within the Evalscript we will filter out results based on a binary forest mask.

The binary forest mask is based on user-specified thresholds for height and cover, and can be adjusted in the Evalscript.

After converting the results to a pandas Dataframe, we will compare the statistics with the previous data that don't account for the presence of forested / non-forested areas. This will provide better insights into the evolution of the forested cover over the parcel of interest.

evalscript_masked_stats = """
//VERSION=3

// Set threshold for height (in meters) and cover (in %) to map forest
const heightThreshold = 5
const coverThreshold = 50

// Set constants
const BIOMASS_CONV = 0.476

function setup() {
    return {
        input: [{
            datasource: "canopyHeight",
            bands: ["CH", "dataMask"]
        },
        {
            datasource: "canopyCover",
            bands: ["CC", "dataMask"]
        },
        {
            datasource: "carbon",
            bands: ["ACD", "dataMask"]
        }
        ],
        output: [
            { id: "canopyHeight", bands: 1, sampleType: "UINT8" },
            { id: "canopyCover", bands: 1, sampleType: "UINT8" },
            { id: "carbon", bands: 1, sampleType: "UINT16" },
            { id: "biomass", bands: 1, sampleType: "FLOAT32" },
            { id: "dataMask", bands: 1, sampleType: "UINT8" }
        ],
        mosaicking: "ORBIT"
    };
}

function evaluatePixel(samples) {
    var sampleHeight = samples.canopyHeight[0].CH
    var nodataHeight = samples.canopyHeight[0].dataMask
    var sampleCover = samples.canopyCover[0].CC
    var nodataCover = samples.canopyCover[0].dataMask
    var sampleCarbon = samples.carbon[0].ACD
    var nodataCarbon = samples.carbon[0].dataMask

    // Biomass simple scalar conversion
    var sampleBiomass = sampleCarbon / BIOMASS_CONV

    // Compute forest mask
    var forestMask = 0
    if (sampleHeight >= heightThreshold && sampleCover >= coverThreshold) {
        // Where the map matches criteria, set to 1
        forestMask = 1
    }

    // Combined noData / Forest mask
    var combinedMask = forestMask * nodataHeight * nodataCover * nodataCarbon

    return {
        canopyHeight: [sampleHeight],
        canopyCover: [sampleCover],
        carbon: [sampleCarbon],
        biomass: [sampleBiomass],
        dataMask: [combinedMask]
    };
}
"""
# Payload
masked_request = SentinelHubStatistical(
    aggregation=SentinelHubStatistical.aggregation(
        evalscript=evalscript_masked_stats,
        time_interval=(f"{start_year}-01-01T00:00:00Z", f"{end_year}-12-31T23:59:59Z"),
        # To fetch the yearly data that is set to 1 date / year we set the aggregation to P1D
        aggregation_interval="P1D",
        # Match the resolution of the PV (30 meters). Beware: resolution should be in the units of the coordinate system used
        resolution=[30, 30],
    ),
    input_data=[
        SentinelHubStatistical.input_data(
            height_collection,
            identifier="canopyHeight",
        ),
        SentinelHubStatistical.input_data(
            cover_collection,
            identifier="canopyCover",
        ),
        SentinelHubStatistical.input_data(
            carbon_collection,
            identifier="carbon",
        ),
    ],
    geometry=parcel_geo,
    config=config,
)
# Run the request and display a sample of the data returned
statistical_response_masked = masked_request.get_data()
statistical_response_masked
Output
    [{'data': [{'interval': {'from': '2013-01-01T00:00:00Z',
         'to': '2013-01-02T00:00:00Z'},
        'outputs': {'biomass': {'bands': {'B0': {'stats': {'min': 54.621849060058594,
             'max': 226.89076232910156,
             'mean': 128.08722994031922,
             'stDev': 41.05462639260904,
             'sampleCount': 51975,
             'noDataCount': 37473}}}},
         'carbon': {'bands': {'B0': {'stats': {'min': 26.0,
             'max': 108.0,
             'mean': 60.9695214453179,
             'stDev': 19.542002065911515,
             'sampleCount': 51975,
             'noDataCount': 37473}}}},
         'canopyHeight': {'bands': {'B0': {'stats': {'min': 5.0,
             'max': 27.0,
             'mean': 15.597848572610692,
             'stDev': 5.350587433451934,
             'sampleCount': 51975,
             'noDataCount': 37473}}}},
         'canopyCover': {'bands': {'B0': {'stats': {'min': 50.0,
             'max': 93.0,
             'mean': 74.26410150324106,
             'stDev': 10.603667601384297,
             'sampleCount': 51975,
             'noDataCount': 37473}}}}}},
       {'interval': {'from': '2014-01-01T00:00:00Z', 'to': '2014-01-02T00:00:00Z'},
        'outputs': {'biomass': {'bands': {'B0': {'stats': {'min': 52.5210075378418,
             'max': 224.7899169921875,
             'mean': 129.29616386274668,
             'stDev': 44.27644451107587,
             'sampleCount': 51975,
             'noDataCount': 40435}}}},
         'carbon': {'bands': {'B0': {'stats': {'min': 25.0,
             'max': 107.0,
             'mean': 61.54497400346619,
             'stDev': 21.075587449728328,
             'sampleCount': 51975,
             'noDataCount': 40435}}}},
         'canopyHeight': {'bands': {'B0': {'stats': {'min': 5.0,
             'max': 27.0,
             'mean': 15.926603119584026,
             'stDev': 5.710512310895107,
             'sampleCount': 51975,
             'noDataCount': 40435}}}},
         'canopyCover': {'bands': {'B0': {'stats': {'min': 50.0,
             'max': 92.0,
             'mean': 73.90433275563224,
             'stDev': 10.642962473141507,
             'sampleCount': 51975,
             'noDataCount': 40435}}}}}},
       {'interval': {'from': '2015-01-01T00:00:00Z', 'to': '2015-01-02T00:00:00Z'},
        'outputs': {'biomass': {'bands': {'B0': {'stats': {'min': 52.5210075378418,
             'max': 224.7899169921875,
             'mean': 131.53803897503903,
             'stDev': 44.783645391552696,
             'sampleCount': 51975,
             'noDataCount': 41204}}}},
         'carbon': {'bands': {'B0': {'stats': {'min': 25.0,
             'max': 107.0,
             'mean': 62.61210658248998,
             'stDev': 21.317015069941753,
             'sampleCount': 51975,
             'noDataCount': 41204}}}},
         'canopyHeight': {'bands': {'B0': {'stats': {'min': 5.0,
             'max': 27.0,
             'mean': 16.242967226812812,
             'stDev': 5.823265573954096,
             'sampleCount': 51975,
             'noDataCount': 41204}}}},
         'canopyCover': {'bands': {'B0': {'stats': {'min': 50.0,
             'max': 92.0,
             'mean': 74.24649521864289,
             'stDev': 10.712455971948414,
             'sampleCount': 51975,
             'noDataCount': 41204}}}}}},
       {'interval': {'from': '2016-01-01T00:00:00Z', 'to': '2016-01-02T00:00:00Z'},
        'outputs': {'biomass': {'bands': {'B0': {'stats': {'min': 52.5210075378418,
             'max': 226.89076232910156,
             'mean': 132.95173141729768,
             'stDev': 45.410074429458604,
             'sampleCount': 51975,
             'noDataCount': 42039}}}},
         'carbon': {'bands': {'B0': {'stats': {'min': 25.0,
             'max': 108.0,
             'mean': 63.28502415458952,
             'stDev': 21.615195299454385,
             'sampleCount': 51975,
             'noDataCount': 42039}}}},
         'canopyHeight': {'bands': {'B0': {'stats': {'min': 5.0,
             'max': 27.0,
             'mean': 16.508655394524943,
             'stDev': 5.918522669669513,
             'sampleCount': 51975,
             'noDataCount': 42039}}}},
         'canopyCover': {'bands': {'B0': {'stats': {'min': 50.0,
             'max': 93.0,
             'mean': 74.26519726248002,
             'stDev': 10.928273162257332,
             'sampleCount': 51975,
             'noDataCount': 42039}}}}}},
       {'interval': {'from': '2017-01-01T00:00:00Z', 'to': '2017-01-02T00:00:00Z'},
        'outputs': {'biomass': {'bands': {'B0': {'stats': {'min': 54.621849060058594,
             'max': 207.9832000732422,
             'mean': 108.69701512160464,
             'stDev': 35.06227180215333,
             'sampleCount': 51975,
             'noDataCount': 45542}}}},
         'carbon': {'bands': {'B0': {'stats': {'min': 26.0,
             'max': 99.0,
             'mean': 51.739779263174356,
             'stDev': 16.6896414084141,
             'sampleCount': 51975,
             'noDataCount': 45542}}}},
         'canopyHeight': {'bands': {'B0': {'stats': {'min': 5.0,
             'max': 25.0,
             'mean': 13.217938753303256,
             'stDev': 4.767180957412843,
             'sampleCount': 51975,
             'noDataCount': 45542}}}},
         'canopyCover': {'bands': {'B0': {'stats': {'min': 50.0,
             'max': 90.0,
             'mean': 69.75376962536926,
             'stDev': 10.075603920119867,
             'sampleCount': 51975,
             'noDataCount': 45542}}}}}}],
      'status': 'OK'}]

Sorting the Data​

Similarly to the previous Statistical API request, in the next cell we will convert the json response into a tabular form using a pandas dataframe.

# Convert to geopandas
masked_forest_df = pd.json_normalize(statistical_response_masked[0]['data'])

# Populate a year column
masked_forest_df['year'] = pd.to_datetime(masked_forest_df['interval.from']).dt.year

# Drop unecessary columns
masked_forest_df = masked_forest_df.drop(columns=['interval.from', 'interval.to'])

# Rename columns
masked_forest_df.columns = masked_forest_df.columns.str.replace('outputs.', '').str.replace('.bands.B0.stats', '')

# Set index
masked_forest_df.set_index('year', inplace=True)

# Show the new Dataframe with the results
masked_forest_df
biomass.minbiomass.maxbiomass.meanbiomass.stDevbiomass.sampleCountbiomass.noDataCountcarbon.mincarbon.maxcarbon.meancarbon.stDev...canopyHeight.meancanopyHeight.stDevcanopyHeight.sampleCountcanopyHeight.noDataCountcanopyCover.mincanopyCover.maxcanopyCover.meancanopyCover.stDevcanopyCover.sampleCountcanopyCover.noDataCount
year
201354.621849226.890762128.08723041.054626519753747326.0108.060.96952119.542002...15.5978495.350587519753747350.093.074.26410210.6036685197537473
201452.521008224.789917129.29616444.276445519754043525.0107.061.54497421.075587...15.9266035.710512519754043550.092.073.90433310.6429625197540435
201552.521008224.789917131.53803944.783645519754120425.0107.062.61210721.317015...16.2429675.823266519754120450.092.074.24649510.7124565197541204
201652.521008226.890762132.95173145.410074519754203925.0108.063.28502421.615195...16.5086555.918523519754203950.093.074.26519710.9282735197542039
201754.621849207.983200108.69701535.062272519754554226.099.051.73977916.689641...13.2179394.767181519754554250.090.069.75377010.0756045197545542

5 rows × 24 columns

In the next cell we will plot the comparison between our two Statistical API queries: one for all pixels located in the parcel of interest, and the other only accounting for forest-covered pixels based on user-defined thresholds. We will plot Canopy Height and aboveground carbon.

The plot shows similar trends for both variables, although it seems that the forested areas were not so impacted by carbon loss over the period 2013 - 2016 as when we account for all types of land. However, in the last year of the time-series, we observe a loss in carbon and a significant drop in average canopy height. Of course, this analysis is for illustrative purposes only and the interpretation of the results would require a more thorough investigation.

# Plot average values for the area of interest over time
fig, ax = plt.subplots(figsize=(10, 8), nrows=2, ncols=1, sharex=True)

# Canopy Height
ax[0].plot(
    statistics_df["canopyHeight.mean"],
    line,
    marker=".",
    color="g",
    linewidth=0.5,
    label="No forest mask applied"
)
ax[0].plot(
    masked_forest_df["canopyHeight.mean"],
    line,
    marker=".",
    color="r",
    linewidth=0.5,
    label= "Forest mask applied"
)
ax[0].set_ylim(2, 20)
ax[0].set_ylabel("/ m")
ax[0].legend()
ax[0].set_title("Average Canopy Height")

# Carbon
ax[1].plot(
    statistics_df["carbon.mean"],
    line,
    marker=".",
    color="g",
    linewidth=0.5,
    label="No forest mask applied"
)
ax[1].plot(
    masked_forest_df["carbon.mean"],
    line,
    marker=".",
    color="r",
    linewidth=0.5,
    label= "Forest mask applied"
)
ax[1].set_ylim(0, 80)
ax[1].set_ylabel("/ Mg⋅ha-1")
ax[1].set_xlabel("Year")
ax[1].legend()
ax[1].set_title("Aboveground carbon density")

plt.show()

Canopy Height and Cover Maps​

In the following cell, we will try to visually interpret the maps of Canopy Cover and Canopy Height, comparing the results in 2013 (start of the time-series), 2016 and 2017.

In a first step, we will run a Processing API request to collect data over the period 2013-2017 for the two variables. Although we could simply request a visualisation of the variables from Sentinel Hub, in this case we will request the raw values for further analysis needs.

In a second step we will plot the images.

evalscript_maps = """
//VERSION=3
function setup() {
    return {
        input: [{
            datasource: "canopyHeight",
            bands: ["CH", "dataMask"]
        },
        {
            datasource: "canopyCover",
            bands: ["CC", "dataMask"]
        }],
        output: [
            { id: "canopyHeight", bands: 5, sampleType: "UINT16" },
            { id: "canopyCover", bands: 5, sampleType: "UINT16" },
        ],
        mosaicking: "ORBIT"
    };
}

function evaluatePixel(samples) {

    var sampleHeight = samples.canopyHeight
    var sampleCover = samples.canopyCover

    // Loop around samples
    let heightValue = []
    let coverValue = []

    sampleHeight.forEach((sample, index) => {
        if (sample.dataMask == 1) {
            heightValue.push(sample.CH)
        } else {
            heightValue.push(255)
        }
    });
    sampleCover.forEach((sample, index) => {
        if (sample.dataMask == 1) {
            coverValue.push(sample.CC)
        } else {
            coverValue.push(255)
        }
    });
    return {
        canopyHeight: heightValue,
        canopyCover: coverValue
    }
}
"""
# Payload
request_raw = SentinelHubRequest(
    evalscript=evalscript_maps,
    input_data=[
        SentinelHubStatistical.input_data(
            height_collection,
            identifier="canopyHeight",
            time_interval=("2013-01-01", "2017-12-31"),
            mosaicking_order="leastRecent",
        ),
        SentinelHubStatistical.input_data(
            cover_collection,
            identifier="canopyCover",
            time_interval=("2013-01-01", "2017-12-31"),
            mosaicking_order="leastRecent",
        ),
    ],
    responses=[
        SentinelHubRequest.output_response("canopyHeight", MimeType.TIFF),
        SentinelHubRequest.output_response("canopyCover", MimeType.TIFF),
    ],
    geometry=parcel_geo,
    resolution=(30, 30),
    data_folder=".",
    config=config,
)
# Run the request and display a sample of the data returned
fcd_maps = request_raw.get_data()

Plot Maps​

In the next three cells we will plot Canopy Cover and height over the time period 2013-2017.

# First we convert 255 (noData) to transparent for plotting purposes.
height_maps = fcd_maps[0]["canopyHeight.tif"]
cover_maps = fcd_maps[0]["canopyCover.tif"]
height_maps = np.ma.masked_where(height_maps == 255, height_maps)
cover_maps = np.ma.masked_where(cover_maps == 255, cover_maps)

Canopy Cover​

By plotting the Canopy Cover over the the entire time-series, we notice the overall decrease of canopy cover over time, with a sudden drop between 2016 and 2017. This confirms the statistical data that we plotted previously.

# Plot average values for the area of interest over time
fig, ax = plt.subplots(figsize=(20, 5), nrows=1, ncols=5, sharey=True)

sy = 2013
for i in range(height_maps.shape[2]):
    im = ax[i].imshow(cover_maps[:, :, i], cmap="Greens", vmin=0, vmax=100)
    ax[i].set_title(str(sy), fontsize=12)
    sy += 1
    ax[i].set_yticks([])
    ax[i].set_xticks([])

fig.subplots_adjust(bottom=0.2, top=0.8, wspace=0.3)
cbar_ax = fig.add_axes([0.2, 0.1, 0.6, 0.03])
cbar = fig.colorbar(im, cax=cbar_ax, orientation="horizontal")
cbar.set_label("Canopy Cover (%)", fontsize=12)

plt.suptitle("Canopy Cover 2013-2017", fontsize=16)
plt.show()

Canopy Height​

The successive losses in Canopy between 2013 and 2017 are also visible when looking at the Canopy Height maps for the parcel of interest. However, it is interesting to note that in the eastern part of the AOI, there is a significant drop in Canopy height without there being deforestation (halving of the average canopy height). More analysis would be needed to understand the causes of this phenomenon.

# Plot average values for the area of interest over time
fig, ax = plt.subplots(figsize=(20, 5), nrows=1, ncols=5, sharey=True)

sy = 2013
for i in range(height_maps.shape[2]):
    im = ax[i].imshow(height_maps[:, :, i], cmap="viridis", vmin=0, vmax=30)
    ax[i].set_title(str(sy), fontsize=12)
    sy += 1
    ax[i].set_yticks([])
    ax[i].set_xticks([])

fig.subplots_adjust(bottom=0.2, top=0.8, wspace=0.3)
cbar_ax = fig.add_axes([0.2, 0.1, 0.6, 0.03])
cbar = fig.colorbar(im, cax=cbar_ax, orientation="horizontal")
cbar.set_label("Canopy Height (m)", fontsize=12)

plt.suptitle("Canopy Height 2013-2017", fontsize=16)
plt.show()