CAMS data#

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'
cams_file = '/data/cams/world/cams_forecast_2022-11.nc'

l1c = GRSdriver.LandsatDriver(image,band_idx=[1,2,4,6,9,10])

l1c.INFO

bandId	1	2	4	6	9	10
ESA	B02	B03	B04	B09	B08	B09
EOREADER	BLUE	GREEN	RED	SWIR_CIRRUS	TIR_1	TIR_2
Wavelength (nm)	490	560	665	1370	11000	12000
Band width (nm)	60	60	33	21	590	1010
Resolution (m)	30	30	30	30	100	100

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:])

Open CAMS data with lazy loading

Get data for the date of the image

date = l1c.prod.time.dt.strftime('%Y-%m-%d').values
date

array('2022-11-20', dtype=object)

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

(cams.t2m-273.15).plot.imshow(col='time',col_wrap=2,cmap=plt.cm.Spectral_r,aspect=1.75,robust=True,subplot_kws=dict(projection=proj),cbar_kwargs={'shrink':0.35}) #,vmin=-5,vmax=20)

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

../_images/391c2e00213c57cf31d35768c3e014325fec796cdd5a5c4ae5300bc2d94112ad.png

cams.tcwv.plot.imshow(col='time',col_wrap=4,cmap=plt.cm.Spectral_r,aspect=1.75,robust=True,subplot_kws=dict(projection=proj),cbar_kwargs={'shrink':0.35})

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

../_images/b90fe485d3a372a4877cced923543ffee46beb45b1e312d6c0ebc290054cec44.png

Now you can subset to the landsat image extent

cams.rio.bounds()

(-180.20000000339462, -90.2, 179.80000610691025, 90.2)

if your subset is small, it might be useful to interpolate within the native raster

cams=cams.rio.reproject(l1c.prod.rio.crs,nodata=np.nan)

cams.tcwv.plot.imshow(col='time',col_wrap=4,cmap=plt.cm.Spectral_r,vmin=-5,vmax=20)

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

../_images/a868b4c9cb6b01a688ebd24d6e7964249c2bd7f8c106b12d3ddf6587c3ae5741.png

Final step is to clip the cams rasters to that of your satellite image

cams.tcwv.plot.imshow(col='time',col_wrap=4,cmap=plt.cm.Spectral_r)

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

../_images/9cdd8724ea70344d8271023b786885c8a925cdc666bfc38fbfa2f93fd39fb0f7.png

Plotting examples:#

Select the closest date

cams_ = cams.sel(time=l1c.prod.time,method='nearest')

wind = np.sqrt(cams_['v10']**2+cams_['u10']**2)

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)

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

l1c.prod.bands.sel(wl=[665,560,490])[:,::5,::5].plot.imshow(rgb='wl', robust=True,subplot_kws=dict(projection=proj))
wind.plot.imshow( cmap=plt.cm.Spectral_r,robust=True,alpha=0.53,subplot_kws=dict(projection=proj),cbar_kwargs={'shrink':0.35})

<matplotlib.image.AxesImage at 0x7f4f58224eb0>

../_images/f89b3d4ab21867598c658e5dec2aa30ee31aec48c0da24dc79faca224b27f7c9.png

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

l1c.prod.bands.isel(wl=-1)[::5,::5].plot.imshow(cmap=plt.cm.Spectral_r,robust=True,subplot_kws=dict(projection=proj),cbar_kwargs={'shrink':0.35})
wind.plot.imshow( cmap=plt.cm.gray,robust=True,alpha=0.53,subplot_kws=dict(projection=proj),cbar_kwargs={'shrink':0.35})

<matplotlib.image.AxesImage at 0x7f4f5800f2b0>

../_images/deef9498c365be20ed017430f5bf30440539d716e5aef1f17507db8e4bf989b5.png

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

l1c.prod.bands.isel(wl=-1)[::5,::5].plot.imshow(cmap=plt.cm.Spectral_r,robust=True,subplot_kws=dict(projection=proj),cbar_kwargs={'shrink':0.35})
(cams_.t2m-273.15).plot.imshow( cmap=plt.cm.gray,robust=True,alpha=0.53,subplot_kws=dict(projection=proj),cbar_kwargs={'shrink':0.35})

<matplotlib.image.AxesImage at 0x7f4f4ff07d30>

../_images/062be897721119b7fdd4b1e0a2b685464807139a17c8a27d743e26a94daa2946.png

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

l1c.prod.bands.isel(wl=-1)[::5,::5].plot.imshow(cmap=plt.cm.Spectral_r,robust=True,subplot_kws=dict(projection=proj),cbar_kwargs={'shrink':0.35})
cams_.tcwv.plot.imshow( cmap=plt.cm.gray,robust=True,alpha=0.753,subplot_kws=dict(projection=proj),cbar_kwargs={'shrink':0.35})

<matplotlib.image.AxesImage at 0x7f4f4fe5dee0>

../_images/1300e9344317b6d33eaab3fd6c1831587207a3fce759956d52f18cdf1ebc26ab.png

h = visual.ViewSpectral(l1c.prod.bands[:,::5,::5],reproject=True)#.assign_coords(time=l1c.datetime).expand_dims('time')

h.visu()

CAMS data

Contents

CAMS data#

Plotting examples:#