#!/usr/bin/env python
# coding: utf-8

# <img src='../../img/dust_banner.png' alt='Training school and workshop on dust' align='center' width='100%'></img>
# 
# <br>

# # CAMS global reanalysis of atmospheric composition (EAC4) - Dust Aerosol Optical Depth

# ### About

# This notebook provides an introduction to [CAMS global reanalysis of atmospheric composition (EAC4)](https://ads.atmosphere.copernicus.eu/cdsapp#!/dataset/cams-global-reanalysis-eac4?tab=overview) and shows you how the variable `Dust Aerosol Optical Depth at 550nm` can be used to better understand past dust events.
# 
# EAC4 (ECMWF Atmospheric Composition Reanalysis 4) is the fourth generation ECMWF global reanalysis of atmospheric composition. CAMS Reanalysis combines model data with observations from across the world into a globally complete and consistent dataset using a model of the atmosphere based on the laws of physics and chemistry. This principle, called data assimilation, is based on the method used by numerical weather prediction centres and air quality forecasting centres, where every so many hours (12 hours at ECMWF) a previous forecast is combined with newly available observations in an optimal way to produce a new best estimate of the state of the atmosphere, called analysis, from which an updated, improved forecast is issued. Reanalysis works in the same way to allow for the provision of a dataset spanning back more than a decade. Reanalysis does not have the constraint of issuing timely forecasts, so there is more time to collect observations, and when going further back in time, to allow for the ingestion of improved versions of the original observations, which all benefit the quality of the reanalysis product.
# 
# The notebook examines the **Saharan Dust event** which occured over Europe at the beginning of February 2021.
# 

# ### Basic Facts

# > **Spatial resolution**: `0.75° x 0.75°` <br>
# > **Spatial coverage**: `Global` <br>
# > **Temporal resolution**:  `3-hourly` <br>
# > **Temporal coverage**: `2003 to 2022` <br>
# > **Data format**: `GRIB` or `zipped NetCDF`
# 

# ### How to access the data

# CAMS global atmospheric composition reanalysis are available for download via the [Copernicus Atmosphere Data Store (ADS)](https://ads.atmosphere.copernicus.eu/). You will need to create an ADS account [here](https://ads.atmosphere.copernicus.eu/user/register).
# 
# Data from the ADS can be downloaded in two ways:
# * `manually` via the [ADS web interface](https://ads.atmosphere.copernicus.eu/cdsapp#!/dataset/cams-global-reanalysis-eac4?tab=form)
# * `programmatically` with a Python package called cdsapi ([more information](https://ads.atmosphere.copernicus.eu/api-how-to))

# ### Module outline
# * [*Optional: Retrieve CAMS global atmospheric reanalysis programmatically*](#cams_cdsapi_global)
# * [1 - Load and browse dust aerosol optical depth (duaod) at 550nm of the CAMS global reanalysis](#load_browse_cams_global)
# * [2 - Retrieve the data variable dust AOD at 550nm as xarray.DataArray](#data_retrieve_cams_global)
# * [3 - Visualize dust AOD at 550nm](#visualize_cams_global)
# * [4 - Create a geographical subset for Europe](#subset_cams_global)
# * [5 - Animate dust AOD at 550nm changes over time](#animate_cams_global)

# <hr>

# #### Load required libraries

# In[1]:


import xarray as xr

import matplotlib.pyplot as plt
import matplotlib.colors
from matplotlib.cm import get_cmap
from matplotlib import animation
from matplotlib.axes import Axes

import cartopy.crs as ccrs
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import cartopy.feature as cfeature
from cartopy.mpl.geoaxes import GeoAxes
GeoAxes._pcolormesh_patched = Axes.pcolormesh

from IPython.display import HTML

import warnings
warnings.simplefilter(action = "ignore", category = RuntimeWarning)


# #### Load helper functions

# In[2]:


get_ipython().run_line_magic('run', '../functions.ipynb')


