{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# MSG SEVIRI RGB composites" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "```{hint} \n", "Execute the notebook on the training platform >>\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This notebook provides you an introduction to Level 1.5 data from the SEVIRI instrument as part of the Meteosat Second Generation 0 Degree Service. It further provides you 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.\n", "\n", "The example features the saharan dust event during 8 April 2024, which brought massive amounts of Saharan dust into central Europe. See a more detailed description of this particular event here.\n", "\n", "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.\n", "\n", "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.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "```{admonition} Basic facts\n", "**Spatial resolution**: `1 km at nadir`
\n", "**Spatial coverage**: `Latitude: -81 to 81 degrees`, `Longitude: -79 to 79 degrees`
\n", "**Revisit time**: `every 15 minutes`
\n", "**Data availability**: `since 2004`\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "```{admonition} How to access the data\n", "High Rate SEVIRI Level 1.5 Image Data is available for download via the EUMETSAT Data Store. For the available time range you can select `2024-04-08` between `10:00 and 16:00 UTC`. Then click on *Show Results* to be able to select and download the data files needed.\n", "\n", "**Note:** You need to register for the EUMETSAT Earth Observation Portal in order to be able to download data from the Data Store.\n", "```\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Load required libraries**" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "import glob\n", "import xarray as xr\n", "import numpy as np\n", "\n", "import matplotlib.pyplot as plt\n", "import matplotlib.colors\n", "from matplotlib.axes import Axes\n", "\n", "import cartopy\n", "import cartopy.crs as ccrs\n", "from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER\n", "from cartopy.mpl.geoaxes import GeoAxes\n", "GeoAxes._pcolormesh_patched = Axes.pcolormesh\n", "\n", "from satpy.scene import Scene\n", "from satpy.composites import GenericCompositor\n", "from satpy.writers import to_image\n", "from satpy.resample import get_area_def\n", "from satpy import available_readers\n", "\n", "import pyresample\n", "\n", "import warnings\n", "warnings.filterwarnings('ignore')\n", "warnings.simplefilter(action = \"ignore\", category = RuntimeWarning)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Load helper functions**" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "%run ../../functions.ipynb" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## MSG SEVIRI True color composite" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Load and browse High Rate SEVIRI Level 1.5 image data" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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`." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "available_readers()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "From the EUMETSAT Data Store, we downloaded for several hours from 8 April 2024 of High Rate SEVIRI Level 1.5 Image data. The data is available in the folder `../../eodata/1_satellite/meteosat/`. Let us load one image scence for 8 April 2024 at 12.27 UTC. First, we specify the file path and create a variable with the name `file_name`." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "['../../eodata/1_satellite/meteosat/MSG3-SEVI-MSG15-0100-NA-20240408122742.149000000Z-NA.nat']" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "file_name = glob.glob('../../eodata/1_satellite/meteosat/*2024040812*.nat')\n", "file_name" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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.\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "scn = Scene(reader=\"seviri_l1b_native\", \n", " filenames=file_name)\n", "scn" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "['HRV',\n", " 'IR_016',\n", " 'IR_039',\n", " 'IR_087',\n", " 'IR_097',\n", " 'IR_108',\n", " 'IR_120',\n", " 'IR_134',\n", " 'VIS006',\n", " 'VIS008',\n", " 'WV_062',\n", " 'WV_073']" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "scn.available_dataset_names()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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`." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray 'reshape-a292807d186e5996fa2d94038f482f80' (y: 3712, x: 3712)> Size: 55MB\n",
       "dask.array<truediv, shape=(3712, 3712), dtype=float32, chunksize=(928, 3712), chunktype=numpy.ndarray>\n",
       "Coordinates:\n",
       "    acq_time  (y) datetime64[ns] 30kB NaT NaT NaT NaT NaT ... NaT NaT NaT NaT\n",
       "    crs       object 8B PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["unknow...\n",
       "  * y         (y) float64 30kB -5.566e+06 -5.563e+06 ... 5.566e+06 5.569e+06\n",
       "  * x         (x) float64 30kB 5.566e+06 5.563e+06 ... -5.566e+06 -5.569e+06\n",
       "Attributes: (12/18)\n",
       "    orbital_parameters:       {'projection_longitude': 0.0, 'projection_latit...\n",
       "    units:                    K\n",
       "    wavelength:               10.8 µm (9.8-11.8 µm)\n",
       "    standard_name:            toa_brightness_temperature\n",
       "    platform_name:            Meteosat-10\n",
       "    sensor:                   seviri\n",
       "    ...                       ...\n",
       "    name:                     IR_108\n",
       "    resolution:               3000.403165817\n",
       "    calibration:              brightness_temperature\n",
       "    modifiers:                ()\n",
       "    _satpy_id:                DataID(name='IR_108', wavelength=WavelengthRang...\n",
       "    ancillary_variables:      []
" ], "text/plain": [ " Size: 55MB\n", "dask.array\n", "Coordinates:\n", " acq_time (y) datetime64[ns] 30kB NaT NaT NaT NaT NaT ... NaT NaT NaT NaT\n", " crs object 8B PROJCRS[\"unknown\",BASEGEOGCRS[\"unknown\",DATUM[\"unknow...\n", " * y (y) float64 30kB -5.566e+06 -5.563e+06 ... 5.566e+06 5.569e+06\n", " * x (x) float64 30kB 5.566e+06 5.563e+06 ... -5.566e+06 -5.569e+06\n", "Attributes: (12/18)\n", " orbital_parameters: {'projection_longitude': 0.0, 'projection_latit...\n", " units: K\n", " wavelength: 10.8 µm (9.8-11.8 µm)\n", " standard_name: toa_brightness_temperature\n", " platform_name: Meteosat-10\n", " sensor: seviri\n", " ... ...\n", " name: IR_108\n", " resolution: 3000.403165817\n", " calibration: brightness_temperature\n", " modifiers: ()\n", " _satpy_id: DataID(name='IR_108', wavelength=WavelengthRang...\n", " ancillary_variables: []" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "scn.load(['IR_108'])\n", "scn['IR_108']" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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()`. " ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "dict_keys(['orbital_parameters', 'units', 'wavelength', 'standard_name', 'platform_name', 'sensor', 'georef_offset_corrected', 'time_parameters', 'start_time', 'end_time', 'reader', 'area', 'name', 'resolution', 'calibration', 'modifiers', '_satpy_id', 'ancillary_variables'])" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "scn['IR_108'].attrs.keys()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "With the `attrs()` function, you can also access individual metadata information, e.g. `start_time` and `end_time`." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(datetime.datetime(2024, 4, 8, 12, 15), datetime.datetime(2024, 4, 8, 12, 30))" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "scn['IR_108'].attrs['start_time'], scn['IR_108'].attrs['end_time']" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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`." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "scn = Scene(reader=\"seviri_l1b_native\", \n", " filenames=file_name)\n", "scn" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Browse and visualize RGB composite IDs" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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`." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[DataID(name='24h_microphysics'),\n", " DataID(name='airmass'),\n", " DataID(name='ash'),\n", " DataID(name='cloud_phase_distinction'),\n", " DataID(name='cloud_phase_distinction_raw'),\n", " DataID(name='cloudtop'),\n", " DataID(name='cloudtop_daytime'),\n", " DataID(name='colorized_ir_clouds'),\n", " DataID(name='convection'),\n", " DataID(name='day_microphysics'),\n", " DataID(name='day_microphysics_winter'),\n", " DataID(name='dust'),\n", " DataID(name='fog'),\n", " DataID(name='green_snow'),\n", " DataID(name='hrv_clouds'),\n", " DataID(name='hrv_fog'),\n", " DataID(name='hrv_severe_storms'),\n", " DataID(name='hrv_severe_storms_masked'),\n", " DataID(name='ir108_3d'),\n", " DataID(name='ir_cloud_day'),\n", " DataID(name='ir_overview'),\n", " DataID(name='ir_sandwich'),\n", " DataID(name='natural_color'),\n", " DataID(name='natural_color_nocorr'),\n", " DataID(name='natural_color_raw'),\n", " DataID(name='natural_color_raw_with_night_ir'),\n", " DataID(name='natural_color_with_night_ir'),\n", " DataID(name='natural_color_with_night_ir_hires'),\n", " DataID(name='natural_enh'),\n", " DataID(name='natural_enh_with_night_ir'),\n", " DataID(name='natural_enh_with_night_ir_hires'),\n", " DataID(name='natural_with_night_fog'),\n", " DataID(name='night_fog'),\n", " DataID(name='night_ir_alpha'),\n", " DataID(name='night_ir_with_background'),\n", " DataID(name='night_ir_with_background_hires'),\n", " DataID(name='night_microphysics'),\n", " DataID(name='overview'),\n", " DataID(name='overview_raw'),\n", " DataID(name='realistic_colors'),\n", " DataID(name='rocket_plume_day'),\n", " DataID(name='rocket_plume_night'),\n", " DataID(name='snow'),\n", " DataID(name='vis_sharpened_ir')]" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "scn.available_composite_ids()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "composite_id = ['natural_color']\n", "scn.load(composite_id, upper_right_corner=\"NE\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A print of the Scene object `scn` shows you that it consists of one `xarray.DataArray` with the standard name `natural_color`." ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " Size: 165MB\n", "dask.array\n", "Coordinates:\n", " crs object 8B PROJCRS[\"unknown\",BASEGEOGCRS[\"unknown\",DATUM[\"unknown...\n", " * y (y) float64 30kB 5.569e+06 5.566e+06 ... -5.563e+06 -5.566e+06\n", " * x (x) float64 30kB -5.569e+06 -5.566e+06 ... 5.563e+06 5.566e+06\n", " * bands (bands) " ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "scn_cropped = scn.crop(ll_bbox=(-5, 31, 20, 51))\n", "scn_cropped.show(\"natural_color\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "From the cropped image, you can load the `area` information that is stored as an attribute entry. The `area` attribute holds information related to the `projection`, `number of columns / rows` and `area extent`. \n", "\n", "**Note**: the `area extent` key provides x and y values in meters from the nadir, which is the original projection unit." ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", " \n", "\n", "\n", "
Area ID: msg_seviri_fes_3km\n",
       "Description: MSG SEVIRI Full Earth Scanning service area definition with 3 km resolution\n",
       "Projection: {'a': '6378169', 'h': '35785831', 'lon_0': '0', 'no_defs': 'None', 'proj': 'geos', 'rf': '295.488065897014', 'type': 'crs', 'units': 'm', 'x_0': '0', 'y_0': '0'}\n",
       "Number of columns: 750\n",
       "Number of rows: 482\n",
       "Area extent: (-454561.0796, 3151923.5257, 1795741.2947, 4598117.8516)