<jupyter_start><jupyter_text>WMO SDS-WAS Barcelona Dust Forecast - MONARCH  About This notebook provides an introduction to dust forecast data from the NMMB/MONARCH model. The notebook introduces you to the variable `Dust Optical Depth` and you will learn how the model has predicted the **Saharan Dust event** which occured over Europe in the first half of February 2021.The [WMO Sand and Dust Storm Warning Advisory and Assessment System (SDS-WAS)](https://community.wmo.int/activity-areas/gaw-sand-and-dust-storm-warning-advisory-and-assessment-system-sds-was) is an international framework linking institutions involved in [Sand and Dust Storm (SDS)](https://public.wmo.int/en/our-mandate/focus-areas/environment/SDS) research, operations and delivery of services. The SDS-WAS has the following objectives:* Provide user communities access to forecasts, observations and information of the SDS through regional centres connected to the WMO Information System (WIS) and the World Wide Web* Identify and improve SDS products through consultation with the operational and user communities* Enhance operational SDS forecasts through technology transfer from research* Improve forecasting and observation technology through coordinated international research and assessment* Build capacity of relevant countries to utilize SDS observations, forecasts and analysis products for meeting societal needs* Build bridges between SDS-WAS and other communities conducting aerosol related studies (air quality, biomass burning, etc.) The framework is organised in several regional centers, which aim to implement SDS-WAS objectives in a specific region. The Barcelona Supercomputing Center (BSC-CNS) and the Meteorological State Agency of Spain (AEMET) are hosting the [SDS-WAS regional center for Northern Africa, Middle East and Europe](https://sds-was.aemet.es/).One of the main activities of the SDS-WAS regional center is to provide daily operational dust forecasts for Northern Africa (north of the equator), Middle East and Europe. The BCS-CNS, in collaboration with NOAA's National Centers for Environmental Prediction (NCEP), the NASA's Goddard Institute for Space Studies and the International Research Institute for Climate Society (IRI), has developed [NMMB/MONARCH](https://dust.aemet.es/methods/the-nmmb-bsc-dust-model), an online multi-scale atmospheric dust model intended to provide short and medium-range dust forecasts for both, regional and global domains. The model provides forecast information up to 72 hours in advance (every 3 hours) of two parameters: `Dust Optical Depth` and `Dust Surface Concentration`.NMMB/MONARCH Forecast data are available in `netCDF` format and are available for download [here](https://sds-was.aemet.es/forecast-products/dust-forecasts/files-download/nmmb-bsc-dust-public/).  Basic Facts > **Spatial resolution**: `0.1° x 0.1°` > **Spatial coverage**: `Northern Africa, Middle East and Europe` > **Temporal resolution**:  `3-hourly up to 72 hours in advance` > **Temporal coverage**: `since February 2012` > **Data format**: `NetCDF`  How to access the data Dust forecast data from the NMMB/MONARCH model are available for download via the website of the [SDS-WAS Regional Center for Northern Africa, Middle East and Europe](https://sds-was.aemet.es/forecast-products/dust-forecasts/files-download).Every participant of the training school and workshop on dust monitoring can use the following credentials for downloading data:* User: `<EMAIL>`* Password: `<PASSWORD>`  Module outline* [1 - Load and browse NMMB/MONARCH dust forecast data](load_file)* [2 - Load the data variables as data arrays](load_variables_sdswas)* [3 - Visualize dust Aerosol Optical Depth](visualize_bcs_dust)* [4 - Animate dust AOD forecasts](animate_bcs_dust)   Load required libraries<jupyter_code>%matplotlib inline
import os
import xarray as xr
import numpy as np
import netCDF4 as nc
import pandas as pd

from IPython.display import HTML

import matplotlib.pyplot as plt
import matplotlib.colors
from matplotlib.cm import get_cmap
from matplotlib import animation
import cartopy.crs as ccrs
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import cartopy.feature as cfeature

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

import warnings
warnings.simplefilter(action = "ignore", category = RuntimeWarning)<jupyter_output><empty_output><jupyter_text>Load helper functions<jupyter_code>%run ../functions.ipynb<jupyter_output><empty_output><jupyter_text>1. Load and browse NMMB/MONARCH dust forecast data The first step is to load a `NMMB/MONARCH` forecast file to better understand its structure. The data is disseminated in the `netCDF` format on a daily basis, with the forecast initialisation at 12:00 UTC. Let us load the `NMMB/MONARCH` dust forecast of 5 February 2021. You can use the function `open_dataset()` from the xarray Python library, which makes reading a `netCDF` file very efficient and easy. The function loads the netCDF file as `xarray.Dataset()`, which is a data collection of multiple variables sharing the same coordinate information.Below, you see that the `NMMB/MONARCH` data has three dimensions, `lat`, `lon` and `time`, and offers two data variables `od550_dust` and `sconc_dust`.<jupyter_code>filepath = '../../../eodata/50_modules/01_dust/03_model/sds_was/2021020512_3H_NMMB-BSC.nc'
file = xr.open_dataset(filepath)
file<jupyter_output><empty_output><jupyter_text>Let us now have a closer look at the dimensions of the data. Let us first inspect the two coordinate dimensions `lat` and `lon`. You can simply access the dimension's `xarray.DataArray` by specifying the name of the dimension. Below you see that the data has a 0.1 x 0.1 degrees resolution and have the following geographical bounds:* `Longitude: [-63., 101.9]`* `Latitude: [-11., 71.4]`<jupyter_code>latitude = file.lat
longitude = file.lon
latitude, longitude<jupyter_output><empty_output><jupyter_text>Now, let us also inspect the time dimension. You see that one daily forecast file has 25 time steps, with three hourly forecast information up to 72 hours (3 days) in advance.<jupyter_code>file.time<jupyter_output><empty_output><jupyter_text>2. Load data variables as data arrays A `xarray.Dataset` is a collection of multiple variables and offers a general overview of the data, but does not offer direct access to the data arrays. For this reason, you want to load a data variable as `xarray.DataArray`. You can access the data array information by simply specifying the name of the variable after the name of the `xarray.Dataset`Let us load both variables as `xarray.DataArray`, `od550_dust` and `sconc_dust`. Below you see that, once the data array is loaded, that each variable is disseminated with additional attributes, e.g. `long_name` and `units`. These attributes give us more information about the data. `od550_dust` e.g. is Aerosol Optical Depth which does not have a unit and `sconc_dust` is the surface dust concentration measured in kg / m3.<jupyter_code>od_dust = file.od550_dust
od_dust
sconc_dust = file.sconc_dust
sconc_dust<jupyter_output><empty_output><jupyter_text>You can define variables for the attributes of a variable. This can be helpful during data visualization, as the attributes `long_name` and `units` can be added as additional information to the plot. From the `xarray.DataArray`, you simply specify the name of the attribute of interest.<jupyter_code>long_name=od_dust.long_name
units= od_dust.units<jupyter_output><empty_output><jupyter_text>3. Visualize dust Aerosol Optical Depth at 550 nm Now we have loaded all necessary information to be able to visualize the dust Aerosol Optical Depth for one specific time during the forecast run. Let us use the function `visualize_pcolormesh()` to visualize the data with the help of the plotting library [matplotlib](https://matplotlib.org/) and [Cartopy](https://scitools.org.uk/cartopy/docs/latest/).You have to specify the following keyword arguments:* `data_array`: the* `longitude`, `latitude`: longitude and latitude variables of the data variable* `projection`: one of [Cartopy's projection options](https://scitools.org.uk/cartopy/docs/v0.15/crs/projections.html)* `color_scale`: one of [matplotlib's colormaps](https://matplotlib.org/stable/tutorials/colors/colormaps.html)* `units`: unit of the data parameter, preferably taken from the data array's attributes* `long_name`: longname of the data parameter is taken as title of the plot, preferably taken from the data array's attributes* `vmin`, `vmax`: minimum and maximum bounds of the color range* `set_global`: False, if data is not global* `lonmin`, `lonmax`, `latmin`, `latmax`: kwargs have to be specified, if `set_global=False` Note: 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(od_dust[,:,:].time.data)`.<jupyter_code>forecast_step = 6
visualize_pcolormesh(data_array=od_dust[forecast_step,:,:],
                     longitude=od_dust.lon,
                     latitude=od_dust.lat,
                     projection=ccrs.PlateCarree(),
                     color_scale='hot_r',
                     unit=units,
                     long_name=long_name + ' ' + str(od_dust[forecast_step,:,:].time.data),
                     vmin=0, 
                     vmax=1.5,
                     set_global=False,
                     lonmin=longitude.min().data,
                     lonmax=longitude.max().data,
                     latmin=latitude.min().data,
                     latmax=latitude.max().data)<jupyter_output><empty_output><jupyter_text>4. Animate dust Aerosol Optical Depth forecasts In the last step, you can animate the `Dust Aerosol Optical Depth` forecasts in order to see how the trace gas develops over a period of 3 days, from 5th Feb 12 UTC 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:** 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:** 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:**  The functions defined before are now combined to build an `animate.FuncAnimation` object.  - **Play the animation as video:** 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<jupyter_code># Setting the initial state:
# 1. Define figure for initial plot
fig, ax = visualize_pcolormesh(data_array=od_dust[0,:,:],
                     longitude=od_dust.lon,
                     latitude=od_dust.lat,
                     projection=ccrs.PlateCarree(),
                     color_scale='hot_r',
                     unit=units,
                     long_name=long_name + ' ' + str(od_dust[0,:,:].time.data),
                     vmin=0, 
                     vmax=1.5,
                     lonmin=longitude.min().data,
                     lonmax=longitude.max().data,
                     latmin=latitude.min().data,
                     latmax=latitude.max().data,
                     set_global=False)

frames = 24

def draw(i):
    img = plt.pcolormesh(od_dust.lon, 
                         od_dust.lat, 
                         od_dust[i,:,:], 
                         cmap='hot_r', 
                         transform=ccrs.PlateCarree(),
                         vmin=0,
                         vmax=1.5,
                         shading='auto')
    
    ax.set_title(long_name + ' '+ str(od_dust.time[i].data), 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)
HTML(ani.to_html5_video())<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08157457908465734, 'fraction_numerical': 0.010196822385582168, 'mean_word_length': 3.9732704402515724, 'pattern_counts': {'":': 0, '<': 39, '<?xml version=': 0, '>': 44, 'https://': 10, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '436378', 'license': 'MIT', 'n_tokens_mistral': 3683, 'n_tokens_neox': 3391, 'n_words': 1401, 'path': '50_modules/01_dust/03_model-based_data/31_SDS-WAS_MONARCH_dust_forecast.ipynb', 'repo_name': 'trivedi-c/atm_Practical3', 'size': '12674'}