Data Client

This notebook shows how to use the dataclient Python library to retrieve timeseries data from the time series service in a very fast manner.

For developers whishing to implement their own web service calls, you can find the REST documentation here: https://proba-v-mep.esa.int/api/timeseries/apidocs/

Defining area of interest

We first need to have a polygon that defines our area of interest. We'll read if from a GeoJson file and show it on the map.

Note that this map may not be visible in the exported version of this notebook!

In [1]:
import geopandas as gpd
import pandas as pd
polygons = gpd.read_file("reserves.geojson")
test_area = polygons[polygons.Reservaat == "De Maat"]
test_area = test_area.iloc(0)[0]
In [2]:
from ipyleaflet import Map,DrawControl,GeoJSON
from shapely.geometry import mapping, shape

m = Map(center=[test_area.geometry.centroid.y,test_area.geometry.centroid.x], zoom=11)
dc = DrawControl()
m.add_control(dc)

m.add_layer(GeoJSON(data=mapping(test_area.geometry)))

m

This notebook uses the MEP dataclient library. It is available at https://bitbucket.org/vitotap/mep-dataclient.

Now we want to retrieve a timeseries for our polygon. Before we can do this, we need to know which 'layer' to query. This list of layers can also be retrieved from the web service:

In [3]:
import requests
response = requests.get("https://proba-v-mep.esa.int/api/timeseries/v1.0/ts/")
if response.status_code == 200:
    layerlist = response.json()['layers']
    layers = [layer['name'] for layer in layerlist if layer['name'].startswith('S2_') or layer['name'].startswith('PROBA') or layer['name'].startswith('BIOPAR')]
else:
    raise IOError(response.text)
layers
Out[3]:
['S2_FAPAR_V102_WEBMERCATOR2',
 'BIOPAR_LAI_V2_GLOBAL',
 'S2_SCENECLASSIFICATION',
 'BIOPAR_FAPAR_V1_GLOBAL',
 'PROBAV_L3_S10_TOC_NDVI_333M',
 'S2_MARKERMEER',
 'BIOPAR_FCOVER300_V1_GLOBAL',
 'BIOPAR_FCOVER_V1_GLOBAL',
 'BIOPAR_LAI300_V1_GLOBAL',
 'S2_FAPAR_CLOUDCOVER',
 'BIOPAR_FAPAR300_V1_GLOBAL',
 'BIOPAR_FAPAR_V2_GLOBAL',
 'BIOPAR_FCOVER_V2_GLOBAL',
 'S2_NDVI',
 'S2_LAI',
 'S2_FAPAR',
 'BIOPAR_LAI_V1_GLOBAL',
 'BIOPAR_NDVI300_V1_GLOBAL',
 'BIOPAR_NDVI_V2_GLOBAL',
 'S2_FAPAR_CLOUDCOVER_SC',
 'BIOPAR_WB300_V1_GLOBAL',
 'S2_FCOVER',
 'BIOPAR_BA300_V1_GLOBAL']

All of these layers can be queried, so let's go multitemporal, and query a PROBA-V and Sentinel 2 timeseries for the same area.

In [4]:
start = "2018-03-01"
end = "2018-12-01"
import dataclient
sentinel2_data = dataclient.get_timeseries(test_area,"S2_NDVI",start,end)
probav_data = dataclient.get_timeseries(test_area,"PROBAV_L3_S10_TOC_NDVI_333M",start,end)

The result is a Pandas Series object, which can be plotted immediately:

In [5]:
%matplotlib inline
import matplotlib.pyplot as plt
sentinel2_data.name="S2 NDVI"
probav_data.name = "PROBA-V NDVI"
sentinel2_data.dropna().plot(figsize=(12,4))
probav_data.dropna().plot(figsize=(12,4))
plt.legend()
plt.show()