# <hr>

# ### <a id='cams_cdsapi_global'></a>*Optional: Retrieve CAMS global atmospheric composition reanalysis programmatically*

# The `CDS Application Program Interface (CDS API)` is a Python library which allows you to access data from the ADS `programmatically`. In order to use the CDS API, follow the steps below:
# 
# * [Self-register](https://ads.atmosphere.copernicus.eu/#!/home) at the ADS registration page (if you do not have an account yet)
# * [Login](https://ads.atmosphere.copernicus.eu/user/login) to the ADS portal and go to the [api-how-to page](https://ads.atmosphere.copernicus.eu/api-how-to)
# * Copy the CDS API key displayed in the black terminal window and replace the `######` of the `KEY` variable below with your individual CDS API key
# 
# **Note:** You find your CDS API key displayed in the black terminal box under the section `Install the CDS API key`. If you do not see a URL or key appear in the black terminal box, please refresh your browser tab. 

# In[3]:


URL='https://ads.atmosphere.copernicus.eu/api/v2'
KEY='#########################'


# <br>

# The next step is then to request the data with a so called `API request`. Via the [ADS web interface](https://ads.atmosphere.copernicus.eu/cdsapp#!/dataset/cams-global-atmospheric-composition-reanalysis?tab=form), you can select the data and at the end of the web interface, you can open the ADS request via `Show API request`. Copy the request displayed there in the cell below. Once you execute the cell, the download of the data starts automatically. You will see a new file called `20210204_duaod.nc` appear in the same folder as this notebook. This is for demonstration purposes only.

# In[ ]:


import cdsapi

c = cdsapi.Client(url=URL, key=KEY)

c.retrieve(
    'cams-global-reanalysis-eac4',
    {
        'format': 'netcdf',
        'variable': 'dust_aerosol_optical_depth_550nm',
        'date': '2021-02-04/2021-02-08',
        'time': [
            '00:00', '03:00', '06:00',
            '09:00', '12:00', '15:00',
            '18:00', '21:00',
        ],
    },
    './20210204_duaod.nc')


# <br>

# <br>

# ### <a id='load_browse_cams_global'></a>1. Load and browse CAMS global atmospheric composition reanalysis

# CAMS global near-real-time forecast data is available either in `GRIB` or `netCDF`. The data for the present example has been downloaded as `netCDF`.
# 
# You can use xarray's function `xr.open_dataset()` to open the netCDF file as `xarray.Dataset`.

# In[3]:


file = xr.open_dataset('../../../eodata/dust/part1/3_model/cams/global_reanalysis/20210204_duaod.nc')
file


# <br>

# The data above has three dimensions (`latitude`, `longitude`, `time`) and one data variable:
# * `duaod550`: Dust Aerosol Optical Depth at 550nm
# 
# Let us inspect the coordinates of the file more in detail.

# Below, you see that the data set consists of 40 time steps, ranging from 4 February 2021 0 UTC to 8 February 2021 21 UTC in a 3-hour timestep.

# In[4]:


file.time


# The latitude values have a 0.75 degrees resolution and have a global N-S coverage.

# In[5]:


file.latitude


# The longitude values have a 0.75 degrees resolution as well, and are disseminated in a [0, 360] grid. 

# In[6]:


file.longitude


# Above, you see that the `longitude` variables are in the range of `[0, 359.6]`. Per default, ECMWF data are on a [0, 360] grid. Let us bring the longitude coordinates to a `[-180, 180]` grid. You can use the xarray function `assign_coords()` to assign new values to coordinates of a xarray data arra. The code below shifts your longitude coordinates from `[0, 359.6]` to `[-180, 179.6]`.

# In[7]:


file = file.assign_coords(longitude=(((file.longitude + 180) % 360) - 180)).sortby('longitude')
file


# Note that data disseminated in 2022 onwards may be already regridded onto a [-180, 180] grid, so you can skip this step if your data is already regridded.

