Metadata

title: "E-TRAINEE Tutorial - Raster time series processing in Python using xarray"
description: "This is a tutorial within the first theme of Module 1 of the E-TRAINEE course."
lastUpdate: 2023-07-20
authors: Andreas Mayr

Raster time series processing in Python using xarray¶

This tutorial provides an introduction to raster time series processing with Python using (mainly) the xarray package (Hoyer and Hamman 2017) for multi-dimensional data. Xarray's data structures (DataArray and Dataset, see here) complement NumPy-like arrays by labels in the form of dimensions, coordinates and attributes, thereby encoding information about how the array values are absolutely related to space, time etc. The labels can also be used to access and analyse the array values. These properties and functionality are very useful if you work with time series of remote sensing images or other raster data.

The Python packages used in this tutorial are contained in the requirements file provided for the course. Please see the instructions on the software page for setting up a Conda environment based on this file.

Data:

As a sample dataset, we use a time series of Sentinel-2 satellite data covering a relatively small area around the village of Obergurgl, located in the Central Alps of Tyrol (Austria), and the years 2018 to 2021. This imagery has been queried from the Google Earth Engine (GEE), then cloud-masked and exported from GEE as 6-band GeoTiffs, one GeoTiff per scene (see theme 2 for a tutorial on exporting cloudless Sentinel-2 imagery from GEE). You find this dataset in a folder named s2 in the data repository of the course.

In [23]:

import pathlib
import glob
import os
from datetime import datetime
import rioxarray
import xarray

import pathlib import glob import os from datetime import datetime import rioxarray import xarray

Loading many files of an image time series¶

We take the set of Sentinel-2 images stored locally as GeoTIFFs (after you downloaded them as the folder s2 from the course data repository). Using rioxarray, we load all these images to a DataArray and construct a DataSet from them. We follow more or less the approaches described here and here.

In [13]:

# Define the path to the s2 directory with the Sentinel-2 image subsets on your local system
in_dir = pathlib.Path('F:/data/etrainee/gee/s2')

# Obtain a list of all files to load from the s2 directory (defined previously as in_dir) on your local system
geotiff_list = glob.glob(str(in_dir / '**TPS.tif')) # We use 'TPS' to constrain the granules with better coverage in our AOI which is covered by TPT and TPS granules.

# Helper function to get the datetime index from the filenames
def paths_to_datetimeindex(paths):
    return [datetime.strptime((os.path.basename(i)[0:8]), '%Y%m%d') for i in paths]

# Create a variable used for the time dimension
time_var = xarray.Variable('time', paths_to_datetimeindex(geotiff_list))

# Load in and concatenate all individual GeoTIFFs
S2_da = xarray.concat([rioxarray.open_rasterio(i, masked=True) for i in geotiff_list], dim=time_var)

# Convert the xarray.DataArray into an xarray.Dataset
S2_ds = S2_da.to_dataset('band')

# Rename the variables to the more useful band names
band_names = {1: 'B2', 2: 'B3', 3: 'B4', 4: 'B8', 5: 'B11', 6: 'B12'}     # A dictionary with the names of the Sentinel-2 bands exported earlier.
S2_ds = S2_ds.rename(band_names)

# Define the path to the s2 directory with the Sentinel-2 image subsets on your local system in_dir = pathlib.Path('F:/data/etrainee/gee/s2') # Obtain a list of all files to load from the s2 directory (defined previously as in_dir) on your local system geotiff_list = glob.glob(str(in_dir / '**TPS.tif')) # We use 'TPS' to constrain the granules with better coverage in our AOI which is covered by TPT and TPS granules. # Helper function to get the datetime index from the filenames def paths_to_datetimeindex(paths): return [datetime.strptime((os.path.basename(i)[0:8]), '%Y%m%d') for i in paths] # Create a variable used for the time dimension time_var = xarray.Variable('time', paths_to_datetimeindex(geotiff_list)) # Load in and concatenate all individual GeoTIFFs S2_da = xarray.concat([rioxarray.open_rasterio(i, masked=True) for i in geotiff_list], dim=time_var) # Convert the xarray.DataArray into an xarray.Dataset S2_ds = S2_da.to_dataset('band') # Rename the variables to the more useful band names band_names = {1: 'B2', 2: 'B3', 3: 'B4', 4: 'B8', 5: 'B11', 6: 'B12'} # A dictionary with the names of the Sentinel-2 bands exported earlier. S2_ds = S2_ds.rename(band_names)

Now S2_ds is a DataSet with three dimensions (x, y, time) and six variables (the spectral bands). We add some attributes and fill them with metadata retrieved by the rio accessor. See the result:

In [14]:

# Fill in some of the missing metadata
S2_ds.attrs["crs"] = S2_ds.rio.crs
S2_ds.attrs["resolution"] = S2_ds.rio.resolution()
S2_ds.attrs["units"] = "meters"

S2_ds

# Fill in some of the missing metadata S2_ds.attrs["crs"] = S2_ds.rio.crs S2_ds.attrs["resolution"] = S2_ds.rio.resolution() S2_ds.attrs["units"] = "meters" S2_ds

Out[14]:

<xarray.Dataset>
Dimensions:      (time: 335, y: 201, x: 251)
Coordinates:
  * x            (x) float64 6.516e+05 6.516e+05 ... 6.566e+05 6.566e+05
  * y            (y) float64 5.194e+06 5.194e+06 5.194e+06 ... 5.19e+06 5.19e+06
    spatial_ref  int32 0
  * time         (time) datetime64[ns] 2018-01-02 2018-01-12 ... 2021-09-08
Data variables:
    B2           (time, y, x) float32 nan nan nan nan ... 1.095e+03 nan nan
    B3           (time, y, x) float32 nan nan nan nan ... 1.351e+03 nan nan
    B4           (time, y, x) float32 nan nan nan nan ... 1.445e+03 nan nan
    B8           (time, y, x) float32 nan nan nan nan ... 2.208e+03 nan nan
    B11          (time, y, x) float32 nan nan nan nan ... 2.162e+03 nan nan
    B12          (time, y, x) float32 nan nan nan nan ... 1.658e+03 nan nan
Attributes:
    AREA_OR_POINT:           Area
    TIFFTAG_RESOLUTIONUNIT:  1 (unitless)
    TIFFTAG_XRESOLUTION:     1
    TIFFTAG_YRESOLUTION:     1
    scale_factor:            1.0
    add_offset:              0.0
    crs:                     EPSG:32632
    resolution:              (20.0, -20.0)
    units:                   meters

xarray.Dataset

• Dimensions:
  • time: 335
  • y: 201
  • x: 251
• Coordinates: (4)
  • x
    (x)
    float64
    6.516e+05 6.516e+05 ... 6.566e+05

    array([651590., 651610., 651630., ..., 656550., 656570., 656590.])

  • y
    (y)
    float64
    5.194e+06 5.194e+06 ... 5.19e+06

    array([5193750., 5193730., 5193710., ..., 5189790., 5189770., 5189750.])

  • spatial_ref
    ()
    int32
    0

    crs_wkt :

    PROJCS["WGS 84 / UTM zone 32N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",9],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32632"]]

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314245179

    inverse_flattening :

    298.257223563

    reference_ellipsoid_name :

    WGS 84

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    WGS 84

    horizontal_datum_name :

    World Geodetic System 1984

    projected_crs_name :

    WGS 84 / UTM zone 32N

    grid_mapping_name :

    transverse_mercator

    latitude_of_projection_origin :

    0.0

    longitude_of_central_meridian :

    9.0

    false_easting :

    500000.0

    false_northing :

    0.0

    scale_factor_at_central_meridian :

    0.9996

    spatial_ref :

    PROJCS["WGS 84 / UTM zone 32N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",9],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32632"]]

    GeoTransform :

    651580.0 20.0 0.0 5193760.0 0.0 -20.0

    array(0)

  • time
    (time)
    datetime64[ns]
    2018-01-02 ... 2021-09-08

    array(['2018-01-02T00:00:00.000000000', '2018-01-12T00:00:00.000000000',
           '2018-01-27T00:00:00.000000000', ..., '2021-09-03T00:00:00.000000000',
           '2021-09-05T00:00:00.000000000', '2021-09-08T00:00:00.000000000'],
          dtype='datetime64[ns]')