The Sentinel 2 data is rather noisy, even though cloud filtering has occured. This cloud filtering is based on the Sen2Cor scene classification, which does not manage to remove all noise. (Mainly undetected cloud shadow.) One classic solution is to apply smoothing, this is also possible through a web service, that currently implements a Whittaker smoothing algorithm. The PROBA-V data is less noisy, because it is based on a 10-daily composite. This means that each data point is a combination of 10 daily images, where only the best value is retained.

In [6]:
import dataclient
sentinel2_smooth = dataclient.get_timeseries(test_area,"S2_NDVI",start,end,endpoint="http://pep.vgt.vito.be:8080/TimeseriesSmooth/v1.0/ts/")
probav_smooth = dataclient.get_timeseries(test_area,"PROBAV_L3_S10_TOC_NDVI_333M",start,end,endpoint="http://pep.vgt.vito.be:8080/TimeseriesSmooth/v1.0/ts/")
In [7]:
sentinel2_smooth.name="S2 NDVI Smooth"
probav_smooth.name = "PROBA-V NDVI Smooth"
sentinel2_smooth.dropna().plot(figsize=(10,4))
probav_smooth.dropna().plot(figsize=(10,4))
plt.legend()
plt.show()

As we can see, the smoothing clearly reduces the effect of the clouds on the Sentinel 2 data.

Timeseries dissection

Up to this point, we have only used the average pixel value for our timeseries analysis. However, we may want to derive more detailed statistics such as percentiles. To help you out, we added support for histograms to our web service and Python client.

The result of the 'dataclient.get_histogram' function is a time series of histograms, represented by a Pandas Dataframe.

In [8]:
import dataclient
avg_timeseries = dataclient.get_timeseries(test_area,"S2_NDVI", start, end)
histogram = dataclient.get_histogram(test_area,"S2_NDVI", start, end)
/opt/rh/rh-python35/root/usr/lib64/python3.5/site-packages/pandas/core/indexes/api.py:57: RuntimeWarning: unorderable types: str() < float(), sort order is undefined for incomparable objects
  union = _union_indexes(indexes)
/opt/rh/rh-python35/root/usr/lib64/python3.5/site-packages/pandas/core/indexes/api.py:57: RuntimeWarning: unorderable types: float() < str(), sort order is undefined for incomparable objects
  union = _union_indexes(indexes)
/opt/rh/rh-python35/root/usr/lib64/python3.5/site-packages/pandas/core/indexes/api.py:87: RuntimeWarning: unorderable types: float() < str(), sort order is undefined for incomparable objects
  result = result.union(other)
/opt/rh/rh-python35/root/usr/lib64/python3.5/site-packages/pandas/core/indexes/api.py:87: RuntimeWarning: unorderable types: str() < float(), sort order is undefined for incomparable objects
  result = result.union(other)