# <br>

# ### <a id='data_retrieve_cams_global'></a>2. Retrieve the variable *Dust Aerosol Optical Depth at 550nm* as data array

# Let us extract from the dataset above the data variable `Dust Aerosol Optical Depth (AOD) at 550nm` as `xarray.DataArray` with the name `du_aod`. You can load a data array from a xarray dataset by specifying the name of the variable (`duaod550`) in square brackets.

# In[8]:


du_aod = file['duaod550']
du_aod


# Above, you see that the variable `du_aod` has two attributes, `units` and `long_name`. Let us define variables for those attributes. The variables can be used later for visualizing the data.

# In[9]:


long_name = du_aod.long_name
units = du_aod.units


# Let us do the same for the coordinates `longitude` and `latitude`.

# In[10]:


latitude = du_aod.latitude
longitude = du_aod.longitude


# <br>

# ### <a id='visualize_cams_global'></a>3. Visualize *Dust Aerosol Optical Depth at 550nm*

# The next step is to visualize the dataset. You can use the function [visualize_pcolormesh](../functions.ipynb#visualize_pcolormesh), which makes use of matploblib's function `pcolormesh` and the [Cartopy](https://scitools.org.uk/cartopy/docs/latest/) library.
# 
# With `?visualize_pcolormesh` you can open the function's docstring to see what keyword arguments are needed to prepare your plot.

# In[11]:


get_ipython().run_line_magic('pinfo', 'visualize_pcolormesh')


# You can make use of the variables we defined above:
# - `units`
# - `long_name`
# - `latitude`
# - `longitude`
# 
# Additionally, you can specify the `color scale` and minimum (`vmin`) and maxium (`vmax`) data values.

# In[12]:


time_index = 6
visualize_pcolormesh(data_array=du_aod[time_index,:,:],
                     longitude=longitude,
                     latitude=latitude,
                     projection=ccrs.PlateCarree(),
                     color_scale='hot_r',
                     unit=units,
                     long_name=long_name + ' ' + str(du_aod[time_index,:,:].time.data)[0:19],
                     vmin=0, 
                     vmax=1.5)


# <br>

# ### <a id='subset_cams_global'></a>4. Create a geographical subset for Europe

# The map above shows Dust Aerosol Optical Depth at 550nm globally. Let us create a geographical subset for Europe, in order to better analyse the Saharan dust event which impacts parts of central and southern Europe.
# 
# For geographical subsetting, you can use the function [generate_geographical_subset](../functions.ipynb#generate_geographical_subset). You can use `?generate_geographical_subset` to open the docstring in order to see the function's keyword arguments.

# In[13]:


get_ipython().run_line_magic('pinfo', 'generate_geographical_subset')


# Define the bounding box information for Europe

# In[14]:


latmin = 28.
latmax = 71.
lonmin = -22.
lonmax = 43


# Now, let us apply the function [generate_geographical_subset](../functions.ipynb#generate_geographcial_subset) to subset the `du_aod` xarray.DataArray. Let us call the new `xarray.DataArray` `du_aod_subset`.

# In[15]:


du_aod_subset = generate_geographical_subset(xarray=du_aod, 
                                             latmin=latmin, 
                                             latmax=latmax, 
                                             lonmin=lonmin, 
                                             lonmax=lonmax)
du_aod_subset


# Let us visualize the subsetted `xarray.DataArray` again. This time, you set the `set_global` kwarg to `False` and you specify the longitude and latitude bounds specified above.
# 
# Additionally, in order to have the time information as part of the title, we add the string of the datetime information to the `long_name` variable: `long_name + ' ' + str(du_aod_subset[##,:,:].time.data)[0:19]`.

# In[16]:


