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. Many of the visualizations in this tutorial could be improved with a few lines of code (e.g. to add axes labels and plot titles) or by using specialized plotting packages, but here we try to keep the code simple.

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.

Getting started with xarray's data structures and basic tools¶

Loading and handling single images¶

Before we deal with a time series of satellite images, we start with importing and exploring an individual image. So we load a single GeoTiff file with six bands, which is part of the image collection subset downloaded from the GEE. A good way to do this is using rioxarray, a geospatial xarray extension powered by rasterio.

In [1]:

import rioxarray
import xarray
import pathlib

import rioxarray import xarray import pathlib

In [2]:

# Open a GeoTiff from a local directory using the rio accessor
in_dir = pathlib.Path('F:/data/etrainee/gee/s2')        # Define the path to the s2 directory with the Sentinel-2 image subsets on your local system
in_file = in_dir / '20180726T102019_20180726T102150_T32TPS.tif' # Specify the path and the filename of an image
xds = rioxarray.open_rasterio(in_file, masked=True)     # Read the mask (NoData in the GeoTiff) and set values to NaN (for proper interpretation by xarray)

# Open a GeoTiff from a local directory using the rio accessor in_dir = pathlib.Path('F:/data/etrainee/gee/s2') # Define the path to the s2 directory with the Sentinel-2 image subsets on your local system in_file = in_dir / '20180726T102019_20180726T102150_T32TPS.tif' # Specify the path and the filename of an image xds = rioxarray.open_rasterio(in_file, masked=True) # Read the mask (NoData in the GeoTiff) and set values to NaN (for proper interpretation by xarray)

Print some information about the xarray.DataArray called 'xds':

In [3]:

xds

xds

Out[3]:

<xarray.DataArray (band: 6, y: 201, x: 251)>
[302706 values with dtype=float32]
Coordinates:
  * band         (band) int32 1 2 3 4 5 6
  * 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
Attributes:
    AREA_OR_POINT:           Area
    TIFFTAG_RESOLUTIONUNIT:  1 (unitless)
    TIFFTAG_XRESOLUTION:     1
    TIFFTAG_YRESOLUTION:     1
    scale_factor:            1.0
    add_offset:              0.0

xarray.DataArray

• band: 6
• y: 201
• x: 251

• ...

  [302706 values with dtype=float32]

• Coordinates: (4)
  • band
    (band)
    int32
    1 2 3 4 5 6

    array([1, 2, 3, 4, 5, 6])

  • 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)

• Indexes: (3)
  • band
    PandasIndex

    PandasIndex(Int64Index([1, 2, 3, 4, 5, 6], dtype='int64', name='band'))

  • 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))

• Attributes: (6)

  AREA_OR_POINT :

  Area

  TIFFTAG_RESOLUTIONUNIT :

  1 (unitless)

  TIFFTAG_XRESOLUTION :

  1

  TIFFTAG_YRESOLUTION :

  1

  scale_factor :

  1.0

  add_offset :

  0.0

Or more specifically:

In [4]:

print(xds.shape)
print(xds.dims)
print(xds.coords)

print(xds.shape) print(xds.dims) print(xds.coords)

(6, 201, 251)
('band', 'y', 'x')
Coordinates:
  * band         (band) int32 1 2 3 4 5 6
  * 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

... and some descriptive statistics for one band (using xarray's built-in methods, but numpy could be used as well):

In [5]:

# Extract a band to a new DataArray (for conveniance)
band_1 =xds.sel(band=1)
print('min: ', band_1.min().values)     # Print only the value(s) of the array; its only one value for one band and one time stamp 
print('max: ', band_1.max().values)
print('mean: ', band_1.mean().values)
print('std: ', band_1.std().values)

# Extract a band to a new DataArray (for conveniance) band_1 =xds.sel(band=1) print('min: ', band_1.min().values) # Print only the value(s) of the array; its only one value for one band and one time stamp print('max: ', band_1.max().values) print('mean: ', band_1.mean().values) print('std: ', band_1.std().values)

min:  138.0
max:  6651.0
mean:  666.3916
std:  358.31906

We can use the extended capabilities of rioxarray to query more metadata, such as coordinate reference system (CRS), bounding coordinates, dimensions, spatial resolution, or defined NoData value:

In [6]:

print(xds.rio.crs)
print(xds.rio.bounds())
print(xds.rio.width)
print(xds.rio.height)
print(xds.rio.resolution())
print(xds.rio.nodata)

print(xds.rio.crs) print(xds.rio.bounds()) print(xds.rio.width) print(xds.rio.height) print(xds.rio.resolution()) print(xds.rio.nodata)

EPSG:32632
(651580.0, 5189740.0, 656600.0, 5193760.0)
251
201
(20.0, -20.0)
nan

Creating a new DataArray with only one selected band is straigthforward and plotting this is easy as well (matplotlib must be installed in the active conda environment). Note the white NoData areas due to cloud masking:

In [7]:

band_1 = xds.sel(band=1)            # Select a band by its label (or by its index with '.isel()'), as already shown above
band_1.plot.imshow(robust=True)     # The 'robust' option stretches the colormap range to 2nd and 98th percentiles.

band_1 = xds.sel(band=1) # Select a band by its label (or by its index with '.isel()'), as already shown above band_1.plot.imshow(robust=True) # The 'robust' option stretches the colormap range to 2nd and 98th percentiles.

Out[7]:

<matplotlib.image.AxesImage at 0x2470d15f580>

Now let's create a histogram for this band. We use xarray.plot.hist(), for more advanced options there are the xhistogram and the boost-histogram packages. You might try different estimators to determine the optimum number of bins (see here; default would be 10 bins; 'auto' takes the maximum of the Sturges and Freedman Diaconis estimators).

In [8]:

xds.sel(band=1).plot.hist(bins='auto'); # The semicolon tells Jupyter Notebook not to display the output of that line of code (i.e. the bin values; but it will show the histogram).

xds.sel(band=1).plot.hist(bins='auto'); # The semicolon tells Jupyter Notebook not to display the output of that line of code (i.e. the bin values; but it will show the histogram).

Similarly, we can plot 3-band combinations as true-color or false-color (color-infrared, CIR) image by selecting the bands to plot as RGB:

In [9]:

xds.sel(band=[3, 2, 1]).plot.imshow(robust=True)

xds.sel(band=[3, 2, 1]).plot.imshow(robust=True)

Out[9]:

<matplotlib.image.AxesImage at 0x2470de47370>

In [10]:

xds.sel(band=[4, 3, 2]).plot.imshow(robust=True)

xds.sel(band=[4, 3, 2]).plot.imshow(robust=True)

Out[10]:

<matplotlib.image.AxesImage at 0x2470e003940>

Instead of using robust=True we could also constrain the color stretch to the respective quantiles using the vmin and vmax keyword arguments. Try also other quantiles or explicit values and see how this affects the image display.

In [11]:

xds.sel(band=[3, 2, 1]).plot.imshow(
    vmin=xds.sel(band=[3, 2, 1]).quantile(0.02),
    vmax=xds.sel(band=[3, 2, 1]).quantile(0.98))

xds.sel(band=[3, 2, 1]).plot.imshow( vmin=xds.sel(band=[3, 2, 1]).quantile(0.02), vmax=xds.sel(band=[3, 2, 1]).quantile(0.98))

Out[11]:

<matplotlib.image.AxesImage at 0x2470e08e6e0>

That's nice, but how can we handle a remote sensing time series?

Loading many files of an image time series¶

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

In [12]:

import glob
import os
from datetime import datetime

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

import glob import os from datetime import datetime # 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 [13]:

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

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

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

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

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

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

Out[15]:

<matplotlib.collections.QuadMesh at 0x247105818a0>

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

In [16]:

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

<matplotlib.collections.PathCollection at 0x2470dce5720>

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

In [17]:

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

<matplotlib.collections.QuadMesh at 0x2471068b910>

To plot an RGB image we create a new DataArray with the respective bands from our DataSet, then plot this.

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 0x2470da11ed0>

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 0x24710659c90>

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 0x2471081b4c0>

Band calculations¶

As an example for xarray's raster calculator functionality we calculate the Normalized Difference Snow Index (NDSI; Hall et al. 199500137-P)) 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 0x2471297a2c0>

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

In [45]:

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

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

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

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

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

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

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

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

Out[32]:

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

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

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

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

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

[<matplotlib.lines.Line2D at 0x247392ae980>,
 <matplotlib.lines.Line2D at 0x24739c8ada0>,
 <matplotlib.lines.Line2D at 0x24739c8ad10>]

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

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