Landsat TIRS

Notebook showing a fast and efficient way to load Sentinel-2 L1C data into xarray

import glob
import os
import numpy as np

import pandas as pd
import geopandas as gpd
import xarray as xr
import rioxarray

import matplotlib.pyplot as plt
import matplotlib as mpl
import cartopy.crs as ccrs
from matplotlib.colors import ListedColormap

import GRSdriver
from GRSdriver import visual

opj = os.path.join

print(f'-GRSdriver: {GRSdriver.__version__}')

-GRSdriver: 1.0.5

Set the path of the images you want to load, here is a example for the former (v3) and current (v4) image formatting

image ='/data/satellite/landsat/LC09_L1TP_196030_20221120_20230321_02_T1'

l1c = GRSdriver.LandsatDriver(image,band_idx=[9,10])

l1c.INFO

bandId	9	10
ESA	B08	B09
EOREADER	TIR_1	TIR_2
Wavelength (nm)	11000	12000
Band width (nm)	590	1010
Resolution (m)	100	100

l1c.load_mask()

l1c.load_product()

l1c.prod

<xarray.Dataset>
Dimensions:      (x: 7671, y: 7791, wl: 2, wl_hr: 2951)
Coordinates:
  * x            (x) float64 5.349e+05 5.349e+05 5.35e+05 ... 7.65e+05 7.65e+05
  * y            (y) float64 4.899e+06 4.899e+06 ... 4.666e+06 4.666e+06
    spatial_ref  int64 0
  * wl           (wl) object 11000 12000
  * wl_hr        (wl_hr) float64 400.0 401.0 402.0 ... 1.35e+04 1.35e+04
    time         datetime64[ns] 2022-11-20T10:24:21
    band         int64 1
Data variables:
    bands        (wl, y, x) float32 nan nan nan nan nan ... nan nan nan nan nan
    SRF          (wl, wl_hr) float64 ...
    wl_true      (wl) float64 1.084e+04 1.203e+04
    sza          (y, x) float32 nan nan nan nan nan nan ... nan nan nan nan nan
    vza          (y, x) float32 nan nan nan nan nan nan ... nan nan nan nan nan
    raa          (y, x) float32 nan nan nan nan nan nan ... nan nan nan nan nan
    flags_l1     (y, x) uint32 1 1 1 1 1 1 1 1 1 1 1 1 ... 1 1 1 1 1 1 1 1 1 1 1
Attributes: (12/17)
    long_name:              TIR_1 TIR_2
    constellation:          Landsat-9
    constellation_id:       L9
    product_path:           /data/satellite/landsat/LC09_L1TP_196030_20221120...
    product_name:           LC09_L1TP_196030_20221120_20230321_02_T1
    product_filename:       LC09_L1TP_196030_20221120_20230321_02_T1
    ...                     ...
    radiometry:             brightness temperature
    cloud_cover:            0.2
    wl_to_process:          [ 443  490  560  590  665  865 1610 2190]
    satellite:              Landsat-9
    solar_irradiance:       NA
    solar_irradiance_unit:  W/m²/µm

xarray.Dataset

• Dimensions:
  • x: 7671
  • y: 7791
  • wl: 2
  • wl_hr: 2951
• Coordinates: (7)
  • x
    (x)
    float64
    5.349e+05 5.349e+05 ... 7.65e+05

    array([534900., 534930., 534960., ..., 764940., 764970., 765000.])

  • y
    (y)
    float64
    4.899e+06 4.899e+06 ... 4.666e+06

    array([4899300., 4899270., 4899240., ..., 4665660., 4665630., 4665600.])

  • spatial_ref
    ()
    int64
    0

    spatial_ref :

    PROJCS["WGS 84 / UTM zone 31N",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",3],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","32631"]]

    crs_wkt :

    PROJCS["WGS 84 / UTM zone 31N",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",3],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","32631"]]

    GeoTransform :

    534885.0 30.0 0.0 4899315.0 0.0 -30.0

    array(0)

  • wl
    (wl)
    object
    11000 12000

    array([11000, 12000], dtype=object)

  • wl_hr
    (wl_hr)
    float64
    400.0 401.0 ... 1.35e+04 1.35e+04

    array([  400.,   401.,   402., ..., 13492., 13496., 13500.])

  • time
    ()
    datetime64[ns]
    2022-11-20T10:24:21

    array('2022-11-20T10:24:21.000000000', dtype='datetime64[ns]')

  • band
    ()
    int64
    1

    array(1)

• Data variables: (7)
  • bands
    (wl, y, x)
    float32
    nan nan nan nan ... nan nan 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, ..., nan, nan, nan],
            [nan, nan, nan, ..., nan, nan, nan],
            [nan, nan, nan, ..., nan, nan, nan],
            ...,
            [nan, nan, nan, ..., nan, nan, nan],
            [nan, nan, nan, ..., nan, nan, nan],
            [nan, nan, nan, ..., nan, nan, nan]]], dtype=float32)

  • SRF
    (wl, wl_hr)
    float64
    ...

    OLI_long_name :

    SRF_Landat_9_OLI

    OLI_file_creation :

    2024-04-12 12:33:43.389548

    OLI_from_file :

    L9_OLI2_Ball_FPM_RSR.v1.0.xlsx

    OLI_description :

    relative spectral response functions for Landat-9 OLI bands

    TIRS_long_name :

    SRF_Landat_9_TIRS

    TIRS_file_creation :

    2024-04-12 12:33:43.464885

    TIRS_from_file :

    L9_TIRS2_Relative_Spectral_Responses.BA.v1.0.xlsx

    TIRS_description :

    relative spectral response functions for Landat-9 TIRS bands

    [5902 values with dtype=float64]

  • wl_true
    (wl)
    float64
    1.084e+04 1.203e+04

    array([10837.17325635, 12025.32314819])

  • sza
    (y, x)
    float32
    nan nan nan nan ... nan nan 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]], dtype=float32)

  • vza
    (y, x)
    float32
    nan nan nan nan ... nan nan 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]], dtype=float32)

  • raa
    (y, x)
    float32
    nan nan nan nan ... nan nan 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]], dtype=float32)

  • flags_l1
    (y, x)
    uint32
    1 1 1 1 1 1 1 1 ... 1 1 1 1 1 1 1 1

    long_name :

    RAW_CLOUDS

    constellation :

    Landsat-9

    constellation_id :

    L9

    product_path :

    /data/satellite/landsat/LC09_L1TP_196030_20221120_20230321_02_T1

    product_name :

    LC09_L1TP_196030_20221120_20230321_02_T1

    product_filename :

    LC09_L1TP_196030_20221120_20230321_02_T1

    instrument :

    OLI-TIRS

    product_type :

    L1

    acquisition_date :

    20221120T102421

    condensed_name :

    20221120T102421_L9_196030_OLI_TIRS

    orbit_direction :

    DESCENDING

    array([[1, 1, 1, ..., 1, 1, 1],
           [1, 1, 1, ..., 1, 1, 1],
           [1, 1, 1, ..., 1, 1, 1],
           ...,
           [1, 1, 1, ..., 1, 1, 1],
           [1, 1, 1, ..., 1, 1, 1],
           [1, 1, 1, ..., 1, 1, 1]], dtype=uint32)