time_index = 14
visualize_pcolormesh(data_array=du_aod_subset[time_index,:,:],
                     longitude=du_aod_subset.longitude,
                     latitude=du_aod_subset.latitude,
                     projection=ccrs.PlateCarree(),
                     color_scale='hot_r',
                     unit=units,
                     long_name=long_name + ' ' + str(du_aod_subset[time_index,:,:].time.data)[0:19],
                     vmin=0, 
                     vmax=1,
                     lonmin=lonmin,
                     lonmax=lonmax,
                     latmin=latmin,
                     latmax=latmax,
                     set_global=False)


# <br>

# ### <a id='animate_cams_global'></a>5. Animate changes of *Dust Aerosol Optical Depth at 550nm* over time

# In the last step, you can animate the `Dust Aerosol Optical Depth at 550nm` in order to see how the trace gas develops over a period of 5 days, from 4th to 8th February 2021.
# 
# You can do animations with matplotlib's function `animation`. Jupyter's function `HTML` can then be used to display HTML and video content.

# The animation function consists of 4 parts:
# - **Setting the initial state:**<br>
#  Here, you define the general plot your animation shall use to initialise the animation. You can also define the number of frames (time steps) your animation shall have.
#  
#  
# - **Functions to animate:**<br>
#  An animation consists of three functions: `draw()`, `init()` and `animate()`. `draw()` is the function where individual frames are passed on and the figure is returned as image. In this example, the function redraws the plot for each time step. `init()` returns the figure you defined for the initial state. `animate()` returns the `draw()` function and animates the function over the given number of frames (time steps).
#  
#  
# - **Create a `animate.FuncAnimation` object:** <br>
#  The functions defined before are now combined to build an `animate.FuncAnimation` object.
#  
#  
# - **Play the animation as video:**<br>
#  As a final step, you can integrate the animation into the notebook with the `HTML` class. You take the generate animation object and convert it to a HTML5 video with the `to_html5_video` function

# In[17]:


# Setting the initial state:
# 1. Define figure for initial plot
fig, ax = visualize_pcolormesh(data_array=du_aod_subset[0,:,:],
                               longitude=du_aod_subset.longitude, 
                               latitude=du_aod_subset.latitude,
                               projection=ccrs.PlateCarree(), 
                               color_scale='hot_r', 
                               unit='-',
                               long_name=long_name + ' '+ str(du_aod_subset.time[0].data)[0:19],
                               vmin=0,
                               vmax=1, 
                               lonmin=lonmin, 
                               lonmax=lonmax, 
                               latmin=latmin, 
                               latmax=latmax,
                               set_global=False)

frames = 30

def draw(i):
    img = plt.pcolormesh(du_aod_subset.longitude, 
                         du_aod_subset.latitude, 
                         du_aod_subset[i,:,:], 
                         cmap='hot_r', 
                         transform=ccrs.PlateCarree(),
                         vmin=0,
                         vmax=1,
                         shading='auto')
    
    ax.set_title(long_name + ' '+ str(du_aod_subset.time[i].data)[0:19], fontsize=20, pad=20.0)
    return img

def init():
    return fig

def animate(i):
    return draw(i)

ani = animation.FuncAnimation(fig, animate, frames, interval=800, blit=False,
                              init_func=init, repeat=True)

HTML(ani.to_html5_video())
plt.close(fig)


# <br>

# #### Play the animation video as HTML5 video

# In[18]:


HTML(ani.to_html5_video())


# <hr>

# ### References
# 
# * Copernicus Service information 2021
# * Generated using Copernicus Atmosphere Monitoring Service Information 2021

# <hr>

# <img src='../../img/copernicus_logo.png' alt='Logo EU Copernicus' align='right' width='20%'><br><br><br><br>
# <p style="text-align:right;">This project is licensed under <a href="../../LICENSE">GNU General Public License v3.0 only</a> and is developed under a Copernicus contract.

# <p style="text-align:right;"> <a href='https://training.eumetsat.int'>EUMETSAT Training</a> | <a href='mailto:training@eumetsat.int'>Contact the training team</a></p>