• Data variables: (6)
  • B2
    (time, y, x)
    float32
    nan nan nan ... 1.095e+03 nan nan

    array([[[   nan,    nan,    nan, ...,    nan,    nan,    nan],
            [   nan,    nan,    nan, ...,    nan,    nan,    nan],
            [   nan,    nan,    nan, ...,    nan,    nan,    nan],
            ...,
            [   nan,    nan,    nan, ...,    nan,    nan,    nan],
            [   nan,    nan,    nan, ...,    nan,    nan,    nan],
            [   nan,    nan,    nan, ...,    nan,    nan,    nan]],
    
           [[   nan,    nan,    nan, ...,  7116.,  6893.,  6975.],
            [   nan,    nan,    nan, ...,  7284.,  6826.,  6941.],
            [   nan,    nan,    nan, ...,  7578.,  7313.,  7179.],
            ...,
            [10715., 10118., 10227., ...,  6981.,  6980.,  5471.],
            [ 9605.,  9503., 10095., ...,  8117.,  8512.,  5604.],
            [ 8790.,  9014.,  9334., ...,  8695.,  8726.,  7408.]],
    
           [[ 8002.,  8192.,  8645., ...,  7557.,  7214.,  7777.],
            [ 9680.,  8141.,  7723., ...,  7704.,  7163.,  7302.],
            [10009.,  9963.,  8712., ...,  8089.,  7683.,  7587.],
            ...,
    ...
            [  497.,   486.,   464., ...,   774.,   913.,   776.],
            [  502.,   454.,   415., ...,   781.,   927.,   694.],
            [  479.,   476.,   412., ...,   823.,   667.,   559.]],
    
           [[  688.,   768.,   801., ...,  1071.,  1124.,   811.],
            [  794.,   757.,   679., ...,   929.,  1109.,  1225.],
            [  838.,   751.,   724., ...,   847.,  1017.,  1118.],
            ...,
            [  512.,   496.,   487., ...,   821.,   905.,   817.],
            [  512.,   463.,   439., ...,   818.,   916.,   734.],
            [  512.,   480.,   430., ...,   832.,   745.,   647.]],
    
           [[   nan,    nan,    nan, ...,  1267.,  1419.,  1130.],
            [   nan,    nan,    nan, ...,  1190.,  1335.,  1522.],
            [   nan,    nan,    nan, ...,  1149.,  1255.,  1380.],
            ...,
            [  554.,   557.,   544., ...,  1122.,    nan,    nan],
            [  583.,   503.,   486., ...,  1059.,    nan,    nan],
            [  532.,   525.,   473., ...,  1095.,    nan,    nan]]],
          dtype=float32)

  • B3
    (time, y, x)
    float32
    nan nan nan ... 1.351e+03 nan nan

    array([[[   nan,    nan,    nan, ...,    nan,    nan,    nan],
            [   nan,    nan,    nan, ...,    nan,    nan,    nan],
            [   nan,    nan,    nan, ...,    nan,    nan,    nan],
            ...,
            [   nan,    nan,    nan, ...,    nan,    nan,    nan],
            [   nan,    nan,    nan, ...,    nan,    nan,    nan],
            [   nan,    nan,    nan, ...,    nan,    nan,    nan]],
    
           [[   nan,    nan,    nan, ...,  5978.,  5810.,  6208.],
            [   nan,    nan,    nan, ...,  6139.,  5621.,  5896.],
            [   nan,    nan,    nan, ...,  6440.,  6103.,  5997.],
            ...,
            [10655., 10096., 10182., ...,  6781.,  6842.,  5474.],
            [ 9626.,  9502., 10141., ...,  7924.,  8409.,  5659.],
            [ 8782.,  8944.,  9354., ...,  8530.,  8645.,  7237.]],
    
           [[ 8117.,  8269.,  8972., ...,  6509.,  6257.,  7102.],
            [ 9813.,  8405.,  7894., ...,  6683.,  6100.,  6382.],
            [10242.,  9839.,  8840., ...,  7019.,  6668.,  6504.],
            ...,
    ...
            [  795.,   776.,   790., ...,  1084.,  1220.,  1073.],
            [  830.,   751.,   638., ...,  1094.,  1223.,   977.],
            [  816.,   784.,   639., ...,  1148.,  1026.,   891.]],
    
           [[  895.,  1016.,  1046., ...,  1476.,  1504.,   973.],
            [ 1051.,   999.,   897., ...,  1318.,  1572.,  1684.],
            [ 1075.,  1032.,  1004., ...,  1093.,  1335.,  1499.],
            ...,
            [  792.,   756.,   788., ...,  1145.,  1242.,  1100.],
            [  832.,   726.,   664., ...,  1126.,  1255.,   993.],
            [  813.,   761.,   608., ...,  1166.,  1110.,  1008.]],
    
           [[   nan,    nan,    nan, ...,  1566.,  1669.,  1264.],
            [   nan,    nan,    nan, ...,  1447.,  1659.,  1835.],
            [   nan,    nan,    nan, ...,  1241.,  1465.,  1669.],
            ...,
            [  817.,   804.,   824., ...,  1349.,    nan,    nan],
            [  848.,   767.,   687., ...,  1326.,    nan,    nan],
            [  836.,   785.,   654., ...,  1351.,    nan,    nan]]],
          dtype=float32)

  • B4
    (time, y, x)
    float32
    nan nan nan ... 1.445e+03 nan nan

    array([[[   nan,    nan,    nan, ...,    nan,    nan,    nan],
            [   nan,    nan,    nan, ...,    nan,    nan,    nan],
            [   nan,    nan,    nan, ...,    nan,    nan,    nan],
            ...,
            [   nan,    nan,    nan, ...,    nan,    nan,    nan],
            [   nan,    nan,    nan, ...,    nan,    nan,    nan],
            [   nan,    nan,    nan, ...,    nan,    nan,    nan]],
    
           [[   nan,    nan,    nan, ...,  5187.,  4992.,  5774.],
            [   nan,    nan,    nan, ...,  5256.,  4782.,  5134.],
            [   nan,    nan,    nan, ...,  5649.,  5427.,  5273.],
            ...,
            [10998., 10420., 10512., ...,  6941.,  7017.,  5465.],
            [ 9933.,  9851., 10476., ...,  8076.,  8508.,  5783.],
            [ 9070.,  9271.,  9623., ...,  8672.,  8696.,  7421.]],
    
           [[ 8208.,  8373.,  8949., ...,  5550.,  5398.,  6620.],
            [ 9820.,  8390.,  8020., ...,  5718.,  5078.,  5425.],
            [10243.,  9897.,  8857., ...,  6095.,  5775.,  5661.],
            ...,
    ...
            [  775.,   769.,   764., ...,  1193.,  1364.,  1175.],
            [  790.,   726.,   679., ...,  1193.,  1369.,  1022.],
            [  738.,   762.,   646., ...,  1250.,  1056.,   955.]],
    
           [[ 1019.,  1111.,  1139., ...,  1410.,  1493.,  1009.],
            [ 1136.,  1104.,   991., ...,  1264.,  1481.,  1698.],
            [ 1210.,  1122.,  1142., ...,  1130.,  1273.,  1535.],
            ...,
            [  810.,   789.,   779., ...,  1277.,  1401.,  1247.],
            [  800.,   740.,   700., ...,  1233.,  1391.,  1065.],
            [  776.,   774.,   651., ...,  1275.,  1190.,  1089.]],
    
           [[   nan,    nan,    nan, ...,  1547.,  1639.,  1208.],
            [   nan,    nan,    nan, ...,  1355.,  1571.,  1820.],
            [   nan,    nan,    nan, ...,  1130.,  1404.,  1637.],
            ...,
            [  840.,   839.,   814., ...,  1475.,    nan,    nan],
            [  855.,   778.,   742., ...,  1444.,    nan,    nan],
            [  802.,   836.,   683., ...,  1445.,    nan,    nan]]],
          dtype=float32)

  • B8
    (time, y, x)
    float32
    nan nan nan ... 2.208e+03 nan nan

    array([[[   nan,    nan,    nan, ...,    nan,    nan,    nan],
            [   nan,    nan,    nan, ...,    nan,    nan,    nan],
            [   nan,    nan,    nan, ...,    nan,    nan,    nan],
            ...,
            [   nan,    nan,    nan, ...,    nan,    nan,    nan],
            [   nan,    nan,    nan, ...,    nan,    nan,    nan],
            [   nan,    nan,    nan, ...,    nan,    nan,    nan]],
    
           [[   nan,    nan,    nan, ...,  4276.,  4089.,  5000.],
            [   nan,    nan,    nan, ...,  4322.,  3919.,  4197.],
            [   nan,    nan,    nan, ...,  4684.,  4451.,  4310.],
            ...,
            [10667., 10119., 10220., ...,  6789.,  6765.,  5422.],
            [ 9642.,  9625., 10164., ...,  7770.,  8218.,  5703.],
            [ 8856.,  9046.,  9407., ...,  8321.,  8372.,  7258.]],
    
           [[ 7688.,  7870.,  8241., ...,  4393.,  4298.,  5567.],
            [ 8961.,  7822.,  7538., ...,  4488.,  4061.,  4483.],
            [ 9254.,  9074.,  8255., ...,  4840.,  4612.,  4475.],
            ...,
    ...
            [ 3299.,  3043.,  3446., ...,  2016.,  1802.,  1492.],
            [ 3497.,  3196.,  2681., ...,  2136.,  1993.,  2335.],
            [ 3557.,  3366.,  2658., ...,  2082.,  2654.,  2835.]],
    
           [[ 1535.,  1779.,  1725., ...,  5187.,  5241.,  3208.],
            [ 1689.,  1687.,  1717., ...,  4914.,  5625.,  5237.],
            [ 1694.,  1749.,  1784., ...,  3391.,  4121.,  4569.],
            ...,
            [ 3162.,  2974.,  3239., ...,  2032.,  1848.,  1540.],
            [ 3349.,  3044.,  2672., ...,  2234.,  1945.,  2323.],
            [ 3347.,  3170.,  2343., ...,  2111.,  2656.,  3014.]],
    
           [[   nan,    nan,    nan, ...,  4832.,  5022.,  3458.],
            [   nan,    nan,    nan, ...,  4571.,  5241.,  5036.],
            [   nan,    nan,    nan, ...,  3000.,  3851.,  4310.],
            ...,
            [ 3175.,  3096.,  3361., ...,  2087.,    nan,    nan],
            [ 3332.,  3154.,  2701., ...,  2206.,    nan,    nan],
            [ 3405.,  3286.,  2539., ...,  2208.,    nan,    nan]]],
          dtype=float32)

  • B11
    (time, y, x)
    float32
    nan nan nan ... 2.162e+03 nan nan

    array([[[  nan,   nan,   nan, ...,   nan,   nan,   nan],
            [  nan,   nan,   nan, ...,   nan,   nan,   nan],
            [  nan,   nan,   nan, ...,   nan,   nan,   nan],
            ...,
            [  nan,   nan,   nan, ...,   nan,   nan,   nan],
            [  nan,   nan,   nan, ...,   nan,   nan,   nan],
            [  nan,   nan,   nan, ...,   nan,   nan,   nan]],
    
           [[  nan,   nan,   nan, ...,   35.,   50.,   86.],
            [  nan,   nan,   nan, ...,   94.,   54.,   61.],
            [  nan,   nan,   nan, ...,  101.,  126.,  130.],
            ...,
            [1356., 1331., 1330., ..., 1038., 1095., 1011.],
            [1263., 1274., 1333., ..., 1047., 1077., 1035.],
            [1193., 1209., 1275., ..., 1077., 1092., 1142.]],
    
           [[ 824.,  810.,  796., ...,   40.,   44.,   69.],
            [ 840.,  803.,  812., ...,   32.,   40.,   65.],
            [ 829.,  822.,  827., ...,   12.,   45.,   53.],
            ...,
    ...
            ...,
            [2406., 2307., 2287., ..., 1990., 1866., 1749.],
            [2454., 2242., 2079., ..., 1934., 1962., 2062.],
            [2454., 2350., 2060., ..., 2048., 2293., 2423.]],
    
           [[2289., 2594., 2779., ..., 4196., 4083., 3178.],
            [2516., 2651., 2776., ..., 3809., 4438., 4386.],
            [2518., 2746., 2777., ..., 3333., 3912., 4490.],
            ...,
            [2509., 2390., 2469., ..., 2143., 2002., 1930.],
            [2540., 2280., 2218., ..., 2073., 1990., 2036.],
            [2640., 2424., 2123., ..., 2080., 2351., 2600.]],
    
           [[  nan,   nan,   nan, ..., 3937., 3844., 3249.],
            [  nan,   nan,   nan, ..., 3667., 4190., 4182.],
            [  nan,   nan,   nan, ..., 2992., 3647., 4200.],
            ...,
            [2477., 2393., 2369., ..., 2132.,   nan,   nan],
            [2496., 2326., 2216., ..., 2071.,   nan,   nan],
            [2542., 2395., 2168., ..., 2162.,   nan,   nan]]], dtype=float32)

  • B12
    (time, y, x)
    float32
    nan nan nan ... 1.658e+03 nan nan

    array([[[  nan,   nan,   nan, ...,   nan,   nan,   nan],
            [  nan,   nan,   nan, ...,   nan,   nan,   nan],
            [  nan,   nan,   nan, ...,   nan,   nan,   nan],
            ...,
            [  nan,   nan,   nan, ...,   nan,   nan,   nan],
            [  nan,   nan,   nan, ...,   nan,   nan,   nan],
            [  nan,   nan,   nan, ...,   nan,   nan,   nan]],
    
           [[  nan,   nan,   nan, ...,   62.,  142.,  100.],
            [  nan,   nan,   nan, ...,   89.,   62.,  100.],
            [  nan,   nan,   nan, ...,  100.,  100.,  145.],
            ...,
            [1512., 1467., 1478., ..., 1123., 1173., 1009.],
            [1411., 1414., 1479., ..., 1188., 1209., 1055.],
            [1322., 1324., 1419., ..., 1227., 1246., 1197.]],
    
           [[ 882.,  878.,  836., ...,  101.,  137.,  123.],
            [ 889.,  870.,  863., ...,   79.,  114.,  106.],
            [ 872.,  867.,  870., ...,   97.,   92.,  128.],
            ...,
    ...
            ...,
            [1350., 1284., 1255., ..., 1509., 1484., 1336.],
            [1372., 1245., 1159., ..., 1440., 1468., 1487.],
            [1390., 1290., 1127., ..., 1463., 1527., 1605.]],
    
           [[1857., 2110., 2247., ..., 2527., 2581., 2188.],
            [2008., 2116., 2225., ..., 2318., 2767., 2799.],
            [2028., 2187., 2198., ..., 2048., 2532., 2998.],
            ...,
            [1450., 1360., 1385., ..., 1683., 1617., 1533.],
            [1439., 1312., 1252., ..., 1542., 1560., 1582.],
            [1512., 1376., 1230., ..., 1562., 1655., 1764.]],
    
           [[  nan,   nan,   nan, ..., 2653., 2550., 2230.],
            [  nan,   nan,   nan, ..., 2424., 2809., 2842.],
            [  nan,   nan,   nan, ..., 2098., 2443., 2885.],
            ...,
            [1460., 1415., 1361., ..., 1757.,   nan,   nan],
            [1483., 1355., 1289., ..., 1645.,   nan,   nan],
            [1480., 1364., 1264., ..., 1658.,   nan,   nan]]], dtype=float32)