In [9]:
histogram = histogram.drop('NaN',axis=1)
histogram.columns = pd.to_numeric(histogram.columns)
histogram
Out[9]:
-0.08 -0.076 -0.068 -0.064 -0.06 -0.056 -0.052000000000000005 -0.048 -0.044 -0.04 ... 0.884 0.888 0.892 0.896 0.9 0.904 0.908 0.912 0.916 0.92
2018-03-02 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
2018-03-04 20.0 2.0 16.0 6.0 19.0 34.0 15.0 23.0 41.0 16.0 ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
2018-03-07 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
2018-03-09 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
2018-03-12 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
2018-03-14 NaN NaN NaN NaN 6.0 6.0 NaN NaN 6.0 NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
2018-03-17 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
2018-03-19 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
2018-03-22 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
2018-03-24 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
2018-03-27 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
2018-03-29 61.0 NaN 2.0 NaN 9.0 NaN 4.0 NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
2018-04-01 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
2018-04-03 9.0 NaN NaN NaN NaN NaN NaN 3.0 2.0 2.0 ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
2018-04-06 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
2018-04-08 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
2018-04-11 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
2018-04-13 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
2018-04-16 69.0 NaN 11.0 NaN 13.0 NaN 6.0 9.0 3.0 NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
2018-04-18 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
2018-04-21 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
2018-04-23 154.0 3.0 24.0 4.0 NaN 11.0 9.0 19.0 3.0 3.0 ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
2018-04-26 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
2018-04-28 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
2018-05-01 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
2018-05-03 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
2018-05-06 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ... 1463.0 1663.0 1561.0 1648.0 1401.0 1290.0 1231.0 871.0 451.0 233.0
2018-05-08 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ... 1406.0 1830.0 1596.0 1789.0 1620.0 1378.0 1191.0 892.0 580.0 350.0
2018-05-11 NaN NaN NaN NaN NaN 3.0 NaN 2.0 6.0 NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
2018-05-13 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
2018-08-19 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
2018-08-21 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
2018-08-24 167.0 5.0 9.0 8.0 20.0 34.0 19.0 6.0 43.0 18.0 ... 12.0 NaN NaN NaN NaN NaN NaN NaN NaN NaN
2018-08-26 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
2018-08-29 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
2018-08-31 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
2018-09-03 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
2018-09-05 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
2018-09-08 75.0 NaN NaN 3.0 NaN 3.0 NaN NaN 3.0 NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
2018-09-10 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
2018-09-13 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
2018-09-18 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ... 447.0 249.0 170.0 86.0 41.0 9.0 NaN NaN NaN NaN
2018-09-20 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
2018-09-23 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
2018-09-25 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
2018-09-28 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
2018-09-30 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ... 685.0 502.0 358.0 294.0 134.0 65.0 37.0 3.0 NaN NaN
2018-10-10 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ... 22.0 NaN 3.0 NaN NaN NaN NaN NaN NaN NaN
2018-10-18 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
2018-10-20 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
2018-10-23 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
2018-10-25 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
2018-10-30 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
2018-11-04 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
2018-11-09 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
2018-11-14 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
2018-11-19 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
2018-11-22 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
2018-11-24 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
2018-11-29 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN

98 rows × 251 columns

Based on these histograms, we can derive statistics that can also be plotted:

In [10]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
def quantile(row,quantile):
    cumsum = row.dropna().cumsum()
    #filter out rows with too little data, these are very often contaminated with remaining cloud or shadow
    if len(row.dropna())/len(row) < 0.4:
        return np.NaN
    return cumsum.index[np.searchsorted(cumsum,cumsum.max()*quantile)[0]]

    
median = histogram.apply(lambda row:quantile(row,0.5),axis=1)
high = histogram.apply(lambda row:quantile(row,0.95),axis=1)
low = histogram.apply(lambda row:quantile(row,0.05),axis=1)

q1 = histogram.apply(lambda row:quantile(row,0.25),axis=1).dropna()
q2 = histogram.apply(lambda row:quantile(row,0.75),axis=1).dropna()

fig, axes = plt.subplots(1, 1, figsize=(16,7))

plt.fill_between(q1.index,q2,q1,color="gray",alpha=0.5)
plt.plot(median.dropna(),label="median")
#plt.plot(high.dropna(),label="0.95 percentile")
#plt.plot(low.dropna(),label="0.05 percentile")

plt.plot(avg_timeseries.filter(median.dropna().index),label="average", linestyle="--")
plt.title("S2 NDVI")
plt.legend()
plt.show()

Timeseries timelapse

