# Solution 1¶

You are working for ENAIRE, the air navigation authority for Spain and western Africa. You know that the Canary Islands are prone to Saharan dust events and for this reason, ENAIRE monitors dust on a daily basis. You are the operations analyst for the week of 21-27 February 2020 and responsible for issuing alerts, and if necessary, to mandate required safety measures.

On 21 February 2020, you are in-charge of analysing the dust forecast and to monitor potential dust events for the coming days. With your new knowledge on aerosol and dust data, you should be able to do this.

1. Brainstorm

• What dust forecasts do you know about?

• How do they differ from each other?

• What satellite data do you know about that can be used for dust nowcasting?

• Which variables can be used to monitor and forecast dust?

3. Download the MSG SEVIRI Level 1.5 data and visualize the Dust RGB composite

• Based on the dust forecast for the next days - which day and hour would you choose for getting a near real-time monitoring of dust from the MSG SEVIRI instrument?

• Hint

4. Interpret the results

• Describe the dust forecast event in comparison with the near real-time monitoring from the satellite.

• What decision as ENAIRE operations analyst do you take? Would you issue an alert? Would you implement some safety measures?

import xarray as xr
import glob

from IPython.display import HTML

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 satpy.scene import Scene
from satpy.composites import GenericCompositor
from satpy.writers import to_image
from satpy.resample import get_area_def

import pyresample

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


%run ../../functions.ipynb


### 1. Brainstorm¶

Model forecasts

Forecast models differ in their spatial and temporal resolution, but also spatial coverage and forecast period.

Satellite observations

• MSG SEVIRI Level 1 RGB composites

• Terra/Aqua MODIS Level 1 RGB composites

The following variables can be used to monitor and forecast dust:

• Dust Optical Depth

• Aerosol Optical Depth

The first step is to load a MONARCH forecast file. The data is disseminated in the netCDF format on a daily basis, with the forecast initialisation at 12:00 UTC. Load the MONARCH dust forecast of 21 February 2020. You can use the function open_dataset() from the xarray Python library.

Once loaded, you see that the data has three dimensions: lat, lon and time; and offers two data variables od550_dust and sconc_dust.

filepath = '../../eodata/case_study/sds_was/2020022112_3H_NMMB-BSC.nc'
file = xr.open_dataset(filepath)
file