• Indexes: (3)
  • x
    PandasIndex

    PandasIndex(Float64Index([651590.0, 651610.0, 651630.0, 651650.0, 651670.0, 651690.0,
                  651710.0, 651730.0, 651750.0, 651770.0,
                  ...
                  656410.0, 656430.0, 656450.0, 656470.0, 656490.0, 656510.0,
                  656530.0, 656550.0, 656570.0, 656590.0],
                 dtype='float64', name='x', length=251))

  • y
    PandasIndex

    PandasIndex(Float64Index([5193750.0, 5193730.0, 5193710.0, 5193690.0, 5193670.0, 5193650.0,
                  5193630.0, 5193610.0, 5193590.0, 5193570.0,
                  ...
                  5189930.0, 5189910.0, 5189890.0, 5189870.0, 5189850.0, 5189830.0,
                  5189810.0, 5189790.0, 5189770.0, 5189750.0],
                 dtype='float64', name='y', length=201))

  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['2018-01-02', '2018-01-12', '2018-01-27', '2018-01-29',
                   '2018-02-08', '2018-02-11', '2018-02-13', '2018-02-16',
                   '2018-02-18', '2018-02-21',
                   ...
                   '2021-08-16', '2021-08-19', '2021-08-21', '2021-08-24',
                   '2021-08-26', '2021-08-31', '2021-09-03', '2021-09-03',
                   '2021-09-05', '2021-09-08'],
                  dtype='datetime64[ns]', name='time', length=335, freq=None))

