<jupyter_start><jupyter_text>31 - Covid-19 case study - GOME-2 anomaly map >>  Copernicus Sentinel-5P TROPOMI Carbonmonoxide (CO) A precursor satellite mission, Sentinel-5P aims to fill in the data gap and provide data continuity between the retirement of the Envisat satellite and NASA's Aura mission and the launch of Sentinel-5. The Copernicus Sentinel-5P mission is being used to closely monitor the changes in air quality and was launched in October 2017.Sentinel-5p Pre-Ops data are disseminated in the `netCDF` format and can be downloaded via the [Copernicus Open Access Hub](https://scihub.copernicus.eu/).Sentinel-5p carries the `TROPOMI` instrument, which is a spectrometer in the UV-VIS-NIR-SWIR spectral range. `TROPOMI` provides measurements on:* `Ozone`* `NO``2`* `SO``2`* `Formaldehyde`* `Aerosol`* `Carbonmonoxide`* `Methane`* `Clouds`  Module outline:* [1 - Load and browse Sentinel-5P data](load_s5p)* [2 - Plotting example - Sentinel-5P data](plotting_s5p)  Load required libraries<jupyter_code>%matplotlib inline
import os
import xarray as xr
import numpy as np
import netCDF4 as nc

import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
import cartopy.crs as ccrs
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER

from matplotlib.axes import Axes
from cartopy.mpl.geoaxes import GeoAxes
GeoAxes._pcolormesh_patched = Axes.pcolormesh
import cartopy.feature as cfeature

import geopandas as gpd

import warnings
warnings.simplefilter(action = "ignore", category = RuntimeWarning)<jupyter_output><empty_output><jupyter_text>Load and browse Sentinel-5P data  Open one individual Sentinel-5P netCDF file with `NetCDF4` The dataset object contains information about the general data structure of the dataset. You can see that the variables of `Sentinel-5P` data are organised in groups, which is analogous to directories in a filesystem.<jupyter_code>s5p_file = nc.Dataset('./eodata/sentinel5p/co/2019/08/19/S5P_OFFL_L2__CO_____20190819T164807_20190819T182937_09581_01_010302_20190825T161022.nc', 'r')
s5p_file.groups<jupyter_output><empty_output><jupyter_text>If you select the `/PRODUCT` group, you get more information on what variables the dataset object contain.<jupyter_code>s5p_file.groups['PRODUCT']<jupyter_output><empty_output><jupyter_text>You see that the object contains the following variables:* `scanline`* `ground_pixel`* `time`* `corner`* `delta_time`* `time_utc`* `ga_value`* `latitude`* `longitude`* `carbonmonoxide_total_column`* `carbonmonoxie_total_column_precision` You can specify one variable of interest and get more detailed information about the variable. E.g. `carbonmonoxide_total_column` is the atmosphere mole content of carbon monoxide, has the unit mol m-2, and is a 3D variable.You can do this for the available variables, but also for the dimensions latitude and longitude.You can see e.g. that the `latitude` coordinates range between -85.9 S and 61.9 S and the `longitude` coordinates range between -124.3 W to 101.9 E.<jupyter_code>co = s5p_file.groups['PRODUCT'].variables['carbonmonoxide_total_column']
lon = s5p_file.groups['PRODUCT'].variables['longitude'][:][0,:,:]
lat = s5p_file.groups['PRODUCT'].variables['latitude'][:][0,:,:]

