L2A GRS image

Contents

L2A GRS image#

Simple example to open and visualize L2A GRS image

import os, glob

import matplotlib
import matplotlib.pyplot as plt

import numpy as np
import xarray as xr
import rioxarray as rio
import cartopy.crs as ccrs

opj = os.path.join

Geoinformation and plotting#

You can get the projection information and then plot the RGB image as follows:

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

plt.figure(figsize=(10,10))
l2a.Rrs.sel(wl=[665,560,490]).where(l2a.mask==0).plot.imshow(rgb='wl', robust=True,subplot_kws=dict(projection=proj))

<matplotlib.image.AxesImage at 0x7f5ae31521f0>

/home/harmel/anaconda3/envs/grstbx/lib/python3.9/site-packages/matplotlib/cm.py:496: RuntimeWarning: invalid value encountered in cast
  xx = (xx * 255).astype(np.uint8)

../_images/28135794bf30ab665ec03e4a01e2da9dec4c7c57c05f1623385006831fe39ac7.png

GRStbx#

Otherwise, I might be interested in using GRSTBX (Tristanovsk/grstbx) that provides the tools for subsettint and visualizing temporal datacubes

import grstbx
from grstbx import visual
import geopandas as gpd
import holoviews as hv
hv.extension('bokeh')

v=visual.ViewSpectral(l2a.Rrs,reproject=True)

The next cell will render the visualization tool. For subsetting, you can choose the polygon drawer tool on the right side and then press long to generate the polygon. Short clicks will add a new segment and long click end the polygon.

v.visu()

The next cell was commented for rendendering in Readthedocs (please uncomment for proper use of the notebook)

#geom_ = v.get_geom(v.aoi_stream,crs=l2a.rio.crs)
#geom_.to_file("roi_grs_visu_simple.json", driver="GeoJSON")

geom_ = gpd.read_file("roi_grs_visu_simple.json", driver="GeoJSON")
Rrs_clipped = l2a.Rrs.sel(wl=slice(400,1000)).where(l2a.mask==0).rio.clip(geom_.geometry.values)

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

<matplotlib.image.AxesImage at 0x7f5c0a526070>

/home/harmel/anaconda3/envs/grstbx/lib/python3.9/site-packages/matplotlib/cm.py:496: RuntimeWarning: invalid value encountered in cast
  xx = (xx * 255).astype(np.uint8)

../_images/3982c2ba77cf50a6fb1b2eaaab5951b48f4c132ce99020a376c8f159bbf0de4d.png

Rrs_clipped.plot.imshow(col='wl',col_wrap=3, robust=True, cmap=plt.cm.Spectral_r,
                                              subplot_kws=dict(projection=proj))

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

../_images/a584811195d07e405c163360704261801b00e89451ca3c6df698a9a0037126ab.png