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: 2025-03-05
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 module1/module1_data/T1/S2 in the ETRAINEE course datasets.

In [1]:

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 [ ]:

# Define the path to the s2 directory with the Sentinel-2 image subsets on your local system
in_dir = pathlib.Path('F:/data/etrainee_m1/module1/module1_data/T1/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)

# Drop potential duplicates along the time dimension (if there are multiple images for a time stamp, keep only the first one)
S2_ds = S2_ds.drop_duplicates(dim='time', keep='first')

# Define the path to the s2 directory with the Sentinel-2 image subsets on your local system in_dir = pathlib.Path('F:/data/etrainee_m1/module1/module1_data/T1/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) # Drop potential duplicates along the time dimension (if there are multiple images for a time stamp, keep only the first one) S2_ds = S2_ds.drop_duplicates(dim='time', keep='first')

Many time series operations (such as facetgrid plots and resampling) 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 (e.g.) an "InvalidIndexError" (Reindexing only valid with uniquely valued Index objects). Therefore, we removed duplicates above.

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 [3]:

# 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[3]:

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

xarray.Dataset

• Dimensions:
  • time: 326
  • 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
    ()
    int64
    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.],
            ...,
    ...
            [  519.,   497.,   483., ...,   793.,   922.,   792.],
            [  513.,   468.,   433., ...,   786.,   936.,   693.],
            [  506.,   488.,   426., ...,   835.,   716.,   610.]],
    
           [[  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.],
            ...,
    ...
            [  813.,   789.,   799., ...,  1090.,  1225.,  1075.],
            [  825.,   756.,   663., ...,  1097.,  1228.,   960.],
            [  831.,   787.,   633., ...,  1152.,  1068.,   928.]],
    
           [[  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.],
            ...,
    ...
            [  787.,   772.,   770., ...,  1206.,  1364.,  1176.],
            [  791.,   732.,   695., ...,  1190.,  1373.,  1004.],
            [  758.,   762.,   646., ...,  1255.,  1099.,   992.]],
    
           [[ 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.],
            ...,
    ...
            [ 3280.,  3100.,  3430., ...,  1959.,  1776.,  1461.],
            [ 3425.,  3175.,  2763., ...,  2140.,  1954.,  2193.],
            [ 3502.,  3333.,  2547., ...,  2069.,  2578.,  2763.]],
    
           [[ 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.],
            ...,
    ...
            ...,
            [2393., 2310., 2296., ..., 1985., 1857., 1723.],
            [2428., 2220., 2095., ..., 1926., 1935., 2013.],
            [2460., 2326., 2036., ..., 2024., 2252., 2377.]],
    
           [[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.],
            ...,
    ...
            ...,
            [1349., 1283., 1255., ..., 1513., 1478., 1324.],
            [1357., 1239., 1174., ..., 1441., 1462., 1459.],
            [1395., 1274., 1118., ..., 1461., 1519., 1597.]],
    
           [[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(Index([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(Index([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-14', '2021-08-16', '2021-08-19', '2021-08-21',
                   '2021-08-24', '2021-08-26', '2021-08-31', '2021-09-03',
                   '2021-09-05', '2021-09-08'],
                  dtype='datetime64[ns]', name='time', length=326, freq=None))

• Attributes: (9)

  TIFFTAG_XRESOLUTION :

  1

  TIFFTAG_YRESOLUTION :

  1

  TIFFTAG_RESOLUTIONUNIT :

  1 (unitless)

  AREA_OR_POINT :

  Area

  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 [4]:

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[4]:

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

xarray.Dataset

• Dimensions:
  • time: 326
  • 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.],
            ...,
    ...
            [  519.,   497.,   483., ...,   793.,   922.,   792.],
            [  513.,   468.,   433., ...,   786.,   936.,   693.],
            [  506.,   488.,   426., ...,   835.,   716.,   610.]],
    
           [[  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.],
            ...,
    ...
            [  813.,   789.,   799., ...,  1090.,  1225.,  1075.],
            [  825.,   756.,   663., ...,  1097.,  1228.,   960.],
            [  831.,   787.,   633., ...,  1152.,  1068.,   928.]],
    
           [[  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.],
            ...,
    ...
            [  787.,   772.,   770., ...,  1206.,  1364.,  1176.],
            [  791.,   732.,   695., ...,  1190.,  1373.,  1004.],
            [  758.,   762.,   646., ...,  1255.,  1099.,   992.]],
    
           [[ 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.],
            ...,
    ...
            [ 3280.,  3100.,  3430., ...,  1959.,  1776.,  1461.],
            [ 3425.,  3175.,  2763., ...,  2140.,  1954.,  2193.],
            [ 3502.,  3333.,  2547., ...,  2069.,  2578.,  2763.]],
    
           [[ 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.],
            ...,
    ...
            ...,
            [2393., 2310., 2296., ..., 1985., 1857., 1723.],
            [2428., 2220., 2095., ..., 1926., 1935., 2013.],
            [2460., 2326., 2036., ..., 2024., 2252., 2377.]],
    
           [[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.],
            ...,
    ...
            ...,
            [1349., 1283., 1255., ..., 1513., 1478., 1324.],
            [1357., 1239., 1174., ..., 1441., 1462., 1459.],
            [1395., 1274., 1118., ..., 1461., 1519., 1597.]],
    
           [[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(Index([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(Index([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-14', '2021-08-16', '2021-08-19', '2021-08-21',
                   '2021-08-24', '2021-08-26', '2021-08-31', '2021-09-03',
                   '2021-09-05', '2021-09-08'],
                  dtype='datetime64[ns]', name='time', length=326, freq=None))

• Attributes: (9)

  TIFFTAG_XRESOLUTION :

  1

  TIFFTAG_YRESOLUTION :

  1

  TIFFTAG_RESOLUTIONUNIT :

  1 (unitless)

  AREA_OR_POINT :

  Area

  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 [5]:

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

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

Out[5]:

<matplotlib.collections.QuadMesh at 0x1e380277950>

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

In [6]:

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[6]:

<matplotlib.collections.PathCollection at 0x1e380670e00>

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

In [7]:

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[7]:

<matplotlib.collections.QuadMesh at 0x1e3806890a0>

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 [8]:

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

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

Out[8]:

<matplotlib.image.AxesImage at 0x1e3807f57f0>

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

In [9]:

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

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

Out[9]:

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

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

In [10]:

# 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[10]:

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

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 [11]:

# 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[11]:

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

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 [12]:

# 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[12]:

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

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 [13]:

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 [14]:

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 no data (in the selected subset) after this date.

In [15]:

# 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[15]:

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

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 [16]:

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[16]:

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

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 [17]:

# 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[17]:

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

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

In [18]:

# 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[18]:

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

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. For an overview of time series groupers/frequencies, see here.

In [22]:

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

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

Out[22]:

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

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 [24]:

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

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

Out[24]:

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

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 (one day) ("upsampling" in the xarrays documentation) and fill daily values by linear interpolation between observations.

In [28]:

# Resample to 1 day temporal resolution
NDSI_interpolated = NDSI.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.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.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.resample(time='1d').interpolate().interpolate_na(dim='time').isel(x=100, y=[50, 100, 150]).plot.line(x='time')

Out[28]:

[<matplotlib.lines.Line2D at 0x1e411dc09b0>,
 <matplotlib.lines.Line2D at 0x1e40f27fce0>,
 <matplotlib.lines.Line2D at 0x1e3eea95df0>]

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 [29]:

# 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()