co, lon, lat<jupyter_output><empty_output><jupyter_text>You can retrieve the array values of the variable object by selecting the `time` dimension and `data`. You can have a look at the `minimum` and `maximum` data value to get an idea of the data range. You see that the data contain negative values. Let's mask the negative values and all values equal to the `_FillValue` and set it to `NaN`. `_FillValue` is used for not significant data. Thus, you want to mask those.<jupyter_code>co_data = co[0,:,:].data
print(co_data.min(), co_data.max())
co_data[co_data <= 0.] = co._FillValue
co_data[co_data == co._FillValue] = np.nan<jupyter_output><empty_output><jupyter_text>Plotting example - Sentinel-5P data  Plot `Dataset` NetCDF library object with `matplotlib` and `cartopy` The retrieved data array from the Dataset NetCDF object is of type `numpy array` and you can plot it with matplotlib's `pcolormesh` function. Due to the nature of the `CO` data values, we apply a logarithmic scale to the color bar with `LogNorm` from `matplotlib.colors`, which facilitates the visualisation of the data.Let's create a function [visualize_pcolormesh](./functions.ipynbvisualize_pcolormesh), where we can specify projection, extent, conversion_factor, color_scale, unit, title and if the plot shall have a global extent.<jupyter_code>def visualize_pcolormesh(data_array, longitude, latitude, projection, color_scale, unit, long_name, vmin, vmax, lonmin, lonmax, latmin, latmax, log=True, set_global=True):
    """ 
    Visualizes a numpy array with matplotlib's 'pcolormesh' function.
    
    Parameters:
        data_array: any numpy MaskedArray, e.g. loaded with the NetCDF library and the Dataset function
        longitude: numpy Array holding longitude information
        latitude: numpy Array holding latitude information
        projection: a projection provided by the cartopy library, e.g. ccrs.PlateCarree()
        color_scale (str): string taken from matplotlib's color ramp reference
        unit (str): the unit of the parameter, taken from the NetCDF file if possible
        long_name (str): long name of the parameter, taken from the NetCDF file if possible
        vmin (int): minimum number on visualisation legend
        vmax (int): maximum number on visualisation legend
        lonmin,lonmax,latmin,latmax: geographic extent of the plot
        log (logical): set True, if the values shall be represented in a logarithmic scale
        set_global (logical): set True, if the plot shall have a global coverage
    """
    fig=plt.figure(figsize=(20, 10))

    ax = plt.axes(projection=projection)

    # define the coordinate system that the grid lons and grid lats are on
    if(log):
        img = plt.pcolormesh(longitude, latitude, np.squeeze(data_array), norm=LogNorm(), 
                             cmap=plt.get_cmap(color_scale), transform=ccrs.PlateCarree(),
                            vmin=vmin,
                            vmax=vmax)
    else:
        img = plt.pcolormesh(longitude, latitude, data_array, 
                        cmap=plt.get_cmap(color_scale), transform=ccrs.PlateCarree(),
                        vmin=vmin,
                        vmax=vmax)

    ax.add_feature(cfeature.BORDERS, edgecolor='black', linewidth=1)
    ax.add_feature(cfeature.COASTLINE, edgecolor='black', linewidth=1)

    if (projection==ccrs.PlateCarree()):
        ax.set_extent([lonmin, lonmax, latmin, latmax], projection)
        gl = ax.gridlines(draw_labels=True, linestyle='--')
        gl.xlabels_top=False
        gl.ylabels_right=False
        gl.xformatter=LONGITUDE_FORMATTER
        gl.yformatter=LATITUDE_FORMATTER
        gl.xlabel_style={'size':14}
        gl.ylabel_style={'size':14}

    if(set_global):
        ax.set_global()
        ax.gridlines()

    cbar = fig.colorbar(img, ax=ax, orientation='horizontal', fraction=0.04, pad=0.1)
    cbar.set_label(unit, fontsize=16)
    cbar.ax.tick_params(labelsize=14)
    ax.set_title(long_name, fontsize=20, pad=20.0)

    return fig, ax<jupyter_output><empty_output><jupyter_text>You can retrieve unit and title information from the load `Dataset`, where the information is stored as attributes. You can now plot the data.<jupyter_code>unit = co.units
long_name = co.long_name

visualize_pcolormesh(co_data, lon, lat, ccrs.Mollweide(), 'viridis', unit, long_name, 0.01,1, lon.min(), lon.max(), lat.min(), lat.max(), log=True, set_global=True)<jupyter_output><empty_output><jupyter_text>You can zoom into a region by specifying a `bounding box` of interest. Let's set the extent to South America, with: `[-100, 0, -80, 40]`. The above plotting function [visualize_pcolormesh](./functions.ipynbvisualize_pcolormesh) allows for setting a specific bounding box. You simply have to set the `set_global` key to False. It is best to adjust the projection to `PlateCarree()`, as this will be more appropriate for a regional subset.<jupyter_code>lonmin=-100
lonmax=0
latmin=-80
latmax=40

visualize_pcolormesh(co_data, lon, lat, ccrs.PlateCarree(), 'viridis', unit, long_name, 0.01,1, lonmin, lonmax, latmin, latmax, log=True, set_global=False)<jupyter_output><empty_output><jupyter_text>Load multiple Sentinel-5p data files with `xarray` and `open_mfdataset` The plots above showed the extent of one Sentinel-5P ground track. You can load multiple ground tracks into a single `xarray` and the `DataArrays` will be concatenated at the `scanline` dimension. This allows you to have a larger region of interest (ROI).<jupyter_code>s5p_mf_19 = xr.open_mfdataset('./eodata/sentinel5p/co/2019/08/19/*.nc', concat_dim='scanline', combine='nested', group='PRODUCT')
s5p_mf_19<jupyter_output><empty_output><jupyter_text>From the `Dataset` object `s5p_mf_19`, you can choose the data variable of interest, e.g. `carbonmonoxide_total_column`. It has three dimensions (`3D`), but the time dimension consisting only of one dimension entry. If you want to reduce it by the dimension time, you can simply select the first dimension and reduce it to a `2D` object. You can again use the function [visualize_pcolormesh](./functions.ipynbvisualize_pcolormesh) to visualize the data.<jupyter_code>co_19 = s5p_mf_19.carbonmonoxide_total_column[0,:,:]
lat_19 = co_19.latitude
lon_19 = co_19.longitude

unit = co_19.units
long_name = co_19.long_name
visualize_pcolormesh(co_19, lon_19, lat_19, ccrs.PlateCarree(), 'viridis', unit, long_name, 0.01, 1.0, lonmin, lonmax, latmin, latmax, log=True, set_global=False)<jupyter_output><empty_output>