Development example for GRS step-by-step application to Landsat images#

'''
Main program
'''

import os
import glob

import matplotlib as mpl
import matplotlib.pyplot as plt

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

import logging

import panel as pn

import cartopy.crs as ccrs
import geopandas as gpd

import GRSdriver
import grs
from grs import Product, AuxData, acutils, CamsProduct, L2aProduct, Rasterization
import grstbx
from grstbx import visual

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

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

-GRSdriver: 1.0.4
-grs: 2.1.6
-grstbx: 2.0.2

Indicate the path where you put the look-up table file#

lut_file = '/data/vrtc/xlut/toa_lut_opac_wind_light.nc'
lut_file = '/data/vrtc/xlut/toa_lut_opac_wind_light_v2.nc'
trans_lut_file = '/data/vrtc/xlut/transmittance_lut_opac_wind_light_v2.nc'

Set the path of your Sentinel-2 image and corresponding CAMS file#

file ='/data/satellite/landsat/LC09_L1TP_196030_20221120_20230321_02_T1'

file_nc = file+'.nc'

cams_file = '/data/cams/world/cams_forecast_2022-11.nc'

l1c = GRSdriver.LandsatDriver(file)
l1c.INFO

bandId	0	1	2	3	4	5	6	7	8
ESA	B01	B02	B03	B08	B04	B05	B09	B06	B07
EOREADER	CA	BLUE	GREEN	PAN	RED	NIR	SWIR_CIRRUS	SWIR_1	SWIR_2
Wavelength (nm)	443	490	560	590	665	865	1370	1610	2190
Band width (nm)	25	60	60	173	33	28	21	95	287
Resolution (m)	30	30	30	15	30	30	30	30	30

l1c.load_mask()

l1c.load_product()

/home/harmel/anaconda3/envs/grstbx/lib/python3.9/site-packages/xarray/core/indexes.py:659: RuntimeWarning: '<' not supported between instances of 'SpectralBandNames' and 'SpectralBandNames', sort order is undefined for incomparable objects.
  new_pd_index = pd_indexes[0].append(pd_indexes[1:])

Load the image and convert it into netcdf for further fatest loading (not necessary)#

#l1c.prod.to_netcdf(file_nc)
#prod = Product(xr.open_dataset(file_nc))

Visualize, interact and subset your region of interest#

Yuo can use the “polygon” tool to subset the image and proceed with further processing

v=visual.view_spectral(prod.raster.bands,reproject=True)

v.title="## Landsat-9 L1C"
v.minmaxvalues=(0,0.3)
v.minmax=[0,0.5]
v.visu()

You can save the area of interest into a json file for reuse (uncomment the next two cells)

#geom_ = v.get_geom(v.aoi_stream,crs=prod.raster.rio.crs)

#geom_.to_file("roi_example.json", driver="GeoJSON")

geom_ = gpd.read_file("roi_example.json", driver="GeoJSON")

raster_clipped = xr.Dataset()
prod.raster.bands.rio.clip(geom_.geometry.values)
for param in prod.raster.keys():
    da = prod.raster[param]
    if 'x' in da.dims and 'y' in da.dims:
        raster_clipped[param]=da.rio.clip(geom_.geometry.values)
    else:
        raster_clipped[param]=da
raster_clipped.attrs = prod.raster.attrs

str_epsg = str(raster_clipped.rio.crs)
zone = str_epsg[-2:]
is_south = str_epsg[2] == 7
proj = ccrs.UTM(zone, is_south)

plt.figure(figsize=(15,15))
raster_clipped.bands.sel(wl=[665,560,490]).plot.imshow(rgb='wl', robust=True,subplot_kws=dict(projection=proj))

<matplotlib.image.AxesImage at 0x7fd7bbd2e550>

../_images/80aa21f73330b0f4667617b38beaa22924773858014e3c546a0cc5ff0cbeaf0c.png

Set aerosol model from CAMS data#

Get spectral aerosol optical thickness from CAMS raster

cams_aot_mean = cams.cams_aod.mean(['x','y'])
cams_aot_ref = cams.cams_aod.interp(wl=550,method='quadratic').compute()
cams_aot_ref_mean =  cams_aot_ref.mean(['x','y'])

Check proximity with tabulated models (LUT) and select representative aerosol model

fig, axs = plt.subplots(1, 2, figsize=(18, 4.5))
aero_lut.aot.sel(model=models,aot_ref=1).plot(ax=axs[0],hue='model')#,add_legend=False)
(cams_aot_mean/cams_aot_ref_mean).plot(ax=axs[0],color='black')
aero_lut.aot.isel(model=[4,2]).mean('model').sel(aot_ref=1).plot(ax=axs[0],color='grey')
lut_aod=aero_lut.aot.sel(model=models,aot_ref=1).interp(wl=cams.cams_aod.wl)
rank = np.abs((cams_aot_mean/cams_aot_ref_mean)-lut_aod).sum('wl') 
axs[1].bar(x=rank.model, height=rank.values)
plt.xticks(rotation=30, ha='right')
plt.show()

idx = np.abs((cams_aot_mean/cams_aot_ref_mean)-lut_aod).sum('wl').argmin()
opac_model = aero_lut.sel(model=models).model.values[idx]
print(opac_model)

../_images/95262bbfa540ac1c1d4686f04f270b8d729d69303acfb0b195bde6b75165a242.png

COAV_rh70

absorbing gases correction#

gas_trans = acutils.GaseousTransmittance(prod, cams)
Tg_raster = gas_trans.get_gaseous_transmittance()
#Tg_raster

logging.info('correct for gaseous absorption')
for wl in prod.raster.wl.values:
    Tg_ = Tg_raster.sel(wl=wl).interp(x=prod.raster.x, y=prod.raster.y)
    prod.raster['bands'].loc[wl] = prod.raster.bands.sel(wl=wl) / Tg_
    del Tg_
prod.raster.bands.attrs['gas_absorption_correction'] = True

INFO:root:correct for gaseous absorption

import gc
gc.collect()

You can check the gaseous transmittance for each spectral band

plt.figure()
Tg_raster.mean('x').mean('y').plot()
fig,axs = plt.subplots(3,4,figsize=(15,9),sharey=True,sharex=True)
axs=axs.ravel()
[axi.set_axis_off() for axi in axs]
for iwl in range(len(prod.wl_process)):
    #axs[iwl].set_axis_on()
    Tg_raster.isel(wl=iwl).plot(ax=axs[iwl],cmap=plt.cm.Spectral_r)
    axs[iwl].set_title( str(Tg_raster.isel(wl=iwl).wl.values))

../_images/868e2ac3aab30a90af644b1df96b55ed84e35b74d7b99011753ba008e06f0820.png

../_images/ce70761bd50ea9aa2851f470d1b2af57ae0547de3cf89494b5e3234b64a39bfc.png

L2A to L2B#

Mask bad quality pixels#

Rrs_blue_avg=0.05
Rrs = l2_prod.Rrs
mask=xr.where((Rrs.sel(wl=490,method='nearest')<0)|(Rrs.sel(wl=565,method='nearest')<0),1,0)
mask=mask.where((Rrs.sel(wl=443,method='nearest')<2*Rrs_blue_avg),2)
mask=mask.where((Rrs.sel(wl=443,method='nearest')<3*Rrs_blue_avg),3)
#mask=mask.where((l2a.l2_prod.wv_band)<0.1,4)


plt.figure(figsize=(15,15))

mask.plot.imshow(cmap=plt.cm.Spectral_r,robust=True,subplot_kws=dict(projection=proj),cbar_kwargs={'shrink':0.35})

<matplotlib.image.AxesImage at 0x7fd7ba917c40>

../_images/3c633407b960b76a002b2471502fdeb4987a6eb78990af51762ea79f26cfcd50.png

Apply mask

Rrs=Rrs.where(mask==0)

plt.figure(figsize=(15,15))
#(dem.shaded+1).plot.imshow(robust=True,subplot_kws=dict(projection=proj),cmap=plt.cm.Greys_r,cbar_kwargs={'shrink':0.4})
Rrs.sel(wl=[665,560,490]).plot.imshow(rgb='wl', robust=True,subplot_kws=dict(projection=proj))

<matplotlib.image.AxesImage at 0x7fd7ba1185b0>

../_images/e661727befe46f844e550236a07792180f9359d11ca535123a45fb27647422f9.png

Check Rrs spectra#

v=visual.view_spectral(l2_prod.Rrs,reproject=True)

v.minmax=[0,0.1]
v.visu()

geom_ = v.get_geom(v.aoi_stream,crs=prod.raster.rio.crs)

raster_clipped = xr.Dataset()
l2_prod.rio.clip(geom_.geometry.values)
for param in l2_prod.keys():
    da = l2_prod[param]
    if 'x' in da.dims and 'y' in da.dims:
        raster_clipped[param]=da.rio.clip(geom_.geometry.values)
    else:
        raster_clipped[param]=da
raster_clipped.attrs = prod.raster.attrs

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
Cell In[55], line 1
----> 1 geom_ = v.get_geom(v.aoi_stream,crs=prod.raster.rio.crs)
      3 raster_clipped = xr.Dataset()
      4 l2_prod.rio.clip(geom_.geometry.values)

File ~/Dropbox/Dropbox/work/git/satellite_app/grstbx/grstbx/visual.py:237, in utils.get_geom(aoi_stream, crs)
    234 @staticmethod
    235 def get_geom(aoi_stream, crs=4326):
    236     geom = aoi_stream.data
--> 237     ys, xs = geom['ys'][-1], geom['xs'][-1]
    238     polygon_geom = Polygon(zip(xs, ys))
    239     polygon = gpd.GeoDataFrame(index=[0], crs=3857, geometry=[polygon_geom])

IndexError: list index out of range

stacked = raster_clipped.Rrs.sel(wl=slice(400,1000)).stack(gridcell=["y", "x"]).dropna('gridcell',thresh=0)

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[56], line 1
----> 1 stacked = raster_clipped.Rrs.sel(wl=slice(400,1000)).stack(gridcell=["y", "x"]).dropna('gridcell',thresh=0)

File ~/anaconda3/envs/grstbx/lib/python3.9/site-packages/xarray/core/common.py:277, in AttrAccessMixin.__getattr__(self, name)
    275         with suppress(KeyError):
    276             return source[name]
--> 277 raise AttributeError(
    278     f"{type(self).__name__!r} object has no attribute {name!r}"
    279 )

AttributeError: 'Dataset' object has no attribute 'Rrs'

group_coord ='wl'
stat_coord='gridcell'
stats = xr.Dataset({'median':stacked.groupby(group_coord).median(stat_coord)})
stats['q25'] = stacked.groupby(group_coord).quantile(0.25,dim=stat_coord)
stats['q75'] = stacked.groupby(group_coord).quantile(0.75,dim=stat_coord)
stats['min'] = stacked.groupby(group_coord).min(stat_coord)
stats['max'] = stacked.groupby(group_coord).max(stat_coord)
stats['mean'] = stacked.groupby(group_coord).mean(stat_coord)
stats['std'] = stacked.groupby(group_coord).std(stat_coord)
stats['pix_num'] = stacked.count(stat_coord)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[57], line 3
      1 group_coord ='wl'
      2 stat_coord='gridcell'
----> 3 stats = xr.Dataset({'median':stacked.groupby(group_coord).median(stat_coord)})
      4 stats['q25'] = stacked.groupby(group_coord).quantile(0.25,dim=stat_coord)
      5 stats['q75'] = stacked.groupby(group_coord).quantile(0.75,dim=stat_coord)

NameError: name 'stacked' is not defined

from mpl_toolkits.axes_grid1.inset_locator import inset_axes

bands=[3,2,1]
fig, axs = plt.subplots(1,1, figsize=(10, 6))#,sharey=True



date = stats.time.dt.date.values
axins = inset_axes(axs, width="45%", height="60%",
               bbox_to_anchor=(.55, .4, 0.9, 0.9),
               bbox_transform=axs.transAxes, loc=3)


raster_clipped.Rrs.isel(wl=bands).plot.imshow(robust=True,ax=axins)
axins.set_title('')
axins.set_axis_off()
axs.axhline(y=0,color='k',lw=1)
axs.plot(stats.wl,stats['median'],marker='o',c='k')
axs.plot(stats.wl,stats['mean'],c='red',ls='--')
axs.fill_between(stats.wl, stats['q25'], stats['q75'],alpha=0.3,color='grey')
axs.set_title(date)
plt.show()

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[58], line 8
      3 bands=[3,2,1]
      4 fig, axs = plt.subplots(1,1, figsize=(10, 6))#,sharey=True
----> 8 date = stats.time.dt.date.values
      9 axins = inset_axes(axs, width="45%", height="60%",
     10                bbox_to_anchor=(.55, .4, 0.9, 0.9),
     11                bbox_transform=axs.transAxes, loc=3)
     14 raster_clipped.Rrs.isel(wl=bands).plot.imshow(robust=True,ax=axins)

NameError: name 'stats' is not defined

../_images/7846670f48a0ac9879147bcc49bcbd54d8239704ca32748c270ea5b8eda22b86.png

Check water quality parameters (e.g., Chl-a concentration from diverse “algorithms”)¶#

Total suspended particulate matter (SPM) from Nechad et al., 2010, 2016 formulation#

spm in mg/L#

a = [610.94*np.pi, 0.2324/np.pi]
Rrs_ = Rrs.isel(wl=3)
spm = a[0] * Rrs_ / (1 - ( Rrs_/ a[1]))
spm.name='SPM_N2016'
spm.attrs['units']='mg/L'
spm.attrs['description']='Concentration of suspended particulate matter from band 665 nm'

a = [428.277*np.pi, 0.3051/np.pi]
Rrs_ = Rrs.isel(wl=3)
spm2 = a[0] * Rrs_ / (1 - ( Rrs_/ a[1]))
spm2.name='SPM_N2016_han1'
spm2.attrs['units']='mg/L'
spm2.attrs['description']='Concentration of suspended particulate matter from band 665 nm'


a = [1444.853*np.pi, 0.3539/np.pi]
Rrs_ = Rrs.isel(wl=3)
spm3 = a[0] * Rrs_ / (1 - ( Rrs_/ a[1]))
spm3.name='SPM_N2016_han2'
spm3.attrs['units']='mg/L'
spm3.attrs['description']='Concentration of suspended particulate matter from band 665 nm'

plt.figure(figsize=(25,18))
colors = ['mediumblue','cadetblue','teal','gold','orange','chocolate','firebrick']
cmap = mpl.colors.LinearSegmentedColormap.from_list('spm',colors)
fig = plt.figure(figsize=(25, 15))
ax = plt.subplot(1, 3, 1, projection=proj)
spm.plot.imshow(cmap=cmap,ax=ax,robust=True,cbar_kwargs={'shrink':0.25},vmin=0)
ax.set_title('Nechad et al., 2016')
ax = plt.subplot(1, 3,2, projection=proj)
spm2.plot.imshow(cmap=cmap,ax=ax,robust=True,cbar_kwargs={'shrink':0.25},vmin=0)
ax.set_title(' Nechad_Han SAA')
ax = plt.subplot(1, 3,3, projection=proj)
spm3.plot.imshow(cmap=cmap,ax=ax,robust=True,cbar_kwargs={'shrink':0.25},vmin=0)
ax.set_title('Nechad_HAN SAA-H')

Text(0.5, 1.0, 'Nechad_HAN SAA-H')

<Figure size 2500x1800 with 0 Axes>

../_images/3bc044cf3d44e23d5731975dc65be733125c11e94306719a460d5a71c97dea67.png

Check blue over green ratio for Chl retrieval with OC2 from NASA#

\(log_{10}(chlor\_a) = a_0 + \sum\limits_{i=1}^4 a_i \left(log_{10}\left(\frac{R_{rs}(\lambda_{blue})}{R_{rs}(\lambda_{green})}\right)\right)^i\)

# NASA OC2 fro MODIS; bands 488, 547 nm
a_modis = [0.2500,-2.4752,1.4061,-2.8233,0.5405]
# NASA OC2 for OCTS; bands 490, 565 nm
a_octs = [0.2236,-1.8296,1.9094,-2.9481,-0.1718]
# Pelevin et al, 2023, Issyk-Kul
a = [0.1977,-1.8117,1.9742,-2.5635,-0.7218]
# mona et al.
a_mona = [0.484,-2.109,2.880,-0.690,-1.040]

ratio = np.log10(Rrs.isel(wl=1)/Rrs.isel(wl=2))
chl={}
for name,a in [('nasa',a_octs),('mona',a_mona)]:
    logchl=0
    for i in range(len(a)):
        logchl+=a[i]*ratio**i
    _chl = 10**(logchl)
    _chl.name='Chla_'+name
    _chl.attrs['units']='mg.m-3'
    _chl.attrs['description']= 'Chl-a concentration from NASA OC2 with OCTS parameterization, bands 490, 565 nm',
    _chl = _chl.where((_chl >= 0) & (_chl <= 2000))
    chl[name]=_chl

import colorcet as cc
ncols=2
nrows=1
vmax=20
fig,axs = plt.subplots(nrows,ncols,figsize=(ncols*5.1,5*nrows+5),sharey=True,sharex=True,subplot_kw={'projection': proj})   
fig.subplots_adjust(bottom=0.08, top=0.9, left=0.086, right=0.98,
                    hspace=0.15, wspace=0.12,)
for ii, (name,_chl) in enumerate(chl.items()):
    _chl.plot.imshow(cmap=cc.cm.CET_D13,robust=True,ax=axs[ii],cbar_kwargs={'shrink':0.35},vmax=vmax)

../_images/805c98c3d8742c1980f2d1592082edb82ed143024d8258742c46e348236b8902.png

CDOM retrieval based on Brezonik et al, 2015#

a = [1.872,-0.83]
acdom = np.exp(a[0] + a[1] * np.log(Rrs.isel(wl=1)/Rrs.isel(wl=5)))
acdom.name='a_cdom_b2015'
acdom.attrs['units']='m-1'
acdom.attrs['description']='CDOM absorption at 440 nm-1'
acdom= acdom.where((acdom >= 0) & (acdom <= 60))

plt.figure(figsize=(15,15))

acdom.plot.imshow(cmap=cc.cm.CET_CBD1,robust=True,subplot_kws=dict(projection=proj),cbar_kwargs={'shrink':0.35},vmax=4)

<matplotlib.image.AxesImage at 0x7fd788389190>

../_images/57c6b8a49f969320565dbed539572e023932dabe4c5985fdd42de9e1eaab9173.png

def cPOC_2(Rrs,p=[2.873,0.945,0.025]):
    ratio1=Rrs.sel(wl=665,method='nearest') / Rrs.sel(wl=490,method='nearest') 
    ratio2=Rrs.sel(wl=665,method='nearest') / Rrs.sel(wl=565,method='nearest') 
    ratio = np.log10(xr.concat([ratio1.assign_coords({'num':1}),ratio1.assign_coords({'num':2})],dim='num').max('num'))
    
    Xpoc = p[0]+p[1]*ratio+p[2]*ratio**2
    return 10**Xpoc

poc = cPOC_2(Rrs)
poc.name = 'cPOC_2'

plt.figure(figsize=(15,15))

poc.plot.imshow(cmap=cc.cm.blues,robust=True,subplot_kws=dict(projection=proj),cbar_kwargs={'shrink':0.35})

<matplotlib.image.AxesImage at 0x7fd788327220>

../_images/0e804d960822cd877f0a0bef7f985618cc5142731365fa3d44f2f6a685c68841.png

Check Optical Water Types (OWT)#

owt_file = ‘/DATA/projet/vrac/owt/Spyrakos_et_al_2018_OWT_inland_mean_standardised.csv’ owt =pd.read_csv(owt_file,index_col=0).stack().to_xarray() owt = owt.rename({‘level_1’:‘wl’}) owt[‘wl’]=owt.wl.astype(float)

import matplotlib.patches as mpatches

