Welcome to EASI

This notebook introduces new users to working with EASI notebooks and the Open Data Cube.

It will demonstrate the following basic functionality:

  • Notebook setup
  • Connect to the OpenDataCube
    • List products
    • List measurements
    • Choose a region of interest
    • Load data
    • Plot data
    • Masking and Scaling
    • Perform a calculation on the data
    • Save the results to a file
  • Summary
  • Further reading

Notebook setup¶

A notebook consists of cells that contain either text descriptions or python code for performing operations on data.

  1. Start by clicking on the cell below to select it.
  2. Execute a selected cell, or each cell in sequence, by clicking the "play" button (in the toolbar above) or pressing Shift+Enter.
  3. Each cell will show an asterisk icon [*]: when it is running. Once this changes to a number, the cell has finished.
  4. This cell below imports packages to use and sets some formatting options.
In [1]:
# Formatting for basic plots
%matplotlib inline
%config InlineBackend.rc = {}
import matplotlib.pyplot as plt
# plt.rcParams['figure.figsize'] = [12, 8]

# Formatting pandas table output
import pandas
pandas.set_option("display.max_rows", None)

# Datacube
import datacube
from datacube.utils import masking  # https://github.com/opendatacube/datacube-core/blob/develop/datacube/utils/masking.py
from odc.algo import enum_to_bool   # https://github.com/opendatacube/odc-tools/blob/develop/libs/algo/odc/algo/_masking.py
from datacube.utils.rio import configure_s3_access

# Notebook helper tools (in dea_tools or in this repo)
import sys
try:
    from dea_tools.plotting import display_map, rgb
except ImportError:
    # Local copy of selected dea_tools
    if 'tools/' not in sys.path:
        sys.path.append('tools/')
    from datacube_utils import display_map
    rgb = None  # Not copied or adapted yet

Connect to the OpenDataCube¶

The Datacube() API provides search, load and information functions for data products indexed in the EASI ODC database. For further information see https://datacube-core.readthedocs.io/en/latest/.

In [2]:
dc = datacube.Datacube()

# Optional: Access AWS "requester-pays" buckets
# This is necessary for Landsat ("landsatN_c2l2_*") and Sentinel-2 ("s2_l2a") products
from datacube.utils.aws import configure_s3_access
configure_s3_access(aws_unsigned=False, requester_pays=True);

List products¶

Get all available products in the ODC database and list them along with selected properties.

  1. View available products and data coverage at the EASI Explorer: https://explorer.asia.easi-eo.solutions