<xarray.Dataset>
Dimensions:     (lon: 307, lat: 211, time: 25)
Coordinates:
* lon         (lon) float64 -31.0 -30.67 -30.33 -30.0 ... 70.33 70.67 71.0
* lat         (lat) float64 -3.0 -2.667 -2.333 -2.0 ... 66.0 66.33 66.67 67.0
* time        (time) datetime64[ns] 2020-02-21T12:00:00 ... 2020-02-24T12:0...
Data variables:
od550_dust  (time, lat, lon) float32 ...
sconc_dust  (time, lat, lon) float32 ...
Attributes:
CDI:                       Climate Data Interface version 1.5.4 (http://c...
Conventions:               CF-1.2
history:                   Fri Feb 21 23:50:54 2020: cdo remapbil,regular...
_FillValue:                -32767.0
missing_value:             -32767.0
title:                     Regional Reanalysis 0.5x0.5 deg NMMB-BSC-Dust ...
History:                   Fri Feb 21 22:12:45 2020: ncrcat -F -O pout_re...
Grid_type:                 B-grid: vectors interpolated to scalar positions
Map_Proj:                  Rotated latitude longitude
NCO:                       4.0.8
CDO:                       Climate Data Operators version 1.5.4 (http://c...

Let us then retrieve the data variable od550_dust, which is the dust optical depth at 550 nm.

od550_dust = file['od550_dust']
od550_dust

<xarray.DataArray 'od550_dust' (time: 25, lat: 211, lon: 307)>
[1619425 values with dtype=float32]
Coordinates:
* lon      (lon) float64 -31.0 -30.67 -30.33 -30.0 ... 70.0 70.33 70.67 71.0
* lat      (lat) float64 -3.0 -2.667 -2.333 -2.0 ... 66.0 66.33 66.67 67.0
* time     (time) datetime64[ns] 2020-02-21T12:00:00 ... 2020-02-24T12:00:00
Attributes:
long_name:  dust optical depth at 550 nm
units:
title:      dust optical depth at 550 nm

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]

latitude = file.lat
longitude = file.lon
latitude, longitude

(<xarray.DataArray 'lat' (lat: 211)>
array([-3.      , -2.666667, -2.333333, ..., 66.333326, 66.66666 , 66.999993])
Coordinates:
* lat      (lat) float64 -3.0 -2.667 -2.333 -2.0 ... 66.0 66.33 66.67 67.0
Attributes:
standard_name:  latitude
long_name:      latitude
units:          degrees_north
axis:           Y,
<xarray.DataArray 'lon' (lon: 307)>
array([-31.      , -30.666667, -30.333333, ...,  70.333323,  70.666657,
70.99999 ])
Coordinates:
* lon      (lon) float64 -31.0 -30.67 -30.33 -30.0 ... 70.0 70.33 70.67 71.0
Attributes:
standard_name:  longitude
long_name:      longitude
units:          degrees_east
axis:           X)


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.

file.time

<xarray.DataArray 'time' (time: 25)>
array(['2020-02-21T12:00:00.000000000', '2020-02-21T15:00:00.000000000',
'2020-02-21T18:00:00.000000000', '2020-02-21T21:00:00.000000000',
'2020-02-22T00:00:00.000000000', '2020-02-22T03:00:00.000000000',
'2020-02-22T06:00:00.000000000', '2020-02-22T09:00:00.000000000',
'2020-02-22T12:00:00.000000000', '2020-02-22T15:00:00.000000000',
'2020-02-22T18:00:00.000000000', '2020-02-22T21:00:00.000000000',
'2020-02-23T00:00:00.000000000', '2020-02-23T03:00:00.000000000',
'2020-02-23T06:00:00.000000000', '2020-02-23T09:00:00.000000000',
'2020-02-23T12:00:00.000000000', '2020-02-23T15:00:00.000000000',
'2020-02-23T18:00:00.000000000', '2020-02-23T21:00:00.000000000',
'2020-02-24T00:00:00.000000000', '2020-02-24T03:00:00.000000000',
'2020-02-24T06:00:00.000000000', '2020-02-24T09:00:00.000000000',
'2020-02-24T12:00:00.000000000'], dtype='datetime64[ns]')
Coordinates:
* time     (time) datetime64[ns] 2020-02-21T12:00:00 ... 2020-02-24T12:00:00
Attributes:
standard_name:  time

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.

long_name=od550_dust.long_name
units= od550_dust.units


#### Visualize dust Aerosol Optical Depth¶

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 and Cartopy.

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

• color_scale: one of matplotlib’s colormaps

• 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).