• Attributes: (9)

  AREA_OR_POINT :

  Area

  TIFFTAG_RESOLUTIONUNIT :

  1 (unitless)

  TIFFTAG_XRESOLUTION :

  1

  TIFFTAG_YRESOLUTION :

  1

  scale_factor :

  1.0

  add_offset :

  0.0

  crs :

  EPSG:32632

  resolution :

  (20.0, -20.0)

  units :

  meters

A spatial_ref coordinate with length one was automatically created but is not needed. Let's get rid of this.

In [15]:

S2_ds = S2_ds.drop_vars(names='spatial_ref')    # Drop the variable with specified name
S2_ds

S2_ds = S2_ds.drop_vars(names='spatial_ref') # Drop the variable with specified name S2_ds

Out[15]:

<xarray.Dataset>
Dimensions:  (time: 335, y: 201, x: 251)
Coordinates:
  * x        (x) float64 6.516e+05 6.516e+05 6.516e+05 ... 6.566e+05 6.566e+05
  * y        (y) float64 5.194e+06 5.194e+06 5.194e+06 ... 5.19e+06 5.19e+06
  * time     (time) datetime64[ns] 2018-01-02 2018-01-12 ... 2021-09-08
Data variables:
    B2       (time, y, x) float32 nan nan nan nan ... 1.095e+03 nan nan
    B3       (time, y, x) float32 nan nan nan nan ... 1.36e+03 1.351e+03 nan nan
    B4       (time, y, x) float32 nan nan nan nan ... 1.445e+03 nan nan
    B8       (time, y, x) float32 nan nan nan nan ... 2.208e+03 nan nan
    B11      (time, y, x) float32 nan nan nan nan ... 1.98e+03 2.162e+03 nan nan
    B12      (time, y, x) float32 nan nan nan nan ... 1.56e+03 1.658e+03 nan nan
Attributes:
    AREA_OR_POINT:           Area
    TIFFTAG_RESOLUTIONUNIT:  1 (unitless)
    TIFFTAG_XRESOLUTION:     1
    TIFFTAG_YRESOLUTION:     1
    scale_factor:            1.0
    add_offset:              0.0
    crs:                     EPSG:32632
    resolution:              (20.0, -20.0)
    units:                   meters

xarray.Dataset

• Dimensions:
  • time: 335
  • y: 201
  • x: 251
• Coordinates: (3)
  • x
    (x)
    float64
    6.516e+05 6.516e+05 ... 6.566e+05

    array([651590., 651610., 651630., ..., 656550., 656570., 656590.])

  • y
    (y)
    float64
    5.194e+06 5.194e+06 ... 5.19e+06

    array([5193750., 5193730., 5193710., ..., 5189790., 5189770., 5189750.])

  • time
    (time)
    datetime64[ns]
    2018-01-02 ... 2021-09-08

    array(['2018-01-02T00:00:00.000000000', '2018-01-12T00:00:00.000000000',
           '2018-01-27T00:00:00.000000000', ..., '2021-09-03T00:00:00.000000000',
           '2021-09-05T00:00:00.000000000', '2021-09-08T00:00:00.000000000'],
          dtype='datetime64[ns]')

