This notebook provides you an introduction to data from the SEVIRI instrument as part of the 0 Degree Service as well an introduction to the Python package satpy which allows for an efficient handling of data from meteorological satellite instruments, including SEVIRI and MODIS. At the end, a specific focus will be put on the SEVIRI Dust RGB, whose primary aim is to detect dust in the atmosphere.
The example features the Saharan dust event during 5./6. February 2021, which brought massive amounts of Saharan dust into central Europe. See a more detailed description of this particular event here.
The Meteosat series have been providing crucial data for weather forecasting since 1977 and is a series of geostationary satellites providing imagery for weather forecasting and climate monitoring. Meteosat Second Genersation (MSG) is the current fleet of operational geostationary satellites. The Spinning Enhanced Visible and InfraRed Imager (SEVIRI) is MSG's primary instrument and has the capacity to observe the Earth in 12 spectral channels. 11 channels provides measurements with a resolution of 3 km at the sub-satellite and one, the High Resolution Visible (HRV) channel, provides measurements with a resolution of 1 km.
The SEVIRI instrument allows for a complete image scan (Full Earth Scan) every 15 minutes. The 0 Degree Service is the main mission of MSG and provides High Rate SEVIRI image data in 12 spectral bands, processed in near real-time to Level 1.5.
Spatial resolution:
1 km at nadir
Spatial coverage:Latitude: -81 to 81 degrees
,Longitude: -79 to 79 degrees
Revisit time:every 15 minutes
Data availability:since 2004
High Rate SEVIRI Level 1.5 Image Data - MSG - 0 degree data is available for download via the EUMETSAT Data Store.
You need to register for the EUMETSAT Earth Observation Portal in order to be able to download data from the Data Store.
import glob
import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors
from matplotlib.axes import Axes
import cartopy
import cartopy.crs as ccrs
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
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
from satpy import available_readers
import pyresample
import warnings
warnings.filterwarnings('ignore')
warnings.simplefilter(action = "ignore", category = RuntimeWarning)
%run ../functions.ipynb
Meteosat Second Generation data are disseminated in the specific MSG Level 1.5 Native format. The Python package satpy supports reading and loading data from many input files. With the function available_readers()
, you can get a list of all available readers satpy offers. For MSG data and the Native format, you can use the reader seviri_l1b_native
.
available_readers()
['abi_l1b', 'abi_l1b_scmi', 'abi_l2_nc', 'acspo', 'agri_l1', 'ahi_hrit', 'ahi_hsd', 'ahi_l1b_gridded_bin', 'ami_l1b', 'amsr2_l1b', 'amsr2_l2', 'amsr2_l2_gaasp', 'amsub_l1c_aapp', 'ascat_l2_soilmoisture_bufr', 'avhrr_l1b_aapp', 'avhrr_l1b_eps', 'avhrr_l1b_hrpt', 'avhrr_l1c_eum_gac_fdr_nc', 'clavrx', 'cmsaf-claas2_l2_nc', 'electrol_hrit', 'fci_l1c_nc', 'fci_l2_nc', 'generic_image', 'geocat', 'glm_l2', 'goes-imager_hrit', 'goes-imager_nc', 'gpm_imerg', 'hy2_scat_l2b_h5', 'iasi_l2', 'iasi_l2_so2_bufr', 'jami_hrit', 'maia', 'mersi2_l1b', 'mhs_l1c_aapp', 'mimicTPW2_comp', 'mirs', 'modis_l1b', 'modis_l2', 'msi_safe', 'mtsat2-imager_hrit', 'mviri_l1b_fiduceo_nc', 'nucaps', 'nwcsaf-geo', 'nwcsaf-msg2013-hdf5', 'nwcsaf-pps_nc', 'olci_l1b', 'olci_l2', 'omps_edr', 'safe_sar_l2_ocn', 'sar-c_safe', 'satpy_cf_nc', 'seadas_l2', 'seviri_l1b_hrit', 'seviri_l1b_icare', 'seviri_l1b_native', 'seviri_l1b_nc', 'seviri_l2_bufr', 'seviri_l2_grib', 'slstr_l1b', 'slstr_l2', 'smos_l2_wind', 'tropomi_l2', 'vaisala_gld360', 'vii_l1b_nc', 'vii_l2_nc', 'viirs_compact', 'viirs_edr_active_fires', 'viirs_edr_flood', 'viirs_l1b', 'viirs_sdr', 'virr_l1b']
From the EUMETSAT Data Store, we downloaded for severals hours from 5./6. February 2021 of High Rate SEVIRI Level 1.5 Image data. The data is available in the folder ../../../eodata/dust/part1/1_satellite/meteosat/event/
. Let us load one image scence for 6 February at 7.27 UTC. First, we specify the file path and create a variable with the name file_name
.
file_name = glob.glob('../../../eodata/dust/part1/1_satellite/meteosat/event/*2021020607*.nat')
file_name
['../../../eodata/dust/part1/1_satellite/meteosat/event/MSG4-SEVI-MSG15-0100-NA-20210206072743.827000000Z-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 0x7f64e63b0910>
A Scene
object loads all band information of a satellite image. With the function available_dataset_names()
, you can see the available bands available from the MSG SEVIRI satellite.
scn.available_dataset_names()
['HRV', 'IR_016', 'IR_039', 'IR_087', 'IR_097', 'IR_108', 'IR_120', 'IR_134', 'VIS006', 'VIS008', 'WV_062', 'WV_073']
The underlying container for data in satpy is the xarray.DataArray
. With the function load()
, you can specify an individual band by name, e.g. IR_108
and load the data. If you then select the loaded band, you see that the band object is a xarray.DataArray
.
scn.load(['IR_108'])
scn['IR_108']
<xarray.DataArray 'reshape-e400215cfba948839ca94f9545f4740f' (y: 3712, x: 3712)> dask.array<truediv, shape=(3712, 3712), dtype=float32, chunksize=(3712, 3712), chunktype=numpy.ndarray> Coordinates: acq_time (y) datetime64[ns] NaT NaT NaT NaT NaT NaT ... NaT NaT NaT NaT NaT crs object PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["unknown",... * y (y) float64 -5.566e+06 -5.563e+06 ... 5.566e+06 5.569e+06 * x (x) float64 5.566e+06 5.563e+06 5.56e+06 ... -5.566e+06 -5.569e+06 Attributes: orbital_parameters: {'projection_longitude': 0.0, 'projection_latit... units: K wavelength: 10.8 µm (9.8-11.8 µm) standard_name: toa_brightness_temperature platform_name: Meteosat-11 sensor: seviri georef_offset_corrected: True start_time: 2021-02-06 07:15:11.201331 end_time: 2021-02-06 07:30:10.197941 reader: seviri_l1b_native area: Area ID: msg_seviri_fes_3km\nDescription: MSG S... name: IR_108 resolution: 3000.403165817 calibration: brightness_temperature modifiers: () _satpy_id: DataID(name='IR_108', wavelength=WavelengthRang... ancillary_variables: []
|
array(['NaT', 'NaT', 'NaT', ..., 'NaT', 'NaT', 'NaT'], dtype='datetime64[ns]')
array(<Projected CRS: PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["unk ...> Name: unknown Axis Info [cartesian]: - E[east]: Easting (metre) - N[north]: Northing (metre) Area of Use: - undefined Coordinate Operation: - name: unknown - method: Geostationary Satellite (Sweep Y) Datum: unknown - Ellipsoid: unknown - Prime Meridian: Greenwich , dtype=object)
array([-5565747.872591, -5562747.469425, -5559747.066259, ..., 5562747.469425, 5565747.872591, 5568748.275757])
array([ 5565747.872591, 5562747.469425, 5559747.066259, ..., -5562747.469425, -5565747.872591, -5568748.275757])
With an xarray data structure, you can handle the object as a xarray.DataArray
. For example, you can print a list of available attributes with the function attrs.keys()
.
scn['IR_108'].attrs.keys()
dict_keys(['orbital_parameters', 'units', 'wavelength', 'standard_name', 'platform_name', 'sensor', 'georef_offset_corrected', 'start_time', 'end_time', 'reader', 'area', 'name', 'resolution', 'calibration', 'modifiers', '_satpy_id', 'ancillary_variables'])
With the attrs()
function, you can also access individual metadata information, e.g. start_time
and end_time
.
scn['IR_108'].attrs['start_time'], scn['IR_108'].attrs['end_time']
(datetime.datetime(2021, 2, 6, 7, 15, 11, 201331), datetime.datetime(2021, 2, 6, 7, 30, 10, 197941))
But for now, let us reload the data file from the beginning. We do not need the IR_108
channel for the next steps, but rather would like to focus on different RGB composites
.
scn = Scene(reader="seviri_l1b_native",
filenames=file_name)
scn
<satpy.scene.Scene at 0x7f64a5788c40>
RGB composites combine multiple window channels of satellite data in order to get e.g. a true-color image of the scene. Depending on the channel combination used, different features can be highlighted in the composite, e.g. dust. SatPy offers several predefined RGB composites options. The function available_composite_ids()
returns a list of available composite IDs. You see that there are predefined composites for e.g. natural_color
, snow
or dust
.
scn.available_composite_ids()
[DataID(name='airmass'), DataID(name='ash'), DataID(name='cloudtop'), DataID(name='cloudtop_daytime'), DataID(name='colorized_ir_clouds'), DataID(name='convection'), DataID(name='day_microphysics'), DataID(name='day_microphysics_winter'), DataID(name='dust'), DataID(name='fog'), DataID(name='green_snow'), DataID(name='hrv_clouds'), DataID(name='hrv_fog'), DataID(name='hrv_severe_storms'), DataID(name='hrv_severe_storms_masked'), DataID(name='ir108_3d'), DataID(name='ir_cloud_day'), DataID(name='ir_overview'), DataID(name='ir_sandwich'), DataID(name='natural_color'), DataID(name='natural_color_nocorr'), DataID(name='natural_color_raw'), DataID(name='natural_color_with_night_ir'), DataID(name='natural_color_with_night_ir_hires'), DataID(name='natural_enh'), DataID(name='natural_enh_with_night_ir'), DataID(name='natural_enh_with_night_ir_hires'), DataID(name='natural_with_night_fog'), DataID(name='night_fog'), DataID(name='night_ir_alpha'), DataID(name='night_ir_with_background'), DataID(name='night_ir_with_background_hires'), DataID(name='night_microphysics'), DataID(name='overview'), DataID(name='overview_raw'), DataID(name='realistic_colors'), DataID(name='snow'), DataID(name='vis_sharpened_ir')]
Let us define a list with the composite ID natural_color
. This list (composite_id
) 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']
scn.load(composite_id, upper_right_corner="NE")
INFO:satpy.readers.yaml_reader:Flipping Dataset unknown_name upsidedown. INFO:satpy.readers.yaml_reader:Flipping Dataset unknown_name leftright. INFO:satpy.readers.yaml_reader:Flipping Dataset unknown_name upsidedown. INFO:satpy.readers.yaml_reader:Flipping Dataset unknown_name leftright. INFO:satpy.readers.yaml_reader:Flipping Dataset unknown_name upsidedown. INFO:satpy.readers.yaml_reader:Flipping Dataset unknown_name leftright.
A print of the Scene object scn
shows you that it consists of one xarray.DataArray
with the standard name natural_color
.
print(scn)
<xarray.DataArray 'where-05559b29beb8ff0ff25046286978f4a2' (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: sun_earth_distance_correction_applied: True sun_earth_distance_correction_factor: 0.972285159690274 standard_name: natural_color start_time: 2021-02-06 07:15:11.201331 ancillary_variables: [] area: Area ID: msg_seviri_fes_3km\nDesc... orbital_parameters: {'projection_longitude': 0.0, 'pr... georef_offset_corrected: True reader: seviri_l1b_native sensor: seviri resolution: 3000.403165817 platform_name: Meteosat-11 end_time: 2021-02-06 07:30:10.197941 wavelength: None optional_datasets: [] name: natural_color _satpy_id: DataID(name='natural_color', reso... prerequisites: [DataQuery(name='IR_016', modifie... optional_prerequisites: [] mode: RGB
The function show()
allows you to visualize a loaded composite or band. Let us e.g. visualize the natural_color
composite. Once visualized, you see that the RGB composite highlights clouds in turquoise, but highlights specific features in their natural color.
scn.show('natural_color')