<jupyter_start><jupyter_text>Accessing Sentinel-5P data on AzureThe [Sentinel-5P](https://sentinel.esa.int/web/sentinel/missions/sentinel-5p) mission provides daily global atmospheric measurements at a resolution of 3.5km x 7km on most bands.  This dataset represents the global archive of Sentinel-5P [Level 2](http://www.tropomi.eu/data-products/level-2-products) products, from 2018 to the present, in NetCDF format.  This notebook demonstrates basic access to Sentinel-5 data on Azure.  Because Sentinel-1 data are in preview, the user needs to provide storage credentials.This dataset is stored in the West Europe Azure region, so this notebook will run most efficiently on Azure compute located in the same region.  If you are using this data for environmental science applications, consider applying for an [AI for Earth grant](http://aka.ms/ai4egrants) to support your compute requirements.This dataset is documented at [aka.ms/ai4edata-sentinel-5p](http://aka.ms/ai4edata-sentinel-5p).Sentinel-5P data on Azure are maintained by [Sinergise](https://sinergise.com/).  Environment setup<jupyter_code>import os
import fsspec
import xarray as xr
import numpy as np
from azure.storage.blob import ContainerClient

# Not used directly, but needs to be installed to read NetCDF files with xarray
import h5netcdf<jupyter_output><empty_output><jupyter_text>Auth files<jupyter_code># A plain-text file with a SAS token (starting with "?sv") on the first line
sas_file = os.path.expanduser('~/tokens/sentinel-5p_sas.txt')<jupyter_output><empty_output><jupyter_text>Constants<jupyter_code># Let's look at ozone concentration from mid-day on Jan 1, 2021
product = 'L2__O3____'
date = '2021/01/01'<jupyter_output><empty_output><jupyter_text>Azure storage constants<jupyter_code>lines = []
with open(sas_file,'r') as f:
    lines = f.readlines()
assert len(lines) >= 1
sas_token = lines[0].strip()
        
storage_account_name = 'sentinel5euwest'
container_name = 'sentinel-5p'
storage_account_url = 'https://' + storage_account_name + '.blob.core.windows.net/'

container_client = ContainerClient(account_url=storage_account_url, 
                                                 container_name=container_name,
                                                 credential=sas_token)<jupyter_output><empty_output><jupyter_text>List products matching our product/date<jupyter_code>prefix = '/'.join(['TROPOMI',product,date])
print('Searching for prefix {}'.format(prefix))
generator = container_client.list_blobs(name_starts_with=prefix)
scene_paths = [blob.name for blob in generator]
print('\nFound {} matching scenes:\n'.format(len(scene_paths)))
for s in scene_paths:
    print(s.split('/')[-1])<jupyter_output>Searching for prefix TROPOMI/L2__O3____/2021/01/01

Found 14 matching scenes:

S5P_OFFL_L2__O3_____20210101T013036_20210101T031206_16680_01_020104_20210102T182728.nc
S5P_OFFL_L2__O3_____20210101T031206_20210101T045336_16681_01_020104_20210102T203532.nc
S5P_OFFL_L2__O3_____20210101T045336_20210101T063506_16682_01_020104_20210102T224247.nc
S5P_OFFL_L2__O3_____20210101T063506_20210101T081636_16683_01_020104_20210103T000719.nc
S5P_OFFL_L2__O3_____20210101T081636_20210101T095806_16684_01_020104_20210103T013023.nc
S5P_OFFL_L2__O3_____20210101T095806_20210101T113936_16685_01_020104_20210103T031528.nc
S5P_OFFL_L2__O3_____20210101T113936_20210101T132106_16686_01_020104_20210103T052522.nc
S5P_OFFL_L2__O3_____20210101T132106_20210101T150236_16687_01_020104_20210103T065107.nc
S5P_OFFL_L2__O3_____20210101T150236_20210101T164406_16688_01_020104_20210103T081248.nc
S5P_OFFL_L2__O3_____20210101T164406_20210101T182536_16689_01_020104_20210103T101745.nc
S5P_OFFL_L2__O3_____20210101T182536_20210101T200706[...]<jupyter_text>Print metadata for one sceneChoose an OFFL scene for this product/date (NRTI products have smaller domains).  In practice the scene we choose impactsthe longitude we're plotting, since S5P has a daily orbit.  We'll choose a mid-day scene.<jupyter_code>offl_scenes = [s for s in scene_paths if 'OFFL' in s]
scene_path = offl_scenes[len(offl_scenes) // 2]
url = storage_account_url + container_name + '/' + scene_path
print('Processing image at URL:\n{}'.format(url))
import warnings; warnings.filterwarnings('ignore')
with fsspec.open(url+sas_token) as f:
    ds = xr.open_dataset(f)
print(ds)<jupyter_output><xarray.Dataset>
Dimensions:  ()
Data variables:
    *empty*
Attributes: (12/56)
    Conventions:                        CF-1.7
    institution:                        DLR-IMF
    source:                             Sentinel 5 precursor, TROPOMI, space-...
    history:                            2021-01-03 06:51:09.090796 f_s5pops u...
    summary:                            TROPOMI/S5P L2 data Swath 5.5x3.5km2 ...
    tracking_id:                        c42431f4-697e-455f-8cea-92f2c6e5449d
    ...                                 ...
    title:                              TROPOMI/S5P Ozone Total Column
    Status_MET_2D:                      Nominal
    Status_NISE__:                      Nominal
    Status_L2__CLOUD_:                  External
    Status_reference_spectrum:          solar
    Status_BG:                          Empty(dtype=dtype('S1'))<jupyter_text>Open the data...which lives in the 'PRODUCT' NetCDF group.<jupyter_code>with fsspec.open(url+sas_token) as f:
    ds = xr.open_dataset(f,group='/PRODUCT')
print(ds)<jupyter_output><xarray.Dataset>
Dimensions:                                (corner: 4, ground_pixel: 450, layer: 14, level: 15, scanline: 4173, time: 1)
Coordinates:
  * scanline                               (scanline) float64 0.0 ... 4.172e+03
  * ground_pixel                           (ground_pixel) float64 0.0 ... 449.0
  * time                                   (time) datetime64[ns] 2021-01-01
  * corner                                 (corner) float64 0.0 1.0 2.0 3.0
  * layer                                  (layer) int32 0 1 2 3 ... 10 11 12 13
  * level                                  (level) int32 0 1 2 3 ... 11 12 13 14
Data variables:
    latitude                               (time, scanline, ground_pixel) float32 ...
    longitude                              (time, scanline, ground_pixel) float32 ...
    delta_time                             (time, scanline, ground_pixel) datetime64[ns] ...
    time_utc                               (time, scanline) object ...
    qa_value           [...]<jupyter_text>Plot the data in its native coordinates<jupyter_code>if 'CH4' in product:
    varname = 'methane_mixing_ratio'
elif 'NO2' in product:
    varname = 'nitrogendioxide_tropospheric_column'
elif 'O3' in product:
    varname = 'ozone_total_vertical_column'
elif 'CO' in product:
    varname = 'carbonmonoxide_total_column'
elif 'SO2' in product:
    varname = 'sulfurdioxide_total_vertical_column'
    
data = ds[varname][0,:,:]
data.plot();<jupyter_output><empty_output><jupyter_text>Extract values as numpy arrays<jupyter_code>z = data.values
lon = ds['longitude'].values.squeeze()
lat = ds['latitude'].values.squeeze()

# Don't plot extreme latitudes; they make it hard to zoom in a nice way
minlat = -60; maxlat = 60

# Zoom to a sensible longitude area by plotting only values that are non-nan
# and above a threshold.
plot_threshold = np.nanpercentile(z,50)

valid_indices = np.argwhere((~np.isnan(z)) & (z > plot_threshold))

minlon = None; maxlon = None
for xy in valid_indices:
    xy_lon = lon[xy[0],xy[1]]
    xy_lat = lat[xy[0],xy[1]]
    if xy_lat > maxlat or xy_lat < minlat:
        continue
    if minlon is None or xy_lon < minlon:
        minlon = xy_lon
    if maxlon is None or xy_lon > maxlon:
        maxlon = xy_lon

plot_extent = (minlon, maxlon, minlat, maxlat)<jupyter_output><empty_output><jupyter_text>Plot on a basemap<jupyter_code>import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
from cartopy.feature import NaturalEarthFeature

figure, ax = plt.subplots(figsize=(15, 10))
ax = plt.axes(projection=ccrs.PlateCarree(central_longitude=0.0))

# Prepare the background and axes
boundaries = cfeature.NaturalEarthFeature(
    category='cultural',name='admin_0_countries',scale='50m',facecolor='none')
ax.add_feature(boundaries, edgecolor='lightgray')
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
                  linewidth=1, color='gray', alpha=0.5, linestyle=':')
gl.xformatter = LONGITUDE_FORMATTER; gl.yformatter = LATITUDE_FORMATTER
ax.set_extent(plot_extent,ccrs.PlateCarree())

# Plot
plt.contourf(lon, lat, z, 50, transform=ccrs.PlateCarree())
plt.colorbar(fraction=0.015, pad=0.08)
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09241379310344827, 'fraction_numerical': 0.09758620689655173, 'mean_word_length': 3.491997934950955, 'pattern_counts': {'":': 0, '<': 42, '<?xml version=': 0, '>': 44, 'https://': 3, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '546129', 'license': 'MIT', 'n_tokens_mistral': 3531, 'n_tokens_neox': 2994, 'n_words': 650, 'path': 'data/sentinel-5p.ipynb', 'repo_name': 'smartlixx/AIforEarthDataSets', 'size': '8700'}