In [3]:
products = dc.list_products()  # Pandas DataFrame
products
Out[3]:
name description license default_crs default_resolution
name
copernicus_dem_30 copernicus_dem_30 Copernicus 30m Digital Elevation Model (GLO-30) None None None
landsat5_c2l2_sr landsat5_c2l2_sr Landsat 5 Collection 2 Level-2 Surface Reflect... None None None
landsat5_c2l2_st landsat5_c2l2_st Landsat 5 Collection 2 Level-2 UTM Surface Tem... CC-BY-4.0 None None
landsat7_c2l2_sr landsat7_c2l2_sr Landsat 7 USGS Collection 2 Surface Reflectanc... None None None
landsat7_c2l2_st landsat7_c2l2_st Landsat 7 Collection 2 Level-2 UTM Surface Tem... CC-BY-4.0 None None
landsat8_c2l1 landsat8_c2l1 Landsat 8 Collection 2 Level-1 (top of atmosph... CC-BY-4.0 None None
landsat8_c2l2_sr landsat8_c2l2_sr Landsat 8 Collection 2 Surface Reflectance, pr... None None None
landsat8_c2l2_st landsat8_c2l2_st Landsat 8 Collection 2 Level-2 UTM Surface Tem... CC-BY-4.0 None None
landsat9_c2l1 landsat9_c2l1 Landsat 9 Collection 2 Level-1 (top of atmosph... CC-BY-4.0 None None
landsat9_c2l2_sr landsat9_c2l2_sr Landsat 9 Collection 2 Surface Reflectance, pr... CC-BY-4.0 None None
landsat9_c2l2_st landsat9_c2l2_st Landsat 9 Collection 2 Level-2 UTM Surface Tem... CC-BY-4.0 None None
lpdaac_nasadem lpdaac_nasadem NASADEM Merged DEM Global 1 arc second V001 None None None
nasa_aqua_l2_oc nasa_aqua_l2_oc NASA MODIS-Aqua L2 Ocean Color, regridded to W... None EPSG:4326 (0.01, 0.01)
nasa_aqua_l2_sst nasa_aqua_l2_sst NASA MODIS-Aqua L2 Sea Surface Temperature, re... None EPSG:4326 (0.01, 0.01)
s2_l2a s2_l2a Sentinel-2a and Sentinel-2b imagery, processed... None None None

List measurements¶

The data arrays available for each product are called "measurements".

  1. List the "measurements" of the Landsat-8 surface reflectance product ("landsat8_c2l2_sr")
  2. The columns are the "attributes" available for each "measurement"
In [4]:
measurements = dc.list_measurements()  # Pandas DataFrame, all products
measurements.loc['landsat8_c2l2_sr']
Out[4]:
name dtype units nodata aliases flags_definition
measurement
coastal coastal uint16 reflectance 0 [SR_B1, band_1, B1, coastal_aerosol] NaN
blue blue uint16 reflectance 0 [SR_B2, band_2, B2] NaN
green green uint16 reflectance 0 [SR_B3, band_3, B3] NaN
red red uint16 reflectance 0 [SR_B4, band_4, B4] NaN
nir08 nir08 uint16 reflectance 0 [SR_B5, band_5, B5, nir] NaN
swir16 swir16 uint16 reflectance 0 [SR_B6, band_6, B6, swir1] NaN
swir22 swir22 uint16 reflectance 0 [SR_B7, band_7, B7, swir2] NaN
qa_pixel qa_pixel uint16 bit_index 1 [pixel_quality, level2_qa, QA_PIXEL, pixel_qa] {'snow': {'bits': 5, 'values': {'0': 'not_high...
qa_aerosol qa_aerosol uint8 bit_index 1 [sr_aerosol, SR_QA_AEROSOL, aerosol_qa] {'water': {'bits': 2, 'values': {'0': 'not_wat...
qa_radsat qa_radsat uint16 bit_index 1 [saturation_qa, QA_RADSAT, radsat_qa] {'qa_radsat': {'bits': [0, 1, 2, 3, 4, 5, 6, 7...

Choose a region of interest¶

See the available latitude/longitude and time ranges in the ODC Explorer.

  1. Feel free to change the latitude/longitude or time ranges of the query below.
In [5]:
# hub.asia.easi-eo.solutions - Lake Tempe, Indonesia
latitude = (-4.2, -3.9)
longitude = (119.8, 120.1)
display_map(longitude, latitude)
Out[5]:
Make this Notebook Trusted to load map: File -> Trust Notebook

Load data¶

Load Landsat-8 surface reflectance data ("landsat8_c2l2_sr") for a given latitude, longitude and time range.

This may take a few minutes while the data are loaded into JupyterLab (so choose a small area and time range!).

The datacube.load() function returns an xarray.Dataset object.

  • The display() view is a convenient way to check the data request size, shape and attributes
  • Once you have an xarray object this can be used with many Python packages.

The output_crs and resolution parameters are dependent on the product chosen.

Further information on xarray: http://xarray.pydata.org/en/stable/user-guide/data-structures.html

In [6]:
# A standard datacube.load() call
# Any SQLAlchemy warnings can be ignored.

data = dc.load(
    product = 'landsat8_c2l2_sr', 
    latitude = latitude,
    longitude = longitude,
    time=('2020-02-01', '2020-04-01'),
    
    output_crs="EPSG:32655",  # Target CRS
    resolution=(30, -30),     # Target resolution
    group_by='solar_day',     # Group by time method
)

display(data)

<xarray.Dataset>
Dimensions:      (time: 4, y: 1286, x: 1291)
Coordinates:
  * time         (time) datetime64[ns] 2020-02-10T02:10:17.788921 ... 2020-03...
  * y            (y) float64 -5.222e+05 -5.221e+05 ... -4.836e+05 -4.836e+05
  * x            (x) float64 -2.601e+06 -2.601e+06 ... -2.64e+06 -2.64e+06
    spatial_ref  int32 32655
Data variables:
    coastal      (time, y, x) uint16 8203 7860 8172 8332 ... 11387 11498 11471
    blue         (time, y, x) uint16 8579 8101 8434 8662 ... 11459 11543 11542
    green        (time, y, x) uint16 9795 9583 9859 10183 ... 12887 12939 12823
    red          (time, y, x) uint16 9550 9115 9580 9974 ... 12978 13125 13015
    nir08        (time, y, x) uint16 21181 19938 19005 ... 20020 20151 20301
    swir16       (time, y, x) uint16 15521 14502 15444 ... 11572 11694 11635
    swir22       (time, y, x) uint16 11529 10801 11726 ... 10941 11000 10982
    qa_pixel     (time, y, x) uint16 21824 21824 21824 ... 22280 22280 22280
    qa_aerosol   (time, y, x) uint8 160 160 160 160 160 ... 224 224 224 224 224
    qa_radsat    (time, y, x) uint16 0 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 0 0
Attributes:
    crs:           EPSG:32655
    grid_mapping:  spatial_ref
xarray.Dataset
    • time: 4
    • y: 1286
    • x: 1291
    • time
      (time)
      datetime64[ns]
      2020-02-10T02:10:17.788921 ... 2...
      units :
      seconds since 1970-01-01 00:00:00
      array(['2020-02-10T02:10:17.788921000', '2020-02-26T02:10:13.879802000',
       '2020-03-13T02:10:07.386121000', '2020-03-29T02:09:58.450496000'],
      dtype='datetime64[ns]')
    • y
      (y)
      float64
      -5.222e+05 ... -4.836e+05
      units :
      metre
      resolution :
      30.0
      crs :
      EPSG:32655
      array([-522165., -522135., -522105., ..., -483675., -483645., -483615.])
    • x
      (x)
      float64
      -2.601e+06 -2.601e+06 ... -2.64e+06
      units :
      metre
      resolution :
      -30.0
      crs :
      EPSG:32655
      array([-2600925., -2600955., -2600985., ..., -2639565., -2639595., -2639625.])
    • spatial_ref
      ()
      int32
      32655
      spatial_ref :
      PROJCS["WGS 84 / UTM zone 55N",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",147],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","32655"]]
      grid_mapping_name :
      transverse_mercator
      array(32655, dtype=int32)
    • coastal
      (time, y, x)
      uint16
      8203 7860 8172 ... 11498 11471
      units :
      reflectance
      nodata :
      0
      crs :
      EPSG:32655
      grid_mapping :
      spatial_ref
      array([[[ 8203,  7860,  8172, ..., 23385, 24062, 23896],
        [ 7903,  8052,  8277, ..., 23956, 24340, 24695],
        [ 7872,  8026,  8035, ..., 23956, 24374, 24641],
        ...,
        [ 8184,  8223,  8125, ...,  8335,  8299,  8317],
        [ 8250,  8272,  8166, ...,  8300,  8324,  8443],
        [ 8162,  8252,  8317, ...,  8344,  8395,  8455]],

       [[ 5731,  4677,  6334, ..., 25484, 25834, 25918],
        [ 2027,  5425,  5112, ..., 20029, 21149, 17172],
        [ 5962,  8546,  7365, ..., 20029, 11914,  9865],
        ...,
        [23256, 23282, 23406, ...,  8772,  8966,  8632],
        [23236, 23304, 23418, ...,  8668,  8897,  8746],
        [23330, 23364, 23456, ...,  8357,  8718,  8759]],

       [[ 8341,  8232,  8517, ...,  7543,  7568,  8258],
        [ 8108,  8424,  8575, ...,  7691,  7632, 10381],
        [ 8091,  8313,  8365, ...,  7691,  7854, 13159],
        ...,
        [ 8181,  8174,  8111, ...,  8140,  8109,  8031],
        [ 8147,  8222,  8164, ...,  8268,  8222,  8119],
        [ 8174,  8266,  8229, ...,  8202,  8146,  8072]],

       [[11745, 12045, 11835, ...,  7450,  7528,  7568],
        [11596, 11743, 11592, ...,  9185,  9348,  8515],
        [11870, 11814, 11509, ...,  9185, 10443,  9642],
        ...,
        [ 8035,  8084,  8141, ..., 11222, 11304, 11263],
        [ 7857,  7930,  7945, ..., 11372, 11459, 11424],
        [ 7728,  7857,  7797, ..., 11387, 11498, 11471]]], dtype=uint16)
    • blue
      (time, y, x)
      uint16
      8579 8101 8434 ... 11543 11542
      units :
      reflectance
      nodata :
      0
      crs :
      EPSG:32655
      grid_mapping :
      spatial_ref
      array([[[ 8579,  8101,  8434, ..., 23614, 24180, 24089],
        [ 8159,  8257,  8577, ..., 23908, 24497, 24774],
        [ 8134,  8280,  8282, ..., 23908, 24339, 24629],
        ...,
        [ 8275,  8267,  8214, ...,  8831,  8836,  8937],
        [ 8303,  8338,  8243, ...,  8869,  8905,  9085],
        [ 8265,  8328,  8375, ...,  8870,  8965,  9065]],

       [[ 6440,  5902,  7068, ..., 25311, 26023, 26147],
        [ 3587,  5831,  6975, ..., 19654, 20793, 16669],
        [ 5895,  8919,  7469, ..., 19654, 11626,  9810],
        ...,
        [22978, 22921, 23041, ...,  8886,  9256,  9018],
        [22969, 23003, 23061, ...,  8834,  9103,  9095],
        [23026, 22977, 23091, ...,  8712,  8951,  8981]],

       [[ 8479,  8275,  8635, ...,  7710,  7806,  8633],
        [ 8167,  8513,  8719, ...,  7921,  7812, 10706],
        [ 8216,  8433,  8473, ...,  7921,  8068, 13523],
        ...,
        [ 8204,  8196,  8125, ...,  8229,  8257,  8105],
        [ 8199,  8270,  8189, ...,  8347,  8300,  8210],
        [ 8192,  8332,  8296, ...,  8309,  8282,  8180]],

       [[11980, 12119, 11942, ...,  9303,  8480,  8264],
        [11667, 11852, 11763, ...,  9731,  9995,  9240],
        [11780, 11720, 11618, ...,  9731, 11009, 10382],
        ...,
        [ 7993,  8028,  8045, ..., 11392, 11485, 11417],
        [ 7831,  7915,  7856, ..., 11405, 11522, 11522],
        [ 7674,  7846,  7764, ..., 11459, 11543, 11542]]], dtype=uint16)
    • green
      (time, y, x)
      uint16
      9795 9583 9859 ... 12939 12823
      units :
      reflectance
      nodata :
      0
      crs :
      EPSG:32655
      grid_mapping :
      spatial_ref
      array([[[ 9795,  9583,  9859, ..., 23635, 23839, 23708],
        [ 9471,  9746, 10035, ..., 24086, 24252, 24225],
        [ 9483,  9669,  9640, ..., 24086, 24609, 24647],
        ...,
        [ 9210,  9195,  9193, ..., 10587, 10605, 10656],
        [ 9388,  9364,  9103, ..., 10580, 10574, 10852],
        [ 9452,  9355,  9303, ..., 10626, 10683, 10814]],

       [[ 9698,  9612, 10338, ..., 27049, 26963, 26998],
        [ 9153, 10689, 11470, ..., 22575, 23097, 20206],
        [10594, 11051, 11923, ..., 22575, 16425, 13676],
        ...,
        [22653, 22721, 22784, ..., 10095, 10522, 10438],
        [22614, 22658, 22744, ..., 10617, 10535, 10354],
        [22672, 22680, 22780, ..., 10216, 10463, 10227]],

       [[ 9422,  9237,  9608, ...,  8782,  8996,  9985],
        [ 9038,  9542,  9749, ...,  9386,  9170, 11951],
        [ 9176,  9472,  9467, ...,  9386,  9297, 14573],
        ...,
        [ 9193,  9164,  9030, ...,  8955,  9165,  8768],
        [ 9184,  9298,  9102, ...,  9074,  9062,  8891],
        [ 9297,  9502,  9239, ...,  9088,  9028,  8909]],

       [[13619, 13604, 13705, ..., 11812,  9562,  9898],
        [13835, 13754, 13633, ..., 11182,  9773,  9977],
        [14158, 13988, 13859, ..., 11182, 11055, 10513],
        ...,
        [11454, 11469, 11481, ..., 12696, 12741, 12744],
        [11214, 11323, 11381, ..., 12791, 12843, 12808],
        [10877, 11136, 11356, ..., 12887, 12939, 12823]]], dtype=uint16)
    • red
      (time, y, x)
      uint16
      9550 9115 9580 ... 13125 13015
      units :
      reflectance
      nodata :
      0
      crs :
      EPSG:32655
      grid_mapping :
      spatial_ref
      array([[[ 9550,  9115,  9580, ..., 23019, 23305, 23138],
        [ 8979,  9417,  9813, ..., 23452, 23724, 23755],
        [ 8838,  9016,  9031, ..., 23452, 23933, 24052],
        ...,
        [ 8443,  8458,  8422, ...,  9852,  9834,  9911],
        [ 8520,  8549,  8411, ...,  9802,  9810, 10007],
        [ 8493,  8560,  8708, ...,  9802,  9897, 10004]],

       [[ 8751,  8660,  9671, ..., 26347, 26357, 26617],
        [ 7850,  9980, 10718, ..., 21402, 21725, 18605],
        [ 9778, 11281, 10870, ..., 21402, 15289, 12893],
        ...,
        [22247, 22282, 22374, ...,  9667,  9959,  9705],
        [22244, 22239, 22322, ...,  9801,  9877,  9545],
        [22284, 22264, 22390, ...,  9359,  9724,  9497]],

       [[ 9044,  8792,  9352, ...,  7957,  8190,  9329],
        [ 8519,  9233,  9545, ...,  8405,  8149, 11198],
        [ 8561,  8966,  8969, ...,  8405,  8389, 14081],
        ...,
        [ 8416,  8387,  8342, ...,  8354,  8521,  8155],
        [ 8439,  8535,  8345, ...,  8434,  8415,  8293],
        [ 8371,  8723,  8790, ...,  8557,  8483,  8351]],

       [[13750, 13895, 13899, ..., 11380,  9657,  9952],
        [13852, 13826, 13776, ..., 11207, 10311, 10342],
        [14198, 14045, 13908, ..., 11207, 11553, 10877],
        ...,
        [10342, 10378, 10407, ..., 12803, 12928, 12898],
        [10279, 10315, 10294, ..., 12927, 13031, 13015],
        [10175, 10348, 10337, ..., 12978, 13125, 13015]]], dtype=uint16)
    • nir08
      (time, y, x)
      uint16
      21181 19938 19005 ... 20151 20301
      units :
      reflectance
      nodata :
      0
      crs :
      EPSG:32655
      grid_mapping :
      spatial_ref
      array([[[21181, 19938, 19005, ..., 29292, 29479, 29201],
        [21015, 19020, 19468, ..., 29585, 29663, 29533],
        [20673, 21338, 21813, ..., 29585, 29735, 29675],
        ...,
        [20889, 21218, 21211, ..., 18704, 19969, 18457],
        [21834, 21490, 20982, ..., 19083, 20432, 20543],
        [23107, 21159, 19469, ..., 19138, 20000, 20737]],

       [[22928, 23593, 21698, ..., 32762, 32504, 32767],
        [23889, 23243, 22184, ..., 29785, 30361, 28775],
        [23564, 23628, 23053, ..., 29785, 28442, 28198],
        ...,
        [26234, 26246, 26347, ..., 17690, 21168, 22110],
        [26263, 26232, 26334, ..., 20379, 21417, 22634],
        [26299, 26300, 26377, ..., 20973, 20792, 22630]],

       [[20984, 19738, 19169, ..., 20581, 21492, 22936],
        [21084, 19271, 19326, ..., 24725, 23242, 23532],
        [20304, 19863, 19880, ..., 24725, 24107, 25458],
        ...,
        [19810, 20553, 20279, ..., 19773, 20154, 19952],
        [21152, 21471, 20724, ..., 19329, 19686, 19349],
        [21992, 20755, 19142, ..., 17781, 18238, 19455]],

       [[23045, 23271, 23113, ..., 20336, 20553, 20528],
        [23102, 23141, 23024, ..., 21425, 21288, 20996],
        [23224, 23167, 23059, ..., 21425, 21947, 21508],
        ...,
        [23334, 23549, 23640, ..., 20326, 20263, 20303],
        [23784, 23540, 23482, ..., 20250, 20331, 20272],
        [24120, 23212, 23184, ..., 20020, 20151, 20301]]], dtype=uint16)
    • swir16
      (time, y, x)
      uint16
      15521 14502 15444 ... 11694 11635
      units :
      reflectance
      nodata :
      0
      crs :
      EPSG:32655
      grid_mapping :
      spatial_ref
      array([[[15521, 14502, 15444, ..., 19840, 19775, 19576],
        [13590, 14546, 15563, ..., 19776, 19705, 19624],
        [13359, 14495, 14692, ..., 19776, 19915, 19783],
        ...,
        [13523, 13591, 13602, ..., 12680, 13338, 12894],
        [14157, 13981, 13381, ..., 12382, 13017, 13313],
        [14309, 13595, 14010, ..., 12601, 13053, 13137]],

       [[15178, 14853, 15269, ..., 22502, 21791, 21643],
        [14994, 15612, 15925, ..., 20200, 20074, 18752],
        [15875, 16451, 16241, ..., 20200, 17536, 16539],
        ...,
        [12444, 12453, 12453, ..., 12373, 14094, 14014],
        [12441, 12451, 12456, ..., 13830, 13831, 13907],
        [12457, 12449, 12448, ..., 13792, 13720, 13799]],

       [[15805, 14744, 15966, ..., 11836, 12037, 14329],
        [13281, 14838, 16095, ..., 13186, 13084, 15428],
        [13535, 14919, 15138, ..., 13186, 13394, 17792],
        ...,
        [13783, 13972, 13705, ..., 12196, 12904, 12232],
        [14161, 14435, 13530, ..., 12149, 12407, 11971],
        [14242, 14645, 14243, ..., 12109, 11961, 11973]],

       [[16493, 16416, 16392, ..., 14507, 14592, 14471],
        [16655, 16586, 16499, ..., 14122, 14192, 14288],
        [16803, 16746, 16732, ..., 14122, 14235, 14173],
        ...,
        [12322, 12319, 12334, ..., 11593, 11686, 11668],
        [12430, 12361, 12245, ..., 11574, 11680, 11660],
        [12470, 12411, 12375, ..., 11572, 11694, 11635]]], dtype=uint16)
    • swir22
      (time, y, x)
      uint16
      11529 10801 11726 ... 11000 10982
      units :
      reflectance
      nodata :
      0
      crs :
      EPSG:32655
      grid_mapping :
      spatial_ref
      array([[[11529, 10801, 11726, ..., 15510, 15565, 15422],
        [10379, 11099, 11770, ..., 15523, 15555, 15442],
        [10106, 10672, 10746, ..., 15523, 15636, 15500],
        ...,
        [ 9760,  9790,  9778, ..., 10334, 10580, 10427],
        [10006,  9948,  9693, ..., 10227, 10402, 10552],
        [ 9963,  9944, 10424, ..., 10277, 10400, 10449]],

       [[11443, 11287, 11625, ..., 17731, 17311, 17186],
        [11480, 12484, 12609, ..., 16320, 16159, 15042],
        [12819, 13608, 13306, ..., 16320, 14002, 13087],
        ...,
        [12758, 12770, 12791, ..., 10269, 10971, 10865],
        [12758, 12771, 12766, ..., 10942, 10876, 10806],
        [12759, 12747, 12768, ..., 10653, 10742, 10559]],

       [[11398, 10736, 11939, ...,  8894,  8998, 10950],
        [ 9961, 11174, 11977, ...,  9452,  9496, 12336],
        [ 9790, 10754, 10921, ...,  9452,  9698, 14818],
        ...,
        [ 9876,  9948,  9738, ...,  9361,  9689,  9356],
        [ 9890, 10064,  9722, ...,  9309,  9383,  9133],
        [ 9958, 10381, 10246, ...,  9361,  9231,  9197]],

       [[15381, 15384, 15394, ..., 13194, 13263, 13317],
        [15566, 15533, 15468, ..., 12719, 12918, 13142],
        [15739, 15695, 15706, ..., 12719, 13331, 13207],
        ...,
        [10567, 10565, 10566, ..., 10923, 10962, 10967],
        [10608, 10582, 10524, ..., 10924, 10967, 10978],
        [10605, 10636, 10662, ..., 10941, 11000, 10982]]], dtype=uint16)
    • qa_pixel
      (time, y, x)
      uint16
      21824 21824 21824 ... 22280 22280
      units :
      bit_index
      nodata :
      1
      flags_definition :
      {'snow': {'bits': 5, 'values': {'0': 'not_high_confidence', '1': 'high_confidence'}}, 'clear': {'bits': 6, 'values': {'0': 'not_clear', '1': 'clear'}}, 'cloud': {'bits': 3, 'values': {'0': 'not_high_confidence', '1': 'high_confidence'}}, 'water': {'bits': 7, 'values': {'0': 'land_or_cloud', '1': 'water'}}, 'cirrus': {'bits': 2, 'values': {'0': 'not_high_confidence', '1': 'high_confidence'}}, 'nodata': {'bits': 0, 'values': {'0': False, '1': True}}, 'qa_pixel': {'bits': [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15], 'values': {'1': 'Fill', '2': 'Dilated Cloud', '4': 'Cirrus', '8': 'Cloud', '16': 'Cloud Shadow', '32': 'Snow', '64': 'Clear', '128': 'Water', '256': 'Cloud Confidence low bit', '512': 'Cloud Confidence high bit', '1024': 'Cloud Shadow Confidence low bit', '2048': 'Cloud Shadow Confidence high bit', '4096': 'Snow Ice Confidence low bit', '8192': 'Snow Ice Confidence high bit', '16384': 'Cirrus Confidence low bit', '32768': 'Cirrus Confidence high bit'}, 'description': 'Level 2 pixel quality'}, 'cloud_shadow': {'bits': 4, 'values': {'0': 'not_high_confidence', '1': 'high_confidence'}}, 'dilated_cloud': {'bits': 1, 'values': {'0': 'not_dilated', '1': 'dilated'}}, 'cloud_confidence': {'bits': [8, 9], 'values': {'0': 'none', '1': 'low', '2': 'medium', '3': 'high'}}, 'cirrus_confidence': {'bits': [14, 15], 'values': {'0': 'none', '1': 'low', '2': 'reserved', '3': 'high'}}, 'snow_ice_confidence': {'bits': [12, 13], 'values': {'0': 'none', '1': 'low', '2': 'reserved', '3': 'high'}}, 'cloud_shadow_confidence': {'bits': [10, 11], 'values': {'0': 'none', '1': 'low', '2': 'reserved', '3': 'high'}}}
      crs :
      EPSG:32655
      grid_mapping :
      spatial_ref
      array([[[21824, 21824, 21824, ..., 55052, 55052, 55052],
        [21824, 21824, 21824, ..., 55052, 55052, 55052],
        [21824, 21824, 21824, ..., 55052, 55052, 55052],
        ...,
        [21824, 21824, 21824, ..., 55052, 55052, 55052],
        [21824, 21824, 21824, ..., 55052, 55052, 55052],
        [21824, 21824, 21824, ..., 55052, 55052, 55052]],

       [[22280, 22280, 22280, ..., 55052, 55052, 55052],
        [22280, 22280, 22280, ..., 55052, 55052, 55052],
        [22280, 22280, 22280, ..., 55052, 55052, 55052],
        ...,
        [55052, 55052, 55052, ..., 22280, 24082, 23826],
        [55052, 55052, 55052, ..., 23826, 24082, 23826],
        [55052, 55052, 55052, ..., 23826, 23826, 23826]],

       [[21824, 21824, 21824, ..., 21762, 21762, 21762],
        [21824, 21824, 21824, ..., 21762, 21762, 22018],
        [21824, 21824, 21824, ..., 21762, 21762, 22280],
        ...,
        [21824, 21824, 21824, ..., 21824, 21824, 21824],
        [21824, 21824, 21824, ..., 21824, 21824, 21824],
        [21824, 21824, 21824, ..., 21824, 21824, 21824]],

       [[55052, 55052, 55052, ..., 55052, 55052, 55052],
        [55052, 55052, 55052, ..., 55052, 55052, 55052],
        [55052, 55052, 55052, ..., 55052, 55052, 55052],
        ...,
        [22280, 22280, 22280, ..., 22280, 22280, 22280],
        [22280, 22280, 22280, ..., 22280, 22280, 22280],
        [22280, 22280, 22280, ..., 22280, 22280, 22280]]], dtype=uint16)
    • qa_aerosol
      (time, y, x)
      uint8
      160 160 160 160 ... 224 224 224 224
      units :
      bit_index
      nodata :
      1
      flags_definition :
      {'water': {'bits': 2, 'values': {'0': 'not_water', '1': 'water'}}, 'nodata': {'bits': 0, 'values': {'0': False, '1': True}}, 'qa_aerosol': {'bits': [0, 1, 2, 3, 4, 5, 6, 7], 'values': {'1': 'Fill', '2': 'Valid aerosol retrieval', '4': 'Water', '8': 'Unused', '16': 'Unused', '32': 'Interpolated Aerosol', '64': 'Aerosol Level low bit', '128': 'Aerosol Level high bit'}, 'description': 'Aerosol quality assessment'}, 'aerosol_level': {'bits': [6, 7], 'values': {'0': 'climatology', '1': 'low', '2': 'medium', '3': 'high'}}, 'valid_retrieval': {'bits': 1, 'values': {'0': 'not_valid', '1': 'valid'}}, 'interp_retrieval': {'bits': 5, 'values': {'0': 'not_aerosol_interpolated', '1': 'aerosol_interpolated'}}}
      crs :
      EPSG:32655
      grid_mapping :
      spatial_ref
      array([[[160, 160, 160, ..., 160, 224, 224],
        [160, 160, 160, ..., 224, 224, 224],
        [130, 160, 160, ..., 224, 224, 192],
        ...,
        [ 96,  66,  96, ..., 224, 224, 224],
        [ 96,  96,  96, ..., 224, 224, 194],
        [ 96,  96,  96, ..., 224, 224, 224]],

       [[224, 224, 224, ..., 224, 224, 224],
        [224, 224, 224, ...,  96,  96, 224],
        [224, 194, 224, ...,  96, 224, 224],
        ...,
        [160, 160, 128, ..., 224, 224, 224],
        [160, 160, 160, ..., 194, 224, 224],
        [160, 160, 160, ..., 224, 224, 224]],

       [[ 96,  96,  96, ..., 160, 160, 224],
        [ 96,  96,  96, ..., 160, 160, 160],
        [ 96,  96,  66, ..., 160, 130, 160],
        ...,
        [ 66,  96,  96, ...,  96,  96,  96],
        [ 96,  96,  96, ...,  96,  66,  96],
        [ 96,  96,  96, ...,  96,  96,  96]],

       [[224, 224, 224, ..., 224, 224, 224],
        [224, 224, 224, ..., 224, 224, 224],
        [224, 194, 224, ..., 224, 224, 224],
        ...,
        [224, 224, 194, ..., 224, 224, 224],
        [224, 224, 224, ..., 194, 224, 224],
        [224, 224, 224, ..., 224, 224, 224]]], dtype=uint8)
    • qa_radsat
      (time, y, x)
      uint16
      0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0
      units :
      bit_index
      nodata :
      1
      flags_definition :
      {'qa_radsat': {'bits': [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11], 'values': {'1': 'Band 1 Data Saturation', '2': 'Band 2 Data Saturation', '4': 'Band 3 Data Saturation', '8': 'Band 4 Data Saturation', '16': 'Band 5 Data Saturation', '32': 'Band 6 Data Saturation', '64': 'Band 7 Data Saturation', '128': 'Unused', '256': 'Band 9 Data Saturation', '512': 'Unused', '1024': 'Unused', '2048': 'Terrain occlusion'}, 'description': 'Radiometric saturation'}, 'b1_saturation': {'bits': 0, 'values': {'0': 'no_saturation', '1': 'saturated_data'}}, 'b2_saturation': {'bits': 1, 'values': {'0': 'no_saturation', '1': 'saturated_data'}}, 'b3_saturation': {'bits': 2, 'values': {'0': 'no_saturation', '1': 'saturated_data'}}, 'b4_saturation': {'bits': 3, 'values': {'0': 'no_saturation', '1': 'saturated_data'}}, 'b5_saturation': {'bits': 4, 'values': {'0': 'no_saturation', '1': 'saturated_data'}}, 'b6_saturation': {'bits': 5, 'values': {'0': 'no_saturation', '1': 'saturated_data'}}, 'b7_saturation': {'bits': 6, 'values': {'0': 'no_saturation', '1': 'saturated_data'}}, 'b9_saturation': {'bits': 8, 'values': {'0': 'no_saturation', '1': 'saturated_data'}}, 'terrain_occlusion': {'bits': 11, 'values': {'0': 'no_terrain_occlusion', '1': 'terrain_occlusion'}}}
      crs :
      EPSG:32655
      grid_mapping :
      spatial_ref
      array([[[0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0],
        ...,
        [0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0]],

       [[0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0],
        ...,
        [0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0]],

       [[0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0],
        ...,
        [0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0]],

       [[0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0],
        ...,
        [0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0]]], dtype=uint16)
  • crs :
    EPSG:32655
    grid_mapping :
    spatial_ref

Plot data¶

Plot the Short-wave Infrared (SWIR) ("swir22") band data at each timestep.

The datacube.load() function does not apply missing values or scaling attributes. These are left to the user's discretion and requirements.

  • The robust option excludes outliers when calculating the colour limts for a more consistent result.
In [7]:
# Xarray operations
band = 'swir22'
data[band].plot(col="time", robust=True, col_wrap=4);

Masking and Scaling¶

Many Earth observation products include a "quality" array that can be used to filter the measurement arrays. For example, most quality layers include a "cloud" confidence quality flag that can be use to remove pixels affected by clouds from further analysis. In addition:

  • Data products usually include a nodata value that defines the "null" or "fill" value used in the array, and;
  • Some data products may also have "scaling" attributes.

More information on these techniques is covered in other notebooks. In particular, see the datasets/*ipynb notebooks for specific examples.

Landsat data requires a scale factor to be applied to convert the data to physical reflectance or temperature values. Once converted, Landsat surface reflectance values will have numbers ranging from 0 to 1 and surface temperature values will be in the units of degrees Kelvin. The scaling factor is different for different Landsat "Collections" and it is different for the Surface Reflectance and Surface Temperature products. See https://www.usgs.gov/faqs/how-do-i-use-a-scale-factor-landsat-level-2-science-products.

In the cell below we:

  1. Apply the nodata value
  2. Define and apply the Landsat "Collection 2" surface reflectance scaling values.
  3. Apply a cloud mask

Note the exclusion of cloud (if present) and the change in the colourbar range from the previous figure.

In [8]:
# Choose bands for further processing
bands = ['red', 'green', 'blue', 'nir08', 'swir22']

# Make a mask array for the nodata value
valid_mask = masking.valid_data_mask(data[bands])

# Define the scaling values (landsat8_c2l2_sr)
scale_factor = 0.0000275
add_offset = -0.2

# Make a scaled data array
scaled_data = data[bands] * scale_factor + add_offset
In [9]:
# Make a cloud mask (landsat8_c2l2_sr)
from datacube.utils import masking

# Optional - Show the flag_definition information
# See also http://explorer.asia.easi-eo.solutions/products/landsat8_c2l2_sr#definition-doc
# display( masking.describe_variable_flags(data.qa_pixel) )

# Multiple flags are combined as logical AND (bitwise)
cloud_mask = masking.make_mask(data['qa_pixel'], 
    clear='clear',
)
In [10]:
# Apply each of the masks
filtered_data = scaled_data.where(valid_mask & cloud_mask)

filtered_data['swir22'].plot(col="time", robust=True, col_wrap=4);

Perform a calculation on the data¶

As a simple example, we calculate the Normalized Difference Vegetation Index (NDVI).

In [11]:
# Calculate the NDVI
band_diff = filtered_data.nir08 - filtered_data.red
band_sum = filtered_data.nir08 + filtered_data.red
ndvi = band_diff / band_sum

# Plot the masked NDVI
ndvi.plot(col="time", robust=True, col_wrap=4, vmin=0, vmax=0.7, cmap='RdYlGn');

Save the results to file¶

After processing the data we can then save the output to a file that can then be imported into other applications for further analysis or publication if required.

The file will be saved to your home directory and appear on the File Browser panel to the left. You may need to select the folder icon to go to the top level (/home/jovyan/).

Download a file by 'right-click' Download.

In [12]:
# Xarray can save the data to a netCDF file

ndvi.time.attrs.pop('units', None)  # Xarray re-applies this
ndvi.to_netcdf("/home/jovyan/landsat8_sr_ndvi.nc")
In [13]:
# Or export to geotiff using rioxarray.

import rioxarray
ndvi.isel(time=0).rio.to_raster("/home/jovyan/landsat8_sr_ndvi.tif")  # Single time slices

Summary¶

This notebook introduced the main steps for querying data (with Datacube), and filtering, plotting, calculating and saving a "cube" of data (with Xarray).

There is plenty of detail and options to explore so please work through the other notebooks to learn more and refer back to these notebooks when required. We encourage you to create or bring your own notebooks, and adapt notebooks from other open-license repositories.

In [14]:
if rgb:
    rgb(filtered_data.isel(time=2), ['red', 'green', 'blue'])

Further reading¶

Open Data Cube¶

  • ODC documentation
  • ODC github

Python¶

Recommended level: Basic Python knowledge and familiarity with array manipulations, numpy and xarray. Familiarity with some plotting libraries (e.g., matplotlib) would also help.

There are many options for learning Python from online resources to in-house or facilitated training. Some examples are offered here with no suggestion that EASI endorses any of them.

  • https://www.python.org/about/gettingstarted
  • Learn Python tutorials: https://www.learnpython.org
  • Data Camp: https://www.datacamp.com
  • David Beazley courses: https://dabeaz-course.github.io/practical-python
  • Python Charmers: https://pythoncharmers.com

Background for selected libraries:

  • Numpy: https://numpy.org/doc/stable/user/quickstart.html
  • Xarray: http://xarray.pydata.org/en/stable/user-guide/data-structures.html
  • Xarray: https://towardsdatascience.com/basic-data-structures-of-xarray-80bab8094efa

JupyterLab¶

Recommended level: Familiarity with notebooks.

The JupyterLab website has excellent documentation including video instructions. We recommend users take a few minutes to orientate themselves with the use and features of JupyterLab.

  • Getting started: https://jupyterlab.readthedocs.io/en/stable/getting_started/overview.html
  • Drag and drop upload of files: https://jupyterlab.readthedocs.io/en/stable/user/files.html

Git¶

Recommended level: Basic understanding of concepts such as clone, add, commit and push would help.

Git is a document version control system. It retains a full history of changes to all files (including deleted ones) by tracking incremental changes and recording a history timeline of changes. Changes you make append to the history timeline. Git allows you to copy ("clone") a repository, make changes to files, and "commit" and "push" these changes back to the source repository.

The best way to learn Git is by practice and incrementally; start with simple, common actions and gain more knowledge when required. Some useful Git links are:

  • Getting started: https://git-scm.com/doc
  • Using the JupyterLab Git extension: https://annefou.github.io/jupyter_publish/02-git/index.html
  • DEA Git guide: https://github.com/GeoscienceAustralia/dea-notebooks/wiki/Guide-to-using-DEA-Notebooks-with-git
  • Undoing things guide: https://git-scm.com/book/en/v2/Git-Basics-Undoing-Things
  • Understanding branches: https://nvie.com/posts/a-successful-git-branching-model
In [ ]: