import os, copy
import glob

import numpy as np
import pandas as pd
import xarray as xr

import datetime as dt

#import netCDF4
#import h5py

import cartopy.crs as ccrs

import matplotlib as mpl
import matplotlib.pyplot as plt
#import colorcet as cc

#import panel as pn

#pn.extension()
opj = os.path.join

file='/data/satellite/enmap/L2A/Trasimeno/ENMAP01-____L1C-DT0000088121_20240817T102235Z_002_V010502_20250128T220420Z.nc'
file='/data/satellite/enmap/L2A/Trasimeno/ENMAP01-____L1C-DT0000082532_20240706T103357Z_002_V010502_20250128T220612Z.nc' #8121_20240817T102235Z_002_V010502_20250128T220420Z.nc'

img =xr.open_dataset(file)


img.load()
l2_prod = img
l2_prod

<xarray.Dataset> Size: 1GB
Dimensions:           (wl: 92, y: 1221, x: 1295, yc: 62, xc: 65)
Coordinates:
  * x                 (x) float64 10kB 2.47e+05 2.47e+05 ... 2.858e+05 2.858e+05
    spatial_ref       int64 8B 0
    time              datetime64[ns] 8B 2024-07-06T10:33:53.032721
  * y                 (y) float64 10kB 4.796e+06 4.796e+06 ... 4.76e+06 4.76e+06
  * wl                (wl) float64 736B 418.4 424.0 ... 1.093e+03 1.105e+03
  * xc                (xc) float64 520B 2.473e+05 2.479e+05 ... 2.856e+05
  * yc                (yc) float64 496B 4.796e+06 4.795e+06 ... 4.76e+06
Data variables: (12/13)
    Rrs               (wl, y, x) float64 1GB nan nan nan nan ... nan nan nan nan
    brdfg_full        (y, x) float64 13MB nan nan nan nan ... nan nan nan nan
    tcwv              (yc, xc) float64 32kB nan nan nan nan ... nan nan nan nan
    tcwv_std          (yc, xc) float64 32kB nan nan nan nan ... nan nan nan nan
    aot_ref           (yc, xc) float64 32kB nan nan nan nan ... nan nan nan nan
    brdfg             (yc, xc) float64 32kB nan nan nan nan ... nan nan nan nan
    ...                ...
    brdfg_std         (yc, xc) float64 32kB nan nan nan nan ... nan nan nan nan
    aot_ref_smoothed  (yc, xc) float64 32kB nan nan nan nan ... nan nan nan nan
    water_pix_prop    (yc, xc) float64 32kB 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0
    pressure          float64 8B 972.4
    to3c              float64 8B 0.006824
    tno2c             float64 8B 3.605e-06
Attributes: (12/41)
    description:                     PRISMA L2A-hGRS cube data
    L1C_product_name:                ENMAP01-____L1C-DT0000082532_20240706T10...
    acquisition_date:                2024-07-06 10:33:53.032721
    vnir_index:                      [  0   1   2   3   4   5   6   7   8   9...
    swir_index:                      [ 79  81  84  86  88  90  93  95  97  99...
    vnir_index_tokeep:               [ 0  1  2  3  4  5  6  7  8  9 10 11 12 ...
    ...                              ...
    to3c:                            0.006824255455285311
    tno2c:                           3.6045623801328475e-06
    tch4c:                           0.0098807318136096
    psl:                             1013
    coef_abs_scat:                   1.0
    altitude:                        0

xarray.Dataset

• Dimensions:
  • wl: 92
  • y: 1221
  • x: 1295
  • yc: 62
  • xc: 65
• Coordinates: (7)
  • x
    (x)
    float64
    2.47e+05 2.47e+05 ... 2.858e+05

    array([246990., 247020., 247050., ..., 285750., 285780., 285810.],
          shape=(1295,))

  • spatial_ref
    ()
    int64
    0

    crs_wkt :

    PROJCS["WGS 84 / UTM zone 33N",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",15],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","32633"]]

    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 33N

    grid_mapping_name :

    transverse_mercator

    latitude_of_projection_origin :

    0.0

    longitude_of_central_meridian :

    15.0

    false_easting :

    500000.0

    false_northing :

    0.0

    scale_factor_at_central_meridian :

    0.9996

    spatial_ref :

    PROJCS["WGS 84 / UTM zone 33N",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",15],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","32633"]]

    GeoTransform :

    246975.0 30.0 0.0 4796205.0 0.0 -30.0

    array(0)

  • time
    ()
    datetime64[ns]
    2024-07-06T10:33:53.032721

    array('2024-07-06T10:33:53.032721000', dtype='datetime64[ns]')

  • y
    (y)
    float64
    4.796e+06 4.796e+06 ... 4.76e+06

    array([4796190., 4796160., 4796130., ..., 4759650., 4759620., 4759590.],
          shape=(1221,))

  • wl
    (wl)
    float64
    418.4 424.0 ... 1.093e+03 1.105e+03

    array([ 418.416,  424.043,  429.457,  434.686,  439.757,  444.699,  449.539,
            454.306,  459.031,  463.73 ,  468.411,  473.08 ,  477.744,  482.411,
            487.087,  491.78 ,  496.497,  501.243,  506.02 ,  510.828,  515.672,
            520.55 ,  525.467,  530.424,  535.422,  540.463,  545.551,  550.687,
            555.873,  561.112,  566.405,  571.756,  577.166,  582.636,  588.171,
            593.773,  599.446,  605.193,  611.017,  616.923,  622.92 ,  628.987,
            635.112,  641.294,  647.537,  653.841,  660.207,  666.637,  673.131,
            679.691,  686.319,  693.014,  699.78 ,  706.617,  713.524,  720.501,
            727.545,  734.654,  741.826,  749.06 ,  756.353,  763.703,  771.108,
            778.567,  786.078,  793.639,  801.248,  808.905,  816.608,  824.355,
            832.145,  839.976,  847.847,  855.757,  863.703,  871.683,  879.692,
            887.729,  931.204,  972.017,  982.524,  993.145, 1003.88 , 1014.72 ,
           1025.66 , 1036.7  , 1047.84 , 1059.07 , 1070.39 , 1081.79 , 1093.26 ,
           1104.81 ])

  • xc
    (xc)
    float64
    2.473e+05 2.479e+05 ... 2.856e+05

    array([247275., 247875., 248475., 249075., 249675., 250275., 250875., 251475.,
           252075., 252675., 253275., 253875., 254475., 255075., 255675., 256275.,
           256875., 257475., 258075., 258675., 259275., 259875., 260475., 261075.,
           261675., 262275., 262875., 263475., 264075., 264675., 265275., 265875.,
           266475., 267075., 267675., 268275., 268875., 269475., 270075., 270675.,
           271275., 271875., 272475., 273075., 273675., 274275., 274875., 275475.,
           276075., 276675., 277275., 277875., 278475., 279075., 279675., 280275.,
           280875., 281475., 282075., 282675., 283275., 283875., 284475., 285075.,
           285600.])

  • yc
    (yc)
    float64
    4.796e+06 4.795e+06 ... 4.76e+06

    array([4795905., 4795305., 4794705., 4794105., 4793505., 4792905., 4792305.,
           4791705., 4791105., 4790505., 4789905., 4789305., 4788705., 4788105.,
           4787505., 4786905., 4786305., 4785705., 4785105., 4784505., 4783905.,
           4783305., 4782705., 4782105., 4781505., 4780905., 4780305., 4779705.,
           4779105., 4778505., 4777905., 4777305., 4776705., 4776105., 4775505.,
           4774905., 4774305., 4773705., 4773105., 4772505., 4771905., 4771305.,
           4770705., 4770105., 4769505., 4768905., 4768305., 4767705., 4767105.,
           4766505., 4765905., 4765305., 4764705., 4764105., 4763505., 4762905.,
           4762305., 4761705., 4761105., 4760505., 4759905., 4759590.])

• Data variables: (13)
  • Rrs
    (wl, y, x)
    float64
    nan nan nan nan ... nan nan nan nan

    unit :

    per steradian

    long_name :

    Remote sensing reflectance

    description :

    Directional water-leaving radiance normalized by downwelling irradiance in the observation geometry

    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]],
    
           [[nan, nan, nan, ..., nan, nan, nan],
            [nan, nan, nan, ..., nan, nan, nan],
            [nan, nan, nan, ..., nan, nan, nan],
            ...,
    ...
            ...,
            [nan, nan, nan, ..., nan, nan, nan],
            [nan, nan, nan, ..., nan, nan, nan],
            [nan, nan, nan, ..., nan, nan, nan]],
    
           [[nan, nan, nan, ..., nan, nan, nan],
            [nan, nan, nan, ..., nan, nan, nan],
            [nan, nan, nan, ..., nan, nan, nan],
            ...,
            [nan, nan, nan, ..., nan, nan, nan],
            [nan, nan, nan, ..., nan, nan, nan],
            [nan, nan, nan, ..., nan, nan, nan]],
    
           [[nan, nan, nan, ..., nan, nan, nan],
            [nan, nan, nan, ..., nan, nan, nan],
            [nan, nan, nan, ..., nan, nan, nan],
            ...,
            [nan, nan, nan, ..., nan, nan, nan],
            [nan, nan, nan, ..., nan, nan, nan],
            [nan, nan, nan, ..., nan, nan, nan]]], shape=(92, 1221, 1295))

  • brdfg_full
    (y, x)
    float64
    nan nan nan nan ... nan nan nan nan

    unit :

    -

    long_name :

    BRDF_sunglint

    description :

    Bidirectional reflectance distribution function estimated from the sunglint in the SWIR for the observation geometry

    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]], shape=(1221, 1295))

  • tcwv
    (yc, xc)
    float64
    nan nan nan nan ... nan nan nan nan

    unit :

    kg m-2

    long_name :

    total_columnar_water_vapor

    description :

    Water vapor integrated over the atmospheric layer

    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]], shape=(62, 65))

  • tcwv_std
    (yc, xc)
    float64
    nan nan nan nan ... nan nan nan nan

    unit :

    kg m-2

    long_name :

    total_columnar_water_vapor_standard_deviation

    description :

    Uncertainty based on optimal estimation procedure

    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]], shape=(62, 65))

  • aot_ref
    (yc, xc)
    float64
    nan nan nan nan ... nan nan nan nan

    unit :

    -

    long_name :

    aerosol_optical_thickness

    description :

    Aerosol optical thickness at the reference wavelength (550nm)

    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]], shape=(62, 65))

  • brdfg
    (yc, xc)
    float64
    nan nan nan nan ... nan nan nan nan

    unit :

    -

    long_name :

    BRDF_sunglint

    description :

    Bidirectional reflectance distribution function estimated from the sunglint in the SWIR for the observation geometry

    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]], shape=(62, 65))

  • aot_ref_std
    (yc, xc)
    float64
    nan nan nan nan ... nan nan nan nan

    unit :

    -

    long_name :

    aerosol_optical_thickness_standard_deviation

    description :

    Uncertainty based on optimal estimation procedure

    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]], shape=(62, 65))

  • brdfg_std
    (yc, xc)
    float64
    nan nan nan nan ... nan nan nan nan

    unit :

    -

    long_name :

    BRDF_sunglint_standard deviation

    description :

    Uncertainty based on optimal estimation procedure

    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]], shape=(62, 65))

  • aot_ref_smoothed
    (yc, xc)
    float64
    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]], shape=(62, 65))

  • water_pix_prop
    (yc, xc)
    float64
    0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0

    unit :

    -

    description :

    Relative number of water pixel within mega-pixel used for inversion

    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.]], shape=(62, 65))

  • pressure
    ()
    float64
    972.4

    unit :

    hPa

    description :

    Atmospheric pressure at the surface level

    source :

    computed from CAMS and DEM (see DEM metadata)

    array(972.37148438)

  • to3c
    ()
    float64
    0.006824

    unit :

    description :

    Total columnar ozone concentration

    source :

    CAMS

    array(0.00682426)

  • tno2c
    ()
    float64
    3.605e-06

    unit :

    description :

    Total columnar Nitrogen dioxide concentration

    source :

    CAMS

    array(3.60456238e-06)

• Indexes: (5)
  • x
    PandasIndex

    PandasIndex(Index([246990.0, 247020.0, 247050.0, 247080.0, 247110.0, 247140.0, 247170.0,
           247200.0, 247230.0, 247260.0,
           ...
           285540.0, 285570.0, 285600.0, 285630.0, 285660.0, 285690.0, 285720.0,
           285750.0, 285780.0, 285810.0],
          dtype='float64', name='x', length=1295))

  • y
    PandasIndex

    PandasIndex(Index([4796190.0, 4796160.0, 4796130.0, 4796100.0, 4796070.0, 4796040.0,
           4796010.0, 4795980.0, 4795950.0, 4795920.0,
           ...
           4759860.0, 4759830.0, 4759800.0, 4759770.0, 4759740.0, 4759710.0,
           4759680.0, 4759650.0, 4759620.0, 4759590.0],
          dtype='float64', name='y', length=1221))

  • wl
    PandasIndex

    PandasIndex(Index([418.416, 424.043, 429.457, 434.686, 439.757, 444.699, 449.539, 454.306,
           459.031,  463.73, 468.411,  473.08, 477.744, 482.411, 487.087,  491.78,
           496.497, 501.243,  506.02, 510.828, 515.672,  520.55, 525.467, 530.424,
           535.422, 540.463, 545.551, 550.687, 555.873, 561.112, 566.405, 571.756,
           577.166, 582.636, 588.171, 593.773, 599.446, 605.193, 611.017, 616.923,
            622.92, 628.987, 635.112, 641.294, 647.537, 653.841, 660.207, 666.637,
           673.131, 679.691, 686.319, 693.014,  699.78, 706.617, 713.524, 720.501,
           727.545, 734.654, 741.826,  749.06, 756.353, 763.703, 771.108, 778.567,
           786.078, 793.639, 801.248, 808.905, 816.608, 824.355, 832.145, 839.976,
           847.847, 855.757, 863.703, 871.683, 879.692, 887.729, 931.204, 972.017,
           982.524, 993.145, 1003.88, 1014.72, 1025.66,  1036.7, 1047.84, 1059.07,
           1070.39, 1081.79, 1093.26, 1104.81],
          dtype='float64', name='wl'))

  • xc
    PandasIndex

    PandasIndex(Index([247275.0, 247875.0, 248475.0, 249075.0, 249675.0, 250275.0, 250875.0,
           251475.0, 252075.0, 252675.0, 253275.0, 253875.0, 254475.0, 255075.0,
           255675.0, 256275.0, 256875.0, 257475.0, 258075.0, 258675.0, 259275.0,
           259875.0, 260475.0, 261075.0, 261675.0, 262275.0, 262875.0, 263475.0,
           264075.0, 264675.0, 265275.0, 265875.0, 266475.0, 267075.0, 267675.0,
           268275.0, 268875.0, 269475.0, 270075.0, 270675.0, 271275.0, 271875.0,
           272475.0, 273075.0, 273675.0, 274275.0, 274875.0, 275475.0, 276075.0,
           276675.0, 277275.0, 277875.0, 278475.0, 279075.0, 279675.0, 280275.0,
           280875.0, 281475.0, 282075.0, 282675.0, 283275.0, 283875.0, 284475.0,
           285075.0, 285600.0],
          dtype='float64', name='xc'))

  • yc
    PandasIndex

    PandasIndex(Index([4795905.0, 4795305.0, 4794705.0, 4794105.0, 4793505.0, 4792905.0,
           4792305.0, 4791705.0, 4791105.0, 4790505.0, 4789905.0, 4789305.0,
           4788705.0, 4788105.0, 4787505.0, 4786905.0, 4786305.0, 4785705.0,
           4785105.0, 4784505.0, 4783905.0, 4783305.0, 4782705.0, 4782105.0,
           4781505.0, 4780905.0, 4780305.0, 4779705.0, 4779105.0, 4778505.0,
           4777905.0, 4777305.0, 4776705.0, 4776105.0, 4775505.0, 4774905.0,
           4774305.0, 4773705.0, 4773105.0, 4772505.0, 4771905.0, 4771305.0,
           4770705.0, 4770105.0, 4769505.0, 4768905.0, 4768305.0, 4767705.0,
           4767105.0, 4766505.0, 4765905.0, 4765305.0, 4764705.0, 4764105.0,
           4763505.0, 4762905.0, 4762305.0, 4761705.0, 4761105.0, 4760505.0,
           4759905.0, 4759590.0],
          dtype='float64', name='yc'))

• Attributes: (41)

  description :

  PRISMA L2A-hGRS cube data

  L1C_product_name :

  ENMAP01-____L1C-DT0000082532_20240706T103357Z_002_V010502_20250128T220612Z

  acquisition_date :

  2024-07-06 10:33:53.032721

  vnir_index :

  [ 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 80 82 83 85 87 89 91 92 94 96 98 100]

  swir_index :

  [ 79 81 84 86 88 90 93 95 97 99 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223]

  vnir_index_tokeep :

  [ 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77]

  swir_index_tokeep :

  [ 86 88 90 93 95 97 99 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223]

  vnir_bands :

  [418.416 424.043 429.457 434.686 439.757 444.699 449.539 454.306 459.031 463.73 468.411 473.08 477.744 482.411 487.087 491.78 496.497 501.243 506.02 510.828 515.672 520.55 525.467 530.424 535.422 540.463 545.551 550.687 555.873 561.112 566.405 571.756 577.166 582.636 588.171 593.773 599.446 605.193 611.017 616.923 622.92 628.987 635.112 641.294 647.537 653.841 660.207 666.637 673.131 679.691 686.319 693.014 699.78 706.617 713.524 720.501 727.545 734.654 741.826 749.06 756.353 763.703 771.108 778.567 786.078 793.639 801.248 808.905 816.608 824.355 832.145 839.976 847.847 855.757 863.703 871.683 879.692 887.729]

  swir_bands :

  [ 931.204 941.218 951.361 961.629 972.017 982.524 993.145 1003.88 1014.72 1025.66 1036.7 1047.84 1059.07 1070.39 1081.79 1093.26 1104.81 1116.43 1128.11 1139.84 1151.62 1163.44 1175.31 1187.2 1199.11 1211.05 1223. 1234.97 1246.95 1258.93 1270.92 1282.92 1294.91 1306.9 1318.88 1330.86 1342.82 1354.76 1366.69 1378.6 1390.48 1449.43 1461.11 1472.74 1484.34 1495.89 1507.4 1518.87 1530.29 1541.68 1553.01 1564.31 1575.55 1586.76 1597.91 1609.02 1620.09 1631.11 1642.08 1653. 1663.87 1674.7 1685.47 1696.2 1706.87 1717.5 1728.08 1738.6 1749.08 1759.51 1769.89 1780.22 1967.66 1977.08 1986.45 1995.79 2005.08 2014.33 2023.54 2032.71 2041.83 2050.92 2059.96 2068.97 2077.93 2086.86 2095.74 2104.59 2113.4 2122.17 2130.9 2139.6 2148.26 2156.88 2165.47 2174.02 2182.53 2191.01 2199.46 2207.86 2216.24 2224.58 2232.89 2241.16 2249.4 2257.61 2265.79 2273.93 2282.04 2290.12 2298.17 2306.19 2314.18 2322.13 2330.05 2337.94 2345.81 2353.64 2361.44 2369.21 2376.95 2384.66 2392.34 2400. 2407.62 2415.21 2422.78 2430.32 2437.83 2445.3 ]

  processing_date :

  2025-11-25 00:06:27.144028

  hgrs_version :

  1.1.0

  DEM :

  not available

  aerosol_model :

  COAV_rh70

  wl_water_vapor :

  slice(800, 1300, None)

  wl_sunglint :

  slice(2150, 2250, None)

  wl_atmo :

  [1000, 1050, 1075, 1100, 1200, 1300, 1600, 1650, 1700, 2150, 2200, 2250]

  wl_to_remove :

  [(935, 967), (1105, 1170), (1320, 1490), (1778, 2033), (2465, 2550)]

  wl_green :

  slice(540, 570, None)

  wl_nir :

  slice(850, 882, None)

  wl_1600 :

  slice(1580, 1650, None)

  wl_rgb :

  [30, 20, 10]

  xcoarsen :

  20

  ycoarsen :

  20

  Npix_per_megapix :

  400

  block_size :

  2

  pixel_percentage :

  20

  pixel_threshold :

  80.0

  ang_resol :

  1

  abs_gas_file :

  /DATA/git/satellite_app/hgrs/data/lut/lut_abs_opt_thickness_normalized.nc

  lut_file :

  /DATA/git/satellite_app/hgrs/data/lut/toa_lut_opac_wind_up_v3.nc

  water_vapor_transmittance_file :

  /DATA/git/satellite_app/hgrs/data/lut/water_vapor_transmittance.nc

  sunglint_threshold :

  0.11

  ndwi_threshold :

  0.01

  green_swir_index_threshold :

  0.1

  pressure :

  972.371484375

  to3c :

  0.006824255455285311

  tno2c :

  3.6045623801328475e-06

  tch4c :

  0.0098807318136096

  psl :

  1013

  coef_abs_scat :

  1.0

  altitude :

  0

# get geographic information 
str_epsg = str(l2_prod.rio.crs)
zone = str_epsg[-2:]
is_south = str_epsg[2] == 7
proj = ccrs.UTM(zone, is_south)        

ncols=2
nrows=2
fig,axs = plt.subplots(nrows,ncols,figsize=(ncols*5.,5*nrows),subplot_kw={'projection': proj}) 
fig.subplots_adjust(bottom=0.08, top=0.9, left=0.086, right=0.98,
                    hspace=0.05, wspace=0.07,)
if (ncols == 1) & (nrows == 1):
    axs = np.array([axs])
axs=axs.ravel()
[axi.set_axis_off() for axi in axs]   

ax=axs[0]
l2_prod.Rrs.sel(wl=[665,560,470],method='nearest').plot.imshow(rgb='wl',robust=True,ax=ax)

ax=axs[1]
l2_prod.brdfg_full.plot.imshow(vmin=0,cmap=plt.cm.Grays_r,robust=True,cbar_kwargs={'shrink':0.45},ax=ax)

ax=axs[2]
l2_prod.aot_ref.plot.imshow(cmap=plt.cm.Spectral_r,cbar_kwargs={'shrink':0.45},robust=True,ax=ax)

ax=axs[3]
l2_prod.aot_ref_smoothed.plot.imshow(cmap=plt.cm.Spectral_r,cbar_kwargs={'shrink':0.45},robust=True,ax=ax)

<matplotlib.image.AxesImage at 0x7cb56700fe30>

Plot and interact

import panel as pn

#hv.extension('bokeh')
pn.extension()
from grstbx import visual

v=visual.ViewSpectral(l2_prod.Rrs.isel(wl=range(0,63,5)) ,reproject=True)
print(v)
v.minmax=[0,0.1]
v.minmaxvalues=(0,0.04)

<grstbx.visual.ViewSpectral object at 0x7cb560c79cd0>

v.visu()

import geopandas as gpd
crs = img.rio.crs
poi_stream = v.poi_stream
geom = poi_stream.data
geom=gpd.points_from_xy(geom['x'],geom['y'])#,crs="EPSG3857")index=[0], c
geom =gpd.GeoDataFrame(crs=3857, geometry=geom).to_crs(crs)
Ndata = len(geom)
cmap = mpl.colors.LinearSegmentedColormap.from_list("",
                                                    ['navy', "blue", 'lightskyblue',
                                                     "grey",   'forestgreen','yellowgreen',
                                                     "khaki", "gold",
                                                     'orangered', "firebrick", 'purple'])
norm = mpl.colors.Normalize(vmin=0.0,vmax=Ndata-1)

wl_range=slice(400,900)
alpha=0.5
lw=1.2
ms=1.5

fig,axs = plt.subplots(1,2,figsize=(15,5))#,sharey=True)   
fig.subplots_adjust(hspace=0.05, wspace=0.25)
ax=axs[0]
indexes=np.dstack([geom.geometry.x,geom.geometry.y])[0]
for ii,index in enumerate(indexes):
    color=cmap(norm(ii))
   
    boa=l2_prod.Rrs.sel(x=index[0],y=index[1],method='nearest') +0.0025
    boa.sel(wl=wl_range).plot(label='pnt-'+str(ii),marker='o',ms=ms,lw=lw,c=color,alpha=alpha,ax=ax)
    boa.rolling(wl=3).mean().sel(wl=wl_range).plot(label='pnt-'+str(ii),marker='o',ms=ms,lw=lw,alpha=alpha,c=color,ax=axs[1])

    if False:
        toa=l1c_prod.Rtoa.sel(x=index[0],y=index[1],method='nearest')
        toa.plot(ms=ms,lw=lw,c=color,alpha=alpha,ax=ax)

for ax in  axs:
    ax.set_ylabel(r'$R_{rs}\ (sr^{-1})\ or\ R_{TOA}\ (-)$')
    ax.axhline(0,color='k')
    ax.set_xlabel(r'$Wavelength\ (nm)$')
    ax.set_title('')
    ax.minorticks_on() 
    #prod.Twv_raster.mean(['x','y']).plot(marker='o',ms=2.5,lw=0.5)
   

axs[0].legend(fontsize=10)

/tmp/ipykernel_104314/3862006790.py:44: UserWarning: No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.
  axs[0].legend(fontsize=10)

<matplotlib.legend.Legend at 0x7cb554549100>

import geopandas as gpd
crs = img.rio.crs
poi_stream = v.poi_stream
geom = poi_stream.data
geom=gpd.points_from_xy(geom['x'],geom['y'])#,crs="EPSG3857")index=[0], c
geom =gpd.GeoDataFrame(crs=3857, geometry=geom).to_crs(crs)
Ndata = len(geom)
cmap = mpl.colors.LinearSegmentedColormap.from_list("",
                                                    ['navy', "blue", 'lightskyblue',
                                                     "grey",   'forestgreen','yellowgreen',
                                                     "khaki", "gold",
                                                     'orangered', "firebrick", 'purple'])
norm = mpl.colors.Normalize(vmin=0.0,vmax=Ndata-1)


wl_range=slice(400,1100)
alpha=0.5
lw=1.2
ms=1.5

fig,axs = plt.subplots(1,2,figsize=(15,5))#,sharey=True)   
fig.subplots_adjust(hspace=0.05, wspace=0.25)
ax=axs[0]
indexes=np.dstack([geom.geometry.x,geom.geometry.y])[0]
for ii,index in enumerate(indexes):
    color=cmap(norm(ii))
   
   
    boa=l2_prod.Rrs.sel(x=index[0],y=index[1],method='nearest')
    
    boa.plot(label='pnt-'+str(ii),marker='o',ms=ms,lw=lw,c=color,alpha=alpha,ax=ax)
    
    boa.sel(wl=wl_range).plot(label='pnt-'+str(ii),marker='o',ms=ms,lw=lw,alpha=alpha,c=color,ax=axs[1])

    if False:
        toa=l1c_prod.Rtoa.sel(x=index[0],y=index[1],method='nearest')
        toa.plot(ms=ms,lw=lw,c=color,alpha=alpha,ax=ax)
        #toa.sel(wl=wl_range).plot(ms=ms,lw=lw,c=color,alpha=alpha,ax=axs[1])
#(0.04*prod.Tg_other).plot(marker='o',ms=2.5,lw=0.5,ax=axs[0])
#(0.04*prod.Tg_other).sel(wl=wl_range).plot(marker='o',color='k',ms=2.5,lw=0.5,ax=axs[1])
#(0.04*Twv).plot(marker='o',ms=2.5,lw=0.5,ax=axs[0])
#(0.04*Twv).sel(wl=wl_range).plot(marker='o',color='b',ms=2.5,lw=0.5,ax=axs[1])

for ax in  axs:
    ax.set_ylabel(r'$R_{rs}\ (sr^{-1})\ or\ R_{TOA}\ (-)$')
    ax.axhline(0,color='k')
    ax.set_xlabel(r'$Wavelength\ (nm)$')
    ax.set_title('')
    ax.minorticks_on() 
    #prod.Twv_raster.mean(['x','y']).plot(marker='o',ms=2.5,lw=0.5)
   

axs[0].legend(fontsize=10)

/tmp/ipykernel_104314/2529270420.py:53: UserWarning: No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.
  axs[0].legend(fontsize=10)

<matplotlib.legend.Legend at 0x7cb4f2194200>

Holoviews tools

from holoviews import streams
import holoviews as hv
import panel as pn
import param
hv.extension('bokeh')
from holoviews import opts
opts.defaults(
    opts.GridSpace(shared_xaxis=True, shared_yaxis=True),
    opts.Image(cmap='binary_r', width=800, height=700),
    opts.Labels(text_color='white', text_font_size='8pt', text_align='left', text_baseline='bottom'),
    opts.Path(color='white'),
    opts.Spread(width=900),
    opts.Overlay(show_legend=True))
# set the parameter for spectra extraction
hv.extension('bokeh')
pn.extension()



param = 'Rrs' #Rtoa'
#img = prod[['Rtoa','Ltoa']] 
raster = img[param]#L2grs #masked[param] 

#param = 'rho'
#raster = dc_l2c[param] 
cmap='Spectral_r'
#cmap='RdBu_r'
third_dim = 'wl'

wl= raster.wl.data
Nwl = len(wl)
ds = hv.Dataset(raster.persist())
im= ds.to(hv.Image, ['x', 'y'], dynamic=True).opts(cmap= cmap,colorbar=True,clim=(0,0.041)).hist(bin_range=(0,0.01)) 

polys = hv.Polygons([])
box_stream = hv.streams.BoxEdit(source=polys)
dmap, dmap_std=[],[]

def roi_curves(data,ds=ds):    
    if not data or not any(len(d) for d in data.values()):
        return hv.NdOverlay({0: hv.Curve([],'Wavelength (nm)', param)})

    curves,envelope = {},{}
    data = zip(data['x0'], data['x1'], data['y0'], data['y1'])
    for i, (x0, x1, y0, y1) in enumerate(data):
        selection = ds.select(x=(x0, x1), y=(y0, y1))
        mean = selection.aggregate(third_dim, np.mean).data
        std = selection.aggregate(third_dim, np.std).data
        wl = mean.wl

        curves[i]= hv.Curve((wl,mean[param]),'Wavelength (nm)', param) 

    return hv.NdOverlay(curves)


# a bit dirty to have two similar function, but holoviews does not like mixing Curve and Spread for the same stream
def roi_spreads(data,ds=ds):    
    if not data or not any(len(d) for d in data.values()):
        return hv.NdOverlay({0: hv.Curve([],'Wavelength (nm)', param)})

    curves,envelope = {},{}
    data = zip(data['x0'], data['x1'], data['y0'], data['y1'])
    for i, (x0, x1, y0, y1) in enumerate(data):
        selection = ds.select(x=(x0, x1), y=(y0, y1))
        mean = selection.aggregate(third_dim, np.mean).data
        std = selection.aggregate(third_dim, np.std).data
        wl = mean.wl

        curves[i]=  hv.Spread((wl,mean[param],std[param]))#.opts(fill_alpha=0.3)

    return hv.NdOverlay(curves)

mean=hv.DynamicMap(roi_curves,streams=[box_stream])
std =hv.DynamicMap(roi_spreads, streams=[box_stream])    
hlines = hv.HoloMap({wl[i]: hv.VLine(wl[i]) for i in range(Nwl)},third_dim )


hv.output(widget_location='top_left')

# visualize and play
graphs = ((mean* std *hlines).relabel(param))
layout = (im * polys +graphs    ).opts(
    opts.Curve(width=750,height=500, framewise=True,xlim=(400,1100)), 
    opts.Polygons(fill_alpha=0.2, color='green',line_color='black'), 
    opts.VLine(color='black')).cols(2)
layout