forecast_step=0
visualize_pcolormesh(data_array=od550_dust[forecast_step,:,:],
longitude=od550_dust.lon,
latitude=od550_dust.lat,
projection=ccrs.PlateCarree(),
color_scale='afmhot_r',
unit=units,
long_name=long_name + ' ' + str(od550_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)

(<Figure size 1440x720 with 2 Axes>,
<GeoAxesSubplot:title={'center':'dust optical depth at 550 nm 2020-02-21T12:00:00.000000000'}>)


#### 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 20th Feb 12 UTC to 24th 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

# Setting the initial state:
# 1. Define figure for initial plot
fig, ax = visualize_pcolormesh(data_array=od550_dust[0,:,:],
longitude=od550_dust.lon,
latitude=od550_dust.lat,
projection=ccrs.PlateCarree(),
color_scale='afmhot_r',
unit=units,
long_name=long_name + ' ' + str(od550_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(od550_dust.lon,
od550_dust.lat,
od550_dust[i,:,:],
cmap='afmhot_r',
transform=ccrs.PlateCarree(),
vmin=0,
vmax=1.5,

ax.set_title(long_name + ' '+ str(od550_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())


### 3. Download the MSG SEVIRI Level 1.5 data and visualize the Dust RGB composite¶

From the EUMETSAT Data Store, we downloaded High Rate SEVIRI Level 1.5 Image data for 23 February 2020 at 18 UTC.

The data is downloaded as a zip archive. We already unzipped the file. Hence, in a next step we can specify the file path and create a variable with the name file_name.

file_name = glob.glob('../../eodata/case_study/meteosat/event/2020/02/23/MSG4-SEVI-MSG15-0100-NA-20200223181244.112000000Z-NA.nat')
file_name

['../../eodata/case_study/meteosat/event/2020/02/23/MSG4-SEVI-MSG15-0100-NA-20200223181244.112000000Z-NA.nat']


In a next step, we use the Scene constructor from the satpy library. Once loaded, a Scene object represents a single geographic region of data, typically at a single continuous time range.

You have to specify the two keyword arguments reader and filenames in order to successfully load a scene. As mentioned above, for MSG SEVIRI data in the Native format, you can use the seviri_l1b_native reader.

scn = Scene(reader="seviri_l1b_native",
filenames=file_name)
scn

<satpy.scene.Scene at 0x7f7c2c5fb340>


Load the RGB composite ID natural color

Let us define load the composite ID natural_color. This list variable can then be passed to the function load(). Per default, scenes are loaded with the north pole facing downwards. You can specify the keyword argument upper_right_corner=NE in order to turn the image around and have the north pole facing upwards.

composite_id = ["natural_color"]


A print of the Scene object scn shows you that one band is available: natural_color

print(scn)

<xarray.DataArray 'where-1eeee2327f1f52a84787d02280871d42' (bands: 3, y: 3712, x: 3712)>
dask.array<where, shape=(3, 3712, 3712), dtype=float32, chunksize=(1, 3712, 3712), chunktype=numpy.ndarray>
Coordinates:
crs      object PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["unknown",E...
* y        (y) float64 5.569e+06 5.566e+06 5.563e+06 ... -5.563e+06 -5.566e+06
* x        (x) float64 -5.569e+06 -5.566e+06 ... 5.563e+06 5.566e+06
* bands    (bands) <U1 'R' 'G' 'B'
Attributes:
end_time:                               2020-02-23 18:15:11.171248
sun_earth_distance_correction_factor:   0.9787042740413612
start_time:                             2020-02-23 18:00:11.057487
area:                                   Area ID: msg_seviri_fes_3km\nDesc...
georef_offset_corrected:                True
orbital_parameters:                     {'projection_longitude': 0.0, 'pr...
sensor:                                 seviri
ancillary_variables:                    []
standard_name:                          natural_color
resolution:                             3000.403165817
platform_name:                          Meteosat-11
sun_earth_distance_correction_applied:  True
wavelength:                             None
optional_datasets:                      []
name:                                   natural_color
_satpy_id:                              DataID(name='natural_color', reso...
prerequisites:                          [DataQuery(name='IR_016', modifie...
optional_prerequisites:                 []
mode:                                   RGB


Generate a geographical subset over Canary Islands

Crop the data with x and y values in original projection unit (x_min, y_min, x_max, y_max). Below we apply the funcion crop() and use the keyword argument xy_bbox=(xmin, ymin, xmax, ymax) in order to crop the file based on an area in original projection unit. For the bbox information we use the following settings: xmin=-20E5, ymin=23E5, xmax=-5E5, ymax=35E5.

Then, let’s visualize the cropped image with the function show().

scn_cropped = scn.crop(xy_bbox=(-20E5, 23E5, -5E5, 35E5))
scn_cropped.show("natural_color")


Let us store the area definition as variable.

area = scn_cropped['natural_color'].attrs['area']
area


Next, use the area definition to resample the loaded Scene object. Afterwards, you can visualize the resampled image with the function show().

scn_resample_nc = scn.resample(area)
scn_resample_nc.show('natural_color')


Load, visualize and interpret the MSG SEVIRI Dust RGB composite

Let us reload the file with the Scene constructor.

scn2 = Scene(reader="seviri_l1b_native",
filenames=file_name)
scn2

<satpy.scene.Scene at 0x7f7d387b6bb0>


Let us define a composite ID for dust in a list. This list variable can then be passed to the function load(). Per default, scenes are loaded with the north pole facing downwards. You can specify the keyword argument upper_right_corner=NE in order to turn the image around and have the north pole facing upwards.

composite_id = ["dust"]

scn2_resample_nc = scn2.resample(area)

Afterwards, you can visualize the resampled image with the function show().
scn2_resample_nc.show('dust')