owt_info={ 1:dict(color=‘olivedrab’,label=‘OWT1: Hypereutrophic waters’), 2:dict(color=‘black’,label=‘OWT2: Common case waters’), 3:dict(color=‘cadetblue’,label=‘OWT3: Clear waters’), 4:dict(color=‘tan’,label=‘OWT4: Turbid waters with organic content’), 5:dict(color=‘chocolate’,label=‘OWT5: Sediment-laden waters’), 6:dict(color=‘teal’,label=‘OWT6: Balanced optical effects at shorter wavelengths’), 7:dict(color=‘blueviolet’,label=‘OWT7: Highly productive cyanobacteria-dominated waters’), 8:dict(color=‘plum’,label=‘OWT8: Productive with cyanobacteria waters’), 9:dict(color=‘red’,label=‘OWT9: OWT2 with higher \(R_{rs}\) at shorter wavelengths’),#‘slategrey’ 10:dict(color=‘orange’,label=‘OWT10: CDOM-rich waters’), 11:dict(color=‘gold’,label=‘OWT11: CDOM-rich with cyanobacteria waters’), 12:dict(color=‘firebrick’,label=‘OWT12: Turbid waters with cyanobacteria’), 13:dict(color=‘mediumblue’,label=‘OWT13: Very clear blue waters’), } colors = [‘olivedrab’,‘black’,‘cadetblue’,‘tan’,‘chocolate’,‘teal’,‘blueviolet’,‘plum’,‘red’,‘orange’,‘gold’,‘firebrick’,‘mediumblue’] cmap_owt = mpl.colors.ListedColormap(colors)

patch = [] for key,info in owt_info.items(): patch.append(mpatches.Patch(color=info[‘color’], label=info[‘label’]))

fig, ax = plt.subplots(nrows=1,ncols=1, sharex=True,figsize=(9, 6)) ax.minorticks_on() for iowt,group in owt.groupby(‘owt’):

group.plot(color=owt_info[iowt]['color'],lw=3)

ax.set_title(‘’ ) ax.set_ylabel(‘\(Standardized\ R_{rs}\ (nm^{-1})\)’,fontsize=20) ax.set_xlabel(‘\(Wavelength\ (nm)\)’,fontsize=20)
plt.legend(handles=patch,fontsize=13,bbox_to_anchor=(1, .5, 0.5, 0.5))

def SAM(R1,R2): denum=(R1*R2).sum(‘wl’) denom = (R1**2).sum(‘wl’)0.5 * (R22).sum(‘wl’)**0.5 return np.arccos(denum/denom)

def SCS(R1,R2): R1_avg = R1.mean(‘wl’) R2_avg = R2.mean(‘wl’) R1_std = R1.std(‘wl’) R2_std = R2.std(‘wl’) Nwl = len(R1.wl) return 1/(Nwl) * ((R1-R1_avg) * (R2-R2_avg)).sum(‘wl’) / (R1_std*R2_std)

Rrs_sat = Rrs.sel(wl=slice(350,800))#.dropna(‘wl’)

Rrs_owt = owt.interp(wl=Rrs_sat.wl)

owt_sam = SAM(Rrs_sat,Rrs_owt) owt_scs =SCS(Rrs_sat,Rrs_owt)

owt_delta = owt_scs + (1-2*owt_sam/np.pi)/2

fig = plt.figure(figsize=(25, 15)) ax = plt.subplot(1, 2, 1, projection=proj)

cmap = plt.get_cmap(‘tab20c’,13) (owt_delta.fillna(-1).argmax(‘owt’)+1).where(Rrs_sat.isel(wl=1)>0).plot.imshow(vmin=0.5,vmax=13.5,cmap=cmap_owt,cbar_kwargs={ ‘ticks’:range(1,14),‘shrink’: 0.4},ax=ax) ax = plt.subplot(1, 2, 2, projection=proj) cmap = plt.get_cmap(‘RdBu’)#,13) (owt_delta.max(‘owt’)).where(Rrs_sat.isel(wl=1)>0).plot.imshow(robust=True,cmap=cmap,ax=ax, cbar_kwargs={ ‘shrink’: 0.4})

#l2b.to_netcdf('/DATA/vrac/test_l2b.nc')

Development example for GRS step-by-step application to Landsat images

Contents