import numpy as np
import pandas as pd
import xarray as xr
import openeo
from utils import (
set_plot_defaults,
plot_harvest_timing_map,
plot_coherence_samples,
plot_interactive_harvest_map,
)Analysing Coherence output generated using openEO for harvest time detection
The main goal of this notebook is to demonstrate how to use openEO to generate Sentinel-1 SAR coherence, and then use Python/xarray to further analyse the resulting coherence time series outside the openEO environment. The notebook is structured as follows:
- openEO processing: We will use openEO to generate a coherence time series from Sentinel-1 SAR data.
- Python/xarray analysis: We will then use Python and xarray to analyse the generated coherence time series, focusing on detecting harvest times.
The workflow is designed to leverage the strengths of both openEO and local Python analysis by letting openEO handle access to large EO data and the expensive SAR-coherence processing close to the data. And then, using the local Python environment for more flexible and interactive analysis of the results.
Let us start by importing the necessary libraries for the analysis.
Connect to an openEO back end
The connection below provides the notebook with access to openEO’s CDSE backend, where we can run the coherence processing workflow.
connection = openeo.connect("openeo.dataspace.copernicus.eu").authenticate_oidc()Authenticated using refresh token.Fetching the coherence time series
Now that we have a connection to the openEO backend, we can call the sentinel1_sar_coherence User-Defined-Process (UDP) to compute the coherence between pairs of Sentinel-1 SAR images.
CWL in openEO
The sentinel1_sar_coherence UDP used below does not reimplement the SAR coherence algorithm in openEO itself. Instead, it wraps an existing Common Workflow Language (CWL) pipeline built on ESA SNAP operators, developed by Eurac Research as part of the CloudInSAR project (ESA project page). CWL lets openEO backends run external, containerized scientific workflows as if they were native processes, so complex tools don’t need to be rewritten as UDFs.
Support for CWL-based processes in openEO is still experimental: it is currently limited to specific backends (e.g. the CDSE federation, via Eurac’s implementation).
This process generates a time series of Sentinel‑1 interferometric coherence for a set of interferometric pairs. By simply calling the UDP, we avoid the need to manually download Sentinel-1 data and perform the complex coherence processing workflow. Moreover, the practical benefit is that all the heavy SAR coherence processing is handled inside openEO, letting researchers focus on the scientific analysis of the results.
In this example, we will request the coherence for the VH polarisation, which is often more sensitive to vegetation changes. The coherence time series will be generated for a specific area of interest and time period. Users can adjust parameters to suit their specific needs, such as the area of interest, time period, or polarisation.
Current limitation: At present, this process can only be used as the final step of an openEO workflow. The workflow returns a STAC result rather than a native openEO datacube. This means it can currently only be used as the final process-graph node
cube = connection.datacube_from_process(
"sentinel1_sar_coherence",
namespace="https://raw.githubusercontent.com/ESA-APEx/apex_algorithms/main/algorithm_catalog/eurac/sentinel1_sar_coherence/openeo_udp/sentinel1_sar_coherence.json",
temporal_extent=["2023-04-01", "2023-09-30"],
temporal_baseline=12,
spatial_extent={
"west": 5.55,
"south": 52.55,
"east": 5.65,
"north": 52.60
},
polarization= "VH",
sub_swath="IW2",
)
cubeThe batch job below submits a coherence processing request to the CDSE backend and saves the resulting NetCDF file, which provides the stacked coherence time series for the selected area and time period.
cube.execute_batch(title="Sentinel-1 Coherence", outputfile="s1_coherence.nc")If you inspect the job information, once completed, you will notice the several indicators of the processing load and resource usage. These are important factors that show all the heavy processing done in the cloud infrastructure. Furthermore, for this particular job, only 50 credits were used out of the 10000 free credits that are replenished every month.
The generated coherence timeseries can be used to detect harvest times by analysing the temporal changes in coherence values. Coherence measures how similar two radar images are over time. A big jump in coherence usually means the field has been harvested (because the field changes from crops to bare soil).
Traditional monitoring uses sensors on harvesters, which is costly and complex. Using a satellite data-based approach helps monitor large farming areas without needing expensive equipment on farm machinery.
At this point, the computationally intensive SAR coherence generation is finished. We now deliberately switch to local Python. The local analysis below is intentionally simple, to answer a higher-level question such as ” when and where might a harvest have occurred?“
Load the processed coherence data
In the previous steps, we saved the coherence time series as a NetCDF file. Now, we will load this file into our local Python environment using xarray for further analysis.
# load the downloaded coherence dataset
ds = xr.open_dataset("s1_coherence.nc")
print(f"Loaded dataset: {dict(ds.sizes)}")Loaded dataset: {'t': 14, 'y': 5076, 'x': 10640}The cube contains 14 acquisitions × 5076 × 10640 spatial pixels.
This is already a useful illustration of the advantage of the previous step: we are working locally with a ready-to-use coherence product instead of having to download and prepare the original Sentinel-1 observations and implement the coherence algorithm ourselves.
Additionally, since the result was saved as a NetCDF file, we can skip the additional step of mosaicking, reprojecting and stacking the individual output.
# select the data variables that have both x and y dimensions
data_vars = [v for v in ds.data_vars if "x" in ds[v].dims and "y" in ds[v].dims]
print(f"Found {len(data_vars)} data variable(s) with spatial dimensions: {data_vars}")
coh_vh = ds[data_vars[0]]Found 1 data variable(s) with spatial dimensions: ['coh_VH']# just to be sure, sort the data by time
coh_vh = coh_vh.sortby("t")
print(f"Loaded VH coherence cube: {dict(coh_vh.sizes)}")
print(f"Time range: {pd.Timestamp(coh_vh.t.min().values).date()} to "
f"{pd.Timestamp(coh_vh.t.max().values).date()} ({coh_vh.sizes['t']} acquisitions)")Loaded VH coherence cube: {'t': 14, 'y': 5076, 'x': 10640}
Time range: 2023-04-12 to 2023-09-15 (14 acquisitions)Calculate harvest change detection
The next step calculates the change in coherence between consecutive acquisitions and interprets potential harvest dates by thresholding the coherence change. A similar idea has also been used in the paper: Monitoring Harvesting by Time Series of Sentinel-1 SAR Data
Thus, in this notebook, we will use a simple thresholding approach to detect potential harvest events based on changes in coherence between consecutive acquisitions. The threshold value (here set to Z_THRESHOLD = 2.0) can be adjusted based on the specific characteristics of the monitored area and the expected changes in coherence due to harvesting activities.
# Cleaner defaults for the figures further down (see utils.py)
set_plot_defaults()#how unusual a coherence jump must be to count as a harvest event
Z_THRESHOLD = 2.0
jump = coh_vh.diff("t")
# z-score formula: (x - mean) / stddev
jump_z = (jump - jump.mean("t")) / (jump.std("t") + 1e-6)A sharp increase in VH coherence between two acquisitions often indicates a change in the surface, such as a crop or canopy being cleared. Therefore, a large z-scored jump is used here as a simple harvest signal. Please note that, for demonstration purposes, only VH polarisation is used here. You can also run the same analysis for VV polarisation, or combine both polarisations to improve detection accuracy.
Next, we convert the change dates to day-of-year values in the following cell that will allow us to visualize the detected harvest events in a more interpretable format. For each pixel, the first change exceeding the threshold is recorded as a possible harvest or clearing date.
# day of the year for each acquisition
doy = jump_z["t"].dt.dayofyear
print(f"DOY range: {doy.min().values} to {doy.max().values} ({doy.size} acquisitions)")
# find the acquisition with the first z-score jump above the threshold for each pixel
harvest_doy = doy.where(jump_z > Z_THRESHOLD).min(dim="t", skipna=True)
print(f"Harvest candidates: {int(np.isfinite(harvest_doy.values).sum())} of {harvest_doy.size} pixels crossed the z-score threshold")DOY range: 114 to 258 (13 acquisitions)
Harvest candidates: 3577936 of 54008640 pixels crossed the z-score thresholdTo further evaluate the detected candidate harvest events, we will calculate two additional metrics: * the maximum jump: the largest detected change in coherence for each pixel, and * the mean coherence: the average coherence value for each pixel over the entire time series.
These metrics help us to distinguish between true harvest events and other changes in the landscape that may not be related to harvesting. For example, a pixel with a large maximum jump but low mean coherence may indicate a temporary change in the surface, while a pixel with both a large maximum jump and high mean coherence is more likely to represent a true harvest event.
# some usual statistics
max_jump_zscore = jump_z.max(dim="t", skipna=True)
mean_vh_coherence = coh_vh.mean(dim="t", skipna=True)
print(f"Max jump z-score range: {max_jump_zscore.min().values:.2f} to {max_jump_zscore.max().values:.2f}")
# also check for the valid pixels in the harvest_doy array
n_flagged = int(np.isfinite(harvest_doy.values).sum())Max jump z-score range: 0.00 to 3.45# the time range of the data
times = pd.to_datetime(coh_vh.t.values)
start_day, end_day = times.dayofyear.min(), times.dayofyear.max()
month_dates = pd.date_range(times.min().replace(day=1), times.max(), freq="MS")
month_days = [max(start_day, date.dayofyear) for date in month_dates]
month_labels = [date.strftime("%b") for date in month_dates]
print(f"Month labels: {month_labels} (DOY {month_days})")
# spatial extent of the data
extent = [float(coh_vh.x.min()), float(coh_vh.x.max()), float(coh_vh.y.min()), float(coh_vh.y.max())]
print(f"Spatial extent: {extent[0]:.3f}, {extent[1]:.3f}, {extent[2]:.3f}, {extent[3]:.3f}")Month labels: ['Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep'] (DOY [np.int32(102), 121, 152, 182, 213, 244])
Spatial extent: 244175.000, 350565.000, 5796655.000, 5847405.000Plot the harvest-timing map
The map below shows where candidate harvest or clearing events were detected and when they occurred. Please note that the detected events are based on a simple thresholding approach and may not precisely represent actual harvest events.
# plot the spatial map of the detected harvest timing (see utils.py)
plot_harvest_timing_map(harvest_doy, extent, start_day, end_day, month_days, month_labels)
(<Figure size 900x600 with 2 Axes>,
<Axes: title={'center': 'Detected Harvest Timing'}, xlabel='X (m)', ylabel='Y (m)'>)Though a map is a good way to visualize the spatial distribution of the event, inspecting pixel-level time series can provide a more detailed view of the changes over time.
Therefore, in the cell below, we try to plot the coherence time series for a few representative pixels. This allows us to visually inspect the coherence evolution and verify whether the detected harvest events correspond to significant changes in the time series.
# Sample up to 3 pixels that have a detected harvest date
rng = np.random.default_rng(0)
ys, xs = np.where(np.isfinite(harvest_doy.values))
if len(ys) == 0:
print("No pixel crossed the z-score threshold - showing the AOI center pixel instead.")
sample_idx = [(coh_vh.sizes["y"] // 2, coh_vh.sizes["x"] // 2)]
else:
pick = rng.choice(len(ys), size=min(3, len(ys)), replace=False)
sample_idx = list(zip(ys[pick], xs[pick]))
# plot the coherence time series for the sampled pixels (see utils.py)
plot_coherence_samples(coh_vh, harvest_doy, sample_idx)
(<Figure size 900x500 with 1 Axes>,
<Axes: title={'center': 'VH Coherence Time Series'}, xlabel='Date', ylabel='VH coherence'>)Also, it can be interesting to explore the interactive map of detected harvest events. The interactive map adds geographic context and hover information. Comparing the detections with a basemap helps assess whether they correspond to plausible agricultural areas.
The plot shows VH coherence for three pixels from April to September 2023. Sharp increases occur on different dates, suggesting possible harvest or vegetation-clearing events at those locations. The dashed vertical lines mark the detected candidate harvest dates based on the z-score threshold.
# interactive harvest-timing map with a basemap (see utils.py)
plot_interactive_harvest_map(harvest_doy, start_day, end_day)
ⓘ