• Indexes: (4)
  • x
    PandasIndex

    PandasIndex(Index([534900.0, 534930.0, 534960.0, 534990.0, 535020.0, 535050.0, 535080.0,
           535110.0, 535140.0, 535170.0,
           ...
           764730.0, 764760.0, 764790.0, 764820.0, 764850.0, 764880.0, 764910.0,
           764940.0, 764970.0, 765000.0],
          dtype='float64', name='x', length=7671))

  • y
    PandasIndex

    PandasIndex(Index([4899300.0, 4899270.0, 4899240.0, 4899210.0, 4899180.0, 4899150.0,
           4899120.0, 4899090.0, 4899060.0, 4899030.0,
           ...
           4665870.0, 4665840.0, 4665810.0, 4665780.0, 4665750.0, 4665720.0,
           4665690.0, 4665660.0, 4665630.0, 4665600.0],
          dtype='float64', name='y', length=7791))

  • wl
    PandasIndex

    PandasIndex(Index([11000, 12000], dtype='object', name='wl'))

  • wl_hr
    PandasIndex

    PandasIndex(Index([  400.0,   401.0,   402.0,   403.0,   404.0,   405.0,   406.0,   407.0,
             408.0,   409.0,
           ...
           13464.0, 13468.0, 13472.0, 13476.0, 13480.0, 13484.0, 13488.0, 13492.0,
           13496.0, 13500.0],
          dtype='float64', name='wl_hr', length=2951))

• Attributes: (17)

  long_name :

  TIR_1 TIR_2

  constellation :

  Landsat-9

  constellation_id :

  L9

  product_path :

  /data/satellite/landsat/LC09_L1TP_196030_20221120_20230321_02_T1

  product_name :

  LC09_L1TP_196030_20221120_20230321_02_T1

  product_filename :

  LC09_L1TP_196030_20221120_20230321_02_T1

  instrument :

  OLI-TIRS

  product_type :

  L1

  acquisition_date :

  2022-11-20 10:24:21

  condensed_name :

  20221120T102421_L9_196030_OLI_TIRS

  orbit_direction :

  DESCENDING

  radiometry :

  brightness temperature

  cloud_cover :

  0.2

  wl_to_process :

  [ 443 490 560 590 665 865 1610 2190]

  satellite :

  Landsat-9

  solar_irradiance :

  NA

  solar_irradiance_unit :

  W/m²/µm

#xr.open_dataset('/DATA/git/satellite_app/GRSdriver/GRSdriver/rsr/data/rsr_landsat_8_oli_tirs.nc').SRF.plot(x='wl_hr',hue='wl',lw=1)

import cycler

n = 12
color = plt.cm.Spectral_r(np.linspace(0, 1,n))
mpl.rcParams['axes.prop_cycle'] = cycler.cycler('color', color)

l1c.prod.SRF.plot(x='wl_hr',hue='wl',lw=1)

[<matplotlib.lines.Line2D at 0x7f3f20b8d430>,
 <matplotlib.lines.Line2D at 0x7f3f20b8d490>]

h = visual.ViewSpectral(l1c.prod.bands.assign_coords(time=l1c.datetime).expand_dims('time'),reproject=True)

h.visu()

import cartopy.crs as ccrs

epsg = l1c.extent.crs.to_epsg()
str_epsg = str(epsg)
zone = str_epsg[-2:]
is_south = str_epsg[2] == 7
proj = ccrs.UTM(zone, is_south)

coarsening=1
cmap = plt.cm.Spectral_r
fig = plt.figure(figsize=(20, 15))
ax = plt.subplot(1, 1, 1, projection=proj)
ax.set_extent(extent_val, proj)
l1c.prod.bands[:, ::coarsening, ::coarsening].isel(wl=0).plot.imshow( cmap=cmap,robust=True,subplot_kws=dict(projection=proj))
ax.set_title('')

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[13], line 5
      3 fig = plt.figure(figsize=(20, 15))
      4 ax = plt.subplot(1, 1, 1, projection=proj)
----> 5 ax.set_extent(extent_val, proj)
      6 l1c.prod.bands[:, ::coarsening, ::coarsening].isel(wl=0).plot.imshow( cmap=cmap,robust=True,subplot_kws=dict(projection=proj))
      7 ax.set_title('')

NameError: name 'extent_val' is not defined