• Data variables: (6)
  • B2
    (time, y, x)
    float32
    nan nan nan ... 1.095e+03 nan nan

    array([[[   nan,    nan,    nan, ...,    nan,    nan,    nan],
            [   nan,    nan,    nan, ...,    nan,    nan,    nan],
            [   nan,    nan,    nan, ...,    nan,    nan,    nan],
            ...,
            [   nan,    nan,    nan, ...,    nan,    nan,    nan],
            [   nan,    nan,    nan, ...,    nan,    nan,    nan],
            [   nan,    nan,    nan, ...,    nan,    nan,    nan]],
    
           [[   nan,    nan,    nan, ...,  7116.,  6893.,  6975.],
            [   nan,    nan,    nan, ...,  7284.,  6826.,  6941.],
            [   nan,    nan,    nan, ...,  7578.,  7313.,  7179.],
            ...,
            [10715., 10118., 10227., ...,  6981.,  6980.,  5471.],
            [ 9605.,  9503., 10095., ...,  8117.,  8512.,  5604.],
            [ 8790.,  9014.,  9334., ...,  8695.,  8726.,  7408.]],
    
           [[ 8002.,  8192.,  8645., ...,  7557.,  7214.,  7777.],
            [ 9680.,  8141.,  7723., ...,  7704.,  7163.,  7302.],
            [10009.,  9963.,  8712., ...,  8089.,  7683.,  7587.],
            ...,
    ...
            [  497.,   486.,   464., ...,   774.,   913.,   776.],
            [  502.,   454.,   415., ...,   781.,   927.,   694.],
            [  479.,   476.,   412., ...,   823.,   667.,   559.]],
    
           [[  688.,   768.,   801., ...,  1071.,  1124.,   811.],
            [  794.,   757.,   679., ...,   929.,  1109.,  1225.],
            [  838.,   751.,   724., ...,   847.,  1017.,  1118.],
            ...,
            [  512.,   496.,   487., ...,   821.,   905.,   817.],
            [  512.,   463.,   439., ...,   818.,   916.,   734.],
            [  512.,   480.,   430., ...,   832.,   745.,   647.]],
    
           [[   nan,    nan,    nan, ...,  1267.,  1419.,  1130.],
            [   nan,    nan,    nan, ...,  1190.,  1335.,  1522.],
            [   nan,    nan,    nan, ...,  1149.,  1255.,  1380.],
            ...,
            [  554.,   557.,   544., ...,  1122.,    nan,    nan],
            [  583.,   503.,   486., ...,  1059.,    nan,    nan],
            [  532.,   525.,   473., ...,  1095.,    nan,    nan]]],
          dtype=float32)

  • B3
    (time, y, x)
    float32
    nan nan nan ... 1.351e+03 nan nan

    array([[[   nan,    nan,    nan, ...,    nan,    nan,    nan],
            [   nan,    nan,    nan, ...,    nan,    nan,    nan],
            [   nan,    nan,    nan, ...,    nan,    nan,    nan],
            ...,
            [   nan,    nan,    nan, ...,    nan,    nan,    nan],
            [   nan,    nan,    nan, ...,    nan,    nan,    nan],
            [   nan,    nan,    nan, ...,    nan,    nan,    nan]],
    
           [[   nan,    nan,    nan, ...,  5978.,  5810.,  6208.],
            [   nan,    nan,    nan, ...,  6139.,  5621.,  5896.],
            [   nan,    nan,    nan, ...,  6440.,  6103.,  5997.],
            ...,
            [10655., 10096., 10182., ...,  6781.,  6842.,  5474.],
            [ 9626.,  9502., 10141., ...,  7924.,  8409.,  5659.],
            [ 8782.,  8944.,  9354., ...,  8530.,  8645.,  7237.]],
    
           [[ 8117.,  8269.,  8972., ...,  6509.,  6257.,  7102.],
            [ 9813.,  8405.,  7894., ...,  6683.,  6100.,  6382.],
            [10242.,  9839.,  8840., ...,  7019.,  6668.,  6504.],
            ...,
    ...
            [  795.,   776.,   790., ...,  1084.,  1220.,  1073.],
            [  830.,   751.,   638., ...,  1094.,  1223.,   977.],
            [  816.,   784.,   639., ...,  1148.,  1026.,   891.]],
    
           [[  895.,  1016.,  1046., ...,  1476.,  1504.,   973.],
            [ 1051.,   999.,   897., ...,  1318.,  1572.,  1684.],
            [ 1075.,  1032.,  1004., ...,  1093.,  1335.,  1499.],
            ...,
            [  792.,   756.,   788., ...,  1145.,  1242.,  1100.],
            [  832.,   726.,   664., ...,  1126.,  1255.,   993.],
            [  813.,   761.,   608., ...,  1166.,  1110.,  1008.]],
    
           [[   nan,    nan,    nan, ...,  1566.,  1669.,  1264.],
            [   nan,    nan,    nan, ...,  1447.,  1659.,  1835.],
            [   nan,    nan,    nan, ...,  1241.,  1465.,  1669.],
            ...,
            [  817.,   804.,   824., ...,  1349.,    nan,    nan],
            [  848.,   767.,   687., ...,  1326.,    nan,    nan],
            [  836.,   785.,   654., ...,  1351.,    nan,    nan]]],
          dtype=float32)

  • B4
    (time, y, x)
    float32
    nan nan nan ... 1.445e+03 nan nan

    array([[[   nan,    nan,    nan, ...,    nan,    nan,    nan],
            [   nan,    nan,    nan, ...,    nan,    nan,    nan],
            [   nan,    nan,    nan, ...,    nan,    nan,    nan],
            ...,
            [   nan,    nan,    nan, ...,    nan,    nan,    nan],
            [   nan,    nan,    nan, ...,    nan,    nan,    nan],
            [   nan,    nan,    nan, ...,    nan,    nan,    nan]],
    
           [[   nan,    nan,    nan, ...,  5187.,  4992.,  5774.],
            [   nan,    nan,    nan, ...,  5256.,  4782.,  5134.],
            [   nan,    nan,    nan, ...,  5649.,  5427.,  5273.],
            ...,
            [10998., 10420., 10512., ...,  6941.,  7017.,  5465.],
            [ 9933.,  9851., 10476., ...,  8076.,  8508.,  5783.],
            [ 9070.,  9271.,  9623., ...,  8672.,  8696.,  7421.]],
    
           [[ 8208.,  8373.,  8949., ...,  5550.,  5398.,  6620.],
            [ 9820.,  8390.,  8020., ...,  5718.,  5078.,  5425.],
            [10243.,  9897.,  8857., ...,  6095.,  5775.,  5661.],
            ...,
    ...
            [  775.,   769.,   764., ...,  1193.,  1364.,  1175.],
            [  790.,   726.,   679., ...,  1193.,  1369.,  1022.],
            [  738.,   762.,   646., ...,  1250.,  1056.,   955.]],
    
           [[ 1019.,  1111.,  1139., ...,  1410.,  1493.,  1009.],
            [ 1136.,  1104.,   991., ...,  1264.,  1481.,  1698.],
            [ 1210.,  1122.,  1142., ...,  1130.,  1273.,  1535.],
            ...,
            [  810.,   789.,   779., ...,  1277.,  1401.,  1247.],
            [  800.,   740.,   700., ...,  1233.,  1391.,  1065.],
            [  776.,   774.,   651., ...,  1275.,  1190.,  1089.]],
    
           [[   nan,    nan,    nan, ...,  1547.,  1639.,  1208.],
            [   nan,    nan,    nan, ...,  1355.,  1571.,  1820.],
            [   nan,    nan,    nan, ...,  1130.,  1404.,  1637.],
            ...,
            [  840.,   839.,   814., ...,  1475.,    nan,    nan],
            [  855.,   778.,   742., ...,  1444.,    nan,    nan],
            [  802.,   836.,   683., ...,  1445.,    nan,    nan]]],
          dtype=float32)

  • B8
    (time, y, x)
    float32
    nan nan nan ... 2.208e+03 nan nan

    array([[[   nan,    nan,    nan, ...,    nan,    nan,    nan],
            [   nan,    nan,    nan, ...,    nan,    nan,    nan],
            [   nan,    nan,    nan, ...,    nan,    nan,    nan],
            ...,
            [   nan,    nan,    nan, ...,    nan,    nan,    nan],
            [   nan,    nan,    nan, ...,    nan,    nan,    nan],
            [   nan,    nan,    nan, ...,    nan,    nan,    nan]],
    
           [[   nan,    nan,    nan, ...,  4276.,  4089.,  5000.],
            [   nan,    nan,    nan, ...,  4322.,  3919.,  4197.],
            [   nan,    nan,    nan, ...,  4684.,  4451.,  4310.],
            ...,
            [10667., 10119., 10220., ...,  6789.,  6765.,  5422.],
            [ 9642.,  9625., 10164., ...,  7770.,  8218.,  5703.],
            [ 8856.,  9046.,  9407., ...,  8321.,  8372.,  7258.]],
    
           [[ 7688.,  7870.,  8241., ...,  4393.,  4298.,  5567.],
            [ 8961.,  7822.,  7538., ...,  4488.,  4061.,  4483.],
            [ 9254.,  9074.,  8255., ...,  4840.,  4612.,  4475.],
            ...,
    ...
            [ 3299.,  3043.,  3446., ...,  2016.,  1802.,  1492.],
            [ 3497.,  3196.,  2681., ...,  2136.,  1993.,  2335.],
            [ 3557.,  3366.,  2658., ...,  2082.,  2654.,  2835.]],
    
           [[ 1535.,  1779.,  1725., ...,  5187.,  5241.,  3208.],
            [ 1689.,  1687.,  1717., ...,  4914.,  5625.,  5237.],
            [ 1694.,  1749.,  1784., ...,  3391.,  4121.,  4569.],
            ...,
            [ 3162.,  2974.,  3239., ...,  2032.,  1848.,  1540.],
            [ 3349.,  3044.,  2672., ...,  2234.,  1945.,  2323.],
            [ 3347.,  3170.,  2343., ...,  2111.,  2656.,  3014.]],
    
           [[   nan,    nan,    nan, ...,  4832.,  5022.,  3458.],
            [   nan,    nan,    nan, ...,  4571.,  5241.,  5036.],
            [   nan,    nan,    nan, ...,  3000.,  3851.,  4310.],
            ...,
            [ 3175.,  3096.,  3361., ...,  2087.,    nan,    nan],
            [ 3332.,  3154.,  2701., ...,  2206.,    nan,    nan],
            [ 3405.,  3286.,  2539., ...,  2208.,    nan,    nan]]],
          dtype=float32)

  • B11
    (time, y, x)
    float32
    nan nan nan ... 2.162e+03 nan nan

    array([[[  nan,   nan,   nan, ...,   nan,   nan,   nan],
            [  nan,   nan,   nan, ...,   nan,   nan,   nan],
            [  nan,   nan,   nan, ...,   nan,   nan,   nan],
            ...,
            [  nan,   nan,   nan, ...,   nan,   nan,   nan],
            [  nan,   nan,   nan, ...,   nan,   nan,   nan],
            [  nan,   nan,   nan, ...,   nan,   nan,   nan]],
    
           [[  nan,   nan,   nan, ...,   35.,   50.,   86.],
            [  nan,   nan,   nan, ...,   94.,   54.,   61.],
            [  nan,   nan,   nan, ...,  101.,  126.,  130.],
            ...,
            [1356., 1331., 1330., ..., 1038., 1095., 1011.],
            [1263., 1274., 1333., ..., 1047., 1077., 1035.],
            [1193., 1209., 1275., ..., 1077., 1092., 1142.]],
    
           [[ 824.,  810.,  796., ...,   40.,   44.,   69.],
            [ 840.,  803.,  812., ...,   32.,   40.,   65.],
            [ 829.,  822.,  827., ...,   12.,   45.,   53.],
            ...,
    ...
            ...,
            [2406., 2307., 2287., ..., 1990., 1866., 1749.],
            [2454., 2242., 2079., ..., 1934., 1962., 2062.],
            [2454., 2350., 2060., ..., 2048., 2293., 2423.]],
    
           [[2289., 2594., 2779., ..., 4196., 4083., 3178.],
            [2516., 2651., 2776., ..., 3809., 4438., 4386.],
            [2518., 2746., 2777., ..., 3333., 3912., 4490.],
            ...,
            [2509., 2390., 2469., ..., 2143., 2002., 1930.],
            [2540., 2280., 2218., ..., 2073., 1990., 2036.],
            [2640., 2424., 2123., ..., 2080., 2351., 2600.]],
    
           [[  nan,   nan,   nan, ..., 3937., 3844., 3249.],
            [  nan,   nan,   nan, ..., 3667., 4190., 4182.],
            [  nan,   nan,   nan, ..., 2992., 3647., 4200.],
            ...,
            [2477., 2393., 2369., ..., 2132.,   nan,   nan],
            [2496., 2326., 2216., ..., 2071.,   nan,   nan],
            [2542., 2395., 2168., ..., 2162.,   nan,   nan]]], dtype=float32)

  • B12
    (time, y, x)
    float32
    nan nan nan ... 1.658e+03 nan nan

    array([[[  nan,   nan,   nan, ...,   nan,   nan,   nan],
            [  nan,   nan,   nan, ...,   nan,   nan,   nan],
            [  nan,   nan,   nan, ...,   nan,   nan,   nan],
            ...,
            [  nan,   nan,   nan, ...,   nan,   nan,   nan],
            [  nan,   nan,   nan, ...,   nan,   nan,   nan],
            [  nan,   nan,   nan, ...,   nan,   nan,   nan]],
    
           [[  nan,   nan,   nan, ...,   62.,  142.,  100.],
            [  nan,   nan,   nan, ...,   89.,   62.,  100.],
            [  nan,   nan,   nan, ...,  100.,  100.,  145.],
            ...,
            [1512., 1467., 1478., ..., 1123., 1173., 1009.],
            [1411., 1414., 1479., ..., 1188., 1209., 1055.],
            [1322., 1324., 1419., ..., 1227., 1246., 1197.]],
    
           [[ 882.,  878.,  836., ...,  101.,  137.,  123.],
            [ 889.,  870.,  863., ...,   79.,  114.,  106.],
            [ 872.,  867.,  870., ...,   97.,   92.,  128.],
            ...,
    ...
            ...,
            [1350., 1284., 1255., ..., 1509., 1484., 1336.],
            [1372., 1245., 1159., ..., 1440., 1468., 1487.],
            [1390., 1290., 1127., ..., 1463., 1527., 1605.]],
    
           [[1857., 2110., 2247., ..., 2527., 2581., 2188.],
            [2008., 2116., 2225., ..., 2318., 2767., 2799.],
            [2028., 2187., 2198., ..., 2048., 2532., 2998.],
            ...,
            [1450., 1360., 1385., ..., 1683., 1617., 1533.],
            [1439., 1312., 1252., ..., 1542., 1560., 1582.],
            [1512., 1376., 1230., ..., 1562., 1655., 1764.]],
    
           [[  nan,   nan,   nan, ..., 2653., 2550., 2230.],
            [  nan,   nan,   nan, ..., 2424., 2809., 2842.],
            [  nan,   nan,   nan, ..., 2098., 2443., 2885.],
            ...,
            [1460., 1415., 1361., ..., 1757.,   nan,   nan],
            [1483., 1355., 1289., ..., 1645.,   nan,   nan],
            [1480., 1364., 1264., ..., 1658.,   nan,   nan]]], dtype=float32)