A final way to visualize the timeseries in full detail, is to retrieve the actual images, so we can plot them in an animated timelapse. To access the web service (WCS) that allows us to download these images, we need to provide our username and password. (The same as used for https://www.vito-eodata.be or https://terrascope.be )

In [ ]:
import getpass
import ipywidgets as widgets

username = widgets.Text(
    value='',
    placeholder='',
    description='Username:',
    disabled=False
)
username
In [ ]:
password = getpass.getpass()
credentials=(username.value, password)

With the following call, we download a set of images based on the bounding box of out test area. They are saved into the specified output directory. ('de_maat_rgb') We try to download an image for each available date. For some dates the image may be fully clouded, and for other dates the bounding box may not intersect the available data, which will result in an error that can be ignored.

In [ ]:
from dateutil.parser import parse
import pandas as pd
import dataclient
dataclient.MepLayer(credentials, "CGS__CGS_S2_10M_BANDS", "CGS", "https://sentineldata.vito.be/ows").download_bbox("de_maat_rgb_2018",test_area.geometry,startdate=start,enddate=end,verbose=False)

To filter our retrieved images by cloudcover, we retrieve a timeseries that contains a cloud percentage.

In [13]:
clouds = dataclient.get_timeseries(test_area,"S2_FAPAR_CLOUDCOVER",start,end)

The retrieved images contain 4 10m bands, with reflectance values from 0 to 10000, we define a function that plots this as RGB.

In [14]:
from rasterio.plot import show
def show_sentinel2_rgb(dataset,axis=None):
    band1 = np.clip((dataset.read(3,masked=True)-200)/1600.,0.0,1.0)
    band2 = np.clip((dataset.read(2,masked=True)-200)/1600.,0.0,1.0)
    band3 = np.clip((dataset.read(1,masked=True)-200)/1600.,0.0,1.0)
    
    rgb = np.dstack([band1,band2,band3])
    return axis.imshow(rgb,animated=True,vmin=0.0,vmax=1.0)
In [15]:
import pandas as pd
date_selection=["2018-05-06 00:00:00","2018-05-11 00:00:00","2018-05-26 00:00:00","2018-06-23 00:00:00"]
pd.to_datetime(date_selection)
Out[15]:
DatetimeIndex(['2018-05-06', '2018-05-11', '2018-05-26', '2018-06-23'], dtype='datetime64[ns]', freq=None)

Now we can convert this into a timelapse, which shows the image timeseries combined with the average timeseries.

In [26]:
%matplotlib inline

import rasterio
from IPython.display import HTML
from os import listdir
from os.path import isfile, join
import re
import matplotlib.animation as anim
import matplotlib.pyplot as plt
import numpy as np

#https://tomroelandts.com/articles/how-to-create-animated-gifs-with-python

image_paths = [f for f in listdir("de_maat_rgb_2018") if f.endswith("tiff") and isfile(join("de_maat_rgb_2018", f))]
image_paths.sort()

size=(800, 300)
fig = plt.figure(figsize=(14, 7))
#fig.set_size_inches(size[0] / 100, size[1] / 100)
ax1=plt.subplot2grid((4, 1), (0, 0), rowspan=3)
ax2=plt.subplot2grid((4, 1), (3, 0))
ax1.set_xticks([])
ax1.set_yticks([])

images = []
prog = re.compile(".*([0-9]{6}).*")
i = 0
for f in image_paths:
    with rasterio.open(join("de_maat_rgb_2018", f)) as dataset:
        fapar = dataset.read(1,masked=True)
        #print(fapar.count())
        if fapar.count() > 10000:                        
            date = pd.to_datetime(prog.match(f).group(1),yearfirst=True)    
            #if date in pd.to_datetime(date_selection):
            if date in sentinel2_smooth.dropna().index and date in clouds and clouds[date]<0.1:                
                plt_txt = ax1.text(0, 0, date.strftime("%Y-%m-%d"), color='red',fontsize=12)
                plt_im = show_sentinel2_rgb(dataset,ax1)
                graph, = ax2.plot(sentinel2_smooth.dropna())
                #point, = ax2.plot(date, filtered[date], 'g*')
                line =ax2.axvline(x=date,animated=True)
                images.append([plt_im,graph,line,plt_txt])
            #else:
            #    images.append([plt_im,graph])
            
animation = anim.ArtistAnimation(fig, images)
HTML(animation.to_jshtml())
Out[26]:


Once Loop Reflect