• Indexes: (3)
  • x
    PandasIndex

    PandasIndex(Float64Index([651590.0, 651610.0, 651630.0, 651650.0, 651670.0, 651690.0,
                  651710.0, 651730.0, 651750.0, 651770.0,
                  ...
                  656410.0, 656430.0, 656450.0, 656470.0, 656490.0, 656510.0,
                  656530.0, 656550.0, 656570.0, 656590.0],
                 dtype='float64', name='x', length=251))

  • y
    PandasIndex

    PandasIndex(Float64Index([5193750.0, 5193730.0, 5193710.0, 5193690.0, 5193670.0, 5193650.0,
                  5193630.0, 5193610.0, 5193590.0, 5193570.0,
                  ...
                  5189930.0, 5189910.0, 5189890.0, 5189870.0, 5189850.0, 5189830.0,
                  5189810.0, 5189790.0, 5189770.0, 5189750.0],
                 dtype='float64', name='y', length=201))

  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['2018-01-02', '2018-01-12', '2018-01-27', '2018-01-29',
                   '2018-02-08', '2018-02-11', '2018-02-13', '2018-02-16',
                   '2018-02-18', '2018-02-21',
                   ...
                   '2021-08-16', '2021-08-19', '2021-08-21', '2021-08-24',
                   '2021-08-26', '2021-08-31', '2021-09-03', '2021-09-03',
                   '2021-09-05', '2021-09-08'],
                  dtype='datetime64[ns]', name='time', length=335, freq=None))

• Attributes: (9)

  AREA_OR_POINT :

  Area

  TIFFTAG_RESOLUTIONUNIT :

  1 (unitless)

  TIFFTAG_XRESOLUTION :

  1

  TIFFTAG_YRESOLUTION :

  1

  scale_factor :

  1.0

  add_offset :

  0.0

  crs :

  EPSG:32632

  resolution :

  (20.0, -20.0)

  units :

  meters

Working with images and time series¶

Plotting¶

As a next step we will explore some of the xarray plotting possibilities with our image time series. If you need more advanced functionality, try out more specialized packages (such as Cartopy for making maps and hvplot for interactive plots, see this tutorial).

First, select one scene and one band and plot this:

In [16]:

S2_ds.B2.sel(time='2020-01-02').plot()

S2_ds.B2.sel(time='2020-01-02').plot()

Out[16]:

<matplotlib.collections.QuadMesh at 0x27b2bcbd060>

Other possibilities for plotting include scatter plots, e.g. with two bands plotted against each other:

In [17]:

S2_ds.sel(time='2020-01-02').plot.scatter(x='B3', y='B4', marker='.', edgecolor='none', alpha=0.01)

S2_ds.sel(time='2020-01-02').plot.scatter(x='B3', y='B4', marker='.', edgecolor='none', alpha=0.01)

Out[17]:

<matplotlib.collections.PathCollection at 0x27b2cf8c400>

Let's see how many valid observations we have for each pixel.

In [18]:

S2_count = S2_ds.B2.count(dim='time')   # Count valid data values of one band along the time dimension
S2_count = S2_count.rename(new_name_or_name_dict='Valid observations')
S2_count.plot()

S2_count = S2_ds.B2.count(dim='time') # Count valid data values of one band along the time dimension S2_count = S2_count.rename(new_name_or_name_dict='Valid observations') S2_count.plot()

Out[18]:

<matplotlib.collections.QuadMesh at 0x27b2cfe88e0>

To plot an RGB image we create a new DataArray with the respective bands from our DataSet, then plot this (plotting functionality of xarray is mainly focused on DataArrays).

In [18]:

S2_RGB_da = S2_ds[['B4', 'B3', 'B2']].sel(time='2020-01-02').to_array()
S2_RGB_da.plot.imshow(robust=True)

S2_RGB_da = S2_ds[['B4', 'B3', 'B2']].sel(time='2020-01-02').to_array() S2_RGB_da.plot.imshow(robust=True)

Out[18]:

<matplotlib.image.AxesImage at 0x258b2a1b7f0>

Plot multiple observations within a defined time period, one band only:

In [19]:

S2_ds.B2.sel(time='2020-05').plot.imshow(col='time')

S2_ds.B2.sel(time='2020-05').plot.imshow(col='time')

Out[19]:

<xarray.plot.facetgrid.FacetGrid at 0x258b2d369b0>

The same for a 3-band combination (here SWIR-1, near-infrared and blue) and a longer time period:

In [20]:

# Create a new DataArray with the respective bands and time period, then plot this
S2_RGB_da = S2_ds[['B11', 'B8', 'B2']].sel(time=slice('2020-05-01', '2020-07-31')).to_array()
S2_RGB_da.plot.imshow(col='time', robust=True)

# Create a new DataArray with the respective bands and time period, then plot this S2_RGB_da = S2_ds[['B11', 'B8', 'B2']].sel(time=slice('2020-05-01', '2020-07-31')).to_array() S2_RGB_da.plot.imshow(col='time', robust=True)

Out[20]:

<xarray.plot.facetgrid.FacetGrid at 0x258b2dd3d90>

Band calculations¶

As an example for xarray's raster calculator functionality we calculate the Normalized Difference Snow Index (NDSI; Hall et al. 1995) for the entire Sentinel-2 time series. We can add it as a new variable to our Sentinel-2 Dataset, thereby keeping our data organised. Alternatively, we can define it as a new DataArray.

In [21]:

# Compute the NDSI from the SWIR-1 (B11) and green (B3) bands and add as a new variable to our Sentinel-2 Dataset
S2_ds['NDSI'] = (S2_ds.B3 - S2_ds.B11) / (S2_ds.B3 + S2_ds.B11)

# For convenience define NDSI as a new DataArray (Note: metadata, such as the CRS, are not propagated to the new DataArray by default)
NDSI = S2_ds.NDSI

# Show the resulting rasters of one month
# for other colormaps see https://matplotlib.org/stable/gallery/color/colormap_reference.html
NDSI.sel(time='2020-05').plot(col='time', cmap='ocean')
NDSI.sel(time='2020-05').plot(col='time', cmap='coolwarm_r', center=0.4, robust=True)    # center cmap around 0.4, a typical threshold for snow discrimination

# Compute the NDSI from the SWIR-1 (B11) and green (B3) bands and add as a new variable to our Sentinel-2 Dataset S2_ds['NDSI'] = (S2_ds.B3 - S2_ds.B11) / (S2_ds.B3 + S2_ds.B11) # For convenience define NDSI as a new DataArray (Note: metadata, such as the CRS, are not propagated to the new DataArray by default) NDSI = S2_ds.NDSI # Show the resulting rasters of one month # for other colormaps see https://matplotlib.org/stable/gallery/color/colormap_reference.html NDSI.sel(time='2020-05').plot(col='time', cmap='ocean') NDSI.sel(time='2020-05').plot(col='time', cmap='coolwarm_r', center=0.4, robust=True) # center cmap around 0.4, a typical threshold for snow discrimination

Out[21]:

<xarray.plot.facetgrid.FacetGrid at 0x258b79eb760>

Time series interpolation¶

For temporal interpolation and aggregation xarray contains some very useful functionality, and much of this works similar as with Pandas. In the rather short time series above we can see the snow disappearing and we also encounter different amounts of valid pixels (white NoData pixels should result from cloud-masking). As a next step, we will fill missing values with different methods.

Can you explain the method-specific differences in the resulting images? Why does the linear interpolation not fill all missing values?

In [22]:

# Plot the original (observed) NDSI values for one month
NDSI.sel(time='2020-05').plot(col='time', cmap='ocean')

# Interpolate missing values in these images and plot (default: linear interpolation)
NDSI.sel(time='2020-05').interpolate_na(dim='time').plot(col='time', cmap='ocean')#(time=new_dates)

# Forward fill missing values in these images and plot (requires the bottleneck package to be installed)
NDSI.sel(time='2020-05').ffill(dim='time').plot(col='time', cmap='ocean')

# Backward fill missing values in these images and plot
NDSI.sel(time='2020-05').bfill(dim='time').plot(col='time', cmap='ocean')

# Plot the original (observed) NDSI values for one month NDSI.sel(time='2020-05').plot(col='time', cmap='ocean') # Interpolate missing values in these images and plot (default: linear interpolation) NDSI.sel(time='2020-05').interpolate_na(dim='time').plot(col='time', cmap='ocean')#(time=new_dates) # Forward fill missing values in these images and plot (requires the bottleneck package to be installed) NDSI.sel(time='2020-05').ffill(dim='time').plot(col='time', cmap='ocean') # Backward fill missing values in these images and plot NDSI.sel(time='2020-05').bfill(dim='time').plot(col='time', cmap='ocean')

Out[22]:

<xarray.plot.facetgrid.FacetGrid at 0x258b9096800>

Via matplotlib we have many options for customizing plots (e.g., create subplots), and here we use xarray's possibility to make line plots directly from a DataArray:

In [23]:

import matplotlib.pyplot as plt                         # for extended plotting functionality

fig, axes = plt.subplots(1, 2, figsize=(15, 4))         # make two subplots next to each other

# Plot NDSI from the original time series
NDSI.isel(x=10, y=[10, 20, 30]).plot.line(ax=axes[0], x='time')

# Forward fill missing values of the entire NDSI time series, select a few locations with constant x and plot the NDSI against time
NDSI.ffill(dim='time').isel(x=10, y=[10, 20, 30]).plot.line(ax=axes[1], x='time')

plt.tight_layout()

import matplotlib.pyplot as plt # for extended plotting functionality fig, axes = plt.subplots(1, 2, figsize=(15, 4)) # make two subplots next to each other # Plot NDSI from the original time series NDSI.isel(x=10, y=[10, 20, 30]).plot.line(ax=axes[0], x='time') # Forward fill missing values of the entire NDSI time series, select a few locations with constant x and plot the NDSI against time NDSI.ffill(dim='time').isel(x=10, y=[10, 20, 30]).plot.line(ax=axes[1], x='time') plt.tight_layout()

... or heatmaps:

In [24]:

fig, axes = plt.subplots(1, 2, figsize=(10, 5))            # make two subplots next to each other

# Plot NDSI heatmap from the original time series (values at locations with constant x and varying y)
NDSI.isel(x=100, y=[10, 20, 30, 40, 50, 60, 70, 80, 90, 100]).plot(ax=axes[0], cmap='coolwarm_r', center=0.4, robust=True)

# Forward fill missing values of the entire NDSI time series, select a few locations with constant x and plot the NDSI heatmap
NDSI.ffill(dim='time').isel(x=100, y=[10, 20, 30, 40, 50, 60, 70, 80, 90, 100]).plot(ax=axes[1], cmap='coolwarm_r', center=0.4, robust=True)

plt.tight_layout()

fig, axes = plt.subplots(1, 2, figsize=(10, 5)) # make two subplots next to each other # Plot NDSI heatmap from the original time series (values at locations with constant x and varying y) NDSI.isel(x=100, y=[10, 20, 30, 40, 50, 60, 70, 80, 90, 100]).plot(ax=axes[0], cmap='coolwarm_r', center=0.4, robust=True) # Forward fill missing values of the entire NDSI time series, select a few locations with constant x and plot the NDSI heatmap NDSI.ffill(dim='time').isel(x=100, y=[10, 20, 30, 40, 50, 60, 70, 80, 90, 100]).plot(ax=axes[1], cmap='coolwarm_r', center=0.4, robust=True) plt.tight_layout()

It is also possible to interpolate other time steps (not covered by a Sentinel-2 overpass). Of course, such processing steps must be used carefully and considered during further analysis and interpretation. Note that the last "new" date cannot be interpolated because there is now data (in the selected subset) after this date.

In [25]:

# Interpolate new dates and create a new DataArray
new_dates = ['2020-05-05', '2020-05-10', '2020-05-15', '2020-05-20', '2020-05-25', '2020-05-30']
NDSI_new_dates_da = NDSI.sel(time='2020-05').interp(time=new_dates)

# Concatenate original data and data for new dates, and sort by 'time' (to maintain a chronological order)
NDSI_all_dates_da = xarray.concat(objs=[NDSI.sel(time='2020-05'), NDSI_new_dates_da], dim='time').sortby('time')

# Plot 
NDSI_all_dates_da.plot(col='time', cmap='ocean')

# Interpolate new dates and create a new DataArray new_dates = ['2020-05-05', '2020-05-10', '2020-05-15', '2020-05-20', '2020-05-25', '2020-05-30'] NDSI_new_dates_da = NDSI.sel(time='2020-05').interp(time=new_dates) # Concatenate original data and data for new dates, and sort by 'time' (to maintain a chronological order) NDSI_all_dates_da = xarray.concat(objs=[NDSI.sel(time='2020-05'), NDSI_new_dates_da], dim='time').sortby('time') # Plot NDSI_all_dates_da.plot(col='time', cmap='ocean')

Out[25]:

<xarray.plot.facetgrid.FacetGrid at 0x258bcbd2170>

Time series aggregation¶

Aggregating over temporal units is straightforward thanks to the time indexing. Here we group all observations of our NDSI time series by month and calculate the averages over these groups, then plot the resulting monthly composite rasters. This temporal aggregation gives a good impression of the typical snow cover distribution in our mountainous area-of-interest across the year.

In [26]:

NDSI.groupby('time.month').mean(dim='time').plot(col='month', cmap='ocean')

NDSI.groupby('time.month').mean(dim='time').plot(col='month', cmap='ocean')

Out[26]:

<xarray.plot.facetgrid.FacetGrid at 0x258bcc3be50>

If we are interested in the variation of snow cover in different months, we can calculate the variance (or the standard deviation) over the months:

In [27]:

# Group by month and compute the variance of NDSI, then plot
NDSI.groupby('time.month').var(dim='time').plot(col='month')

# Group by month and compute the variance of NDSI, then plot NDSI.groupby('time.month').var(dim='time').plot(col='month')

Out[27]:

<xarray.plot.facetgrid.FacetGrid at 0x258b931c6d0>

... or over seasons (named by the months' first letters):

In [28]:

# Group by seasons and compute the variance of NDSI, then plot
# I would prefer a chronological order of seasons but cannot get them sorted easily
NDSI.groupby('time.season').var(dim='time').plot(col='season')

# Group by seasons and compute the variance of NDSI, then plot # I would prefer a chronological order of seasons but cannot get them sorted easily NDSI.groupby('time.season').var(dim='time').plot(col='season')

Out[28]:

<xarray.plot.facetgrid.FacetGrid at 0x258c1207310>

Now let's have a look at monthly NDSI aggregates of different years. We resample to one month temporal resolution by computing the mean, then make a facet plot with 12 columns, so that we get each year in a separate row.

In [29]:

NDSI.resample(time='1M').mean().plot(col='time', col_wrap=12, cmap='coolwarm_r', center=0.4, robust=True)

NDSI.resample(time='1M').mean().plot(col='time', col_wrap=12, cmap='coolwarm_r', center=0.4, robust=True)

Out[29]:

<xarray.plot.facetgrid.FacetGrid at 0x258c9b535e0>

There seem to be no images from May 2018 (resulting in a NoData aggregate), and we may be interested in the number of observations for all other months as well. Using count() instead of mean(), we can inspect the number of valid observations per month:

In [30]:

NDSI.resample(time='1M').count().plot(col='time', col_wrap=12)

NDSI.resample(time='1M').count().plot(col='time', col_wrap=12)

Out[30]:

<xarray.plot.facetgrid.FacetGrid at 0x258d32b9b70>

Time series smoothing¶

A couple of other operations related to time series processing are straightforward with xarrays, such as smoothing and filtering with rolling windows.

To test such methods, we first resample to a higher temporal resolution ("upsampling" in the xarrays documentation) by linear interpolation. This will not work if there are duplicates in the time index (e.g. due to S2 products sensed on the same day but available in two versions with different product generation time stamps). Then you will get an "InvalidIndexError" (Reindexing only valid with uniquely valued Index objects). So we treat the duplicates first, e.g. by searching, printing and removing them.

In [31]:

# We use pandas' .duplicated() to check for duplicates
import pandas as pd

# Print the dates with duplicates
for i,j in zip(pd.Index(NDSI['time']).duplicated(), NDSI['time'].values):
    if i == True:
        print(j)

# We use pandas' .duplicated() to check for duplicates import pandas as pd # Print the dates with duplicates for i,j in zip(pd.Index(NDSI['time']).duplicated(), NDSI['time'].values): if i == True: print(j)

2019-03-28T00:00:00.000000000
2020-10-18T00:00:00.000000000
2020-10-28T00:00:00.000000000
2020-11-17T00:00:00.000000000
2020-11-27T00:00:00.000000000
2020-12-17T00:00:00.000000000
2021-01-16T00:00:00.000000000
2021-02-15T00:00:00.000000000
2021-09-03T00:00:00.000000000

In [32]:

# Remove duplicates and keep only the first occurrences (default: keep='first')) in a new DataArray with unique time index
NDSI_unique = NDSI.drop_duplicates('time')

# Compare the shape of the two DataArrays
print(NDSI.shape)
print(NDSI_unique.shape)

# Remove duplicates and keep only the first occurrences (default: keep='first')) in a new DataArray with unique time index NDSI_unique = NDSI.drop_duplicates('time') # Compare the shape of the two DataArrays print(NDSI.shape) print(NDSI_unique.shape)

(335, 201, 251)
(326, 201, 251)

Now we resample to one day temporal resolution and fill daily values by linear interpolation between observations.

In [33]:

# Resample to 1 day temporal resolution
NDSI_interpolated = NDSI_unique.resample(time='1d').interpolate()

# Interpolate the missing values along the time dimension
NDSI_interpolated = NDSI_interpolated.interpolate_na(dim='time')

# Plot the resulting time series at three points
NDSI_interpolated.isel(x=100, y=[50, 100, 150]).plot.line(x='time')

# Or in one line of code
#NDSI_unique.resample(time='1d').interpolate().interpolate_na(dim='time').isel(x=100, y=[50, 100, 150]).plot.line(x='time')

# Resample to 1 day temporal resolution NDSI_interpolated = NDSI_unique.resample(time='1d').interpolate() # Interpolate the missing values along the time dimension NDSI_interpolated = NDSI_interpolated.interpolate_na(dim='time') # Plot the resulting time series at three points NDSI_interpolated.isel(x=100, y=[50, 100, 150]).plot.line(x='time') # Or in one line of code #NDSI_unique.resample(time='1d').interpolate().interpolate_na(dim='time').isel(x=100, y=[50, 100, 150]).plot.line(x='time')

Out[33]:

[<matplotlib.lines.Line2D at 0x258deee70d0>,
 <matplotlib.lines.Line2D at 0x258def2f7f0>,
 <matplotlib.lines.Line2D at 0x258def2f760>]

Compute the 30-day rolling average of the NDSI time series and plot for one point (aggregation results are assigned the coordinate at the end of each window by default, but can be centered by passing center=True when constructing the rolling object):

In [34]:

# Set up the plot
fig, axes = plt.subplots(1, 1, figsize=(12, 4))

# The daily NDSI values
NDSI_interpolated.isel(x=100, y=50).plot.line(label='Interpolated NDSI')

# Take the mean value over the NDSI values in a 30-day rolling window and plot the time series for one specific point
NDSI_interpolated.rolling(time=30, center=True).mean().isel(x=100, y=50).plot.line(label='Interpolated NDSI, 30-days rolling average')

# Show the legend and show the plot
fig.legend()
plt.tight_layout()

# Set up the plot fig, axes = plt.subplots(1, 1, figsize=(12, 4)) # The daily NDSI values NDSI_interpolated.isel(x=100, y=50).plot.line(label='Interpolated NDSI') # Take the mean value over the NDSI values in a 30-day rolling window and plot the time series for one specific point NDSI_interpolated.rolling(time=30, center=True).mean().isel(x=100, y=50).plot.line(label='Interpolated NDSI, 30-days rolling average') # Show the legend and show the plot fig.legend() plt.tight_layout()