{ "cells": [ { "cell_type": "markdown", "id": "c587fbf9", "metadata": {}, "source": "# SEVIRI Level 1.5" }, { "cell_type": "markdown", "id": "4e468e0f-d850-420e-9d78-872164e7e8a0", "metadata": {}, "source": [ "## Context\n", "### Purpose\n", "Explore [SEVIRI](https://www.eumetsat.int/seviri) satellite imagery and wildfire data that is open and free for scientific use.\n", "\n", "### Sensor description\n", "The [SEVIRI Level 1.5 Image Data](https://data.eumetsat.int/data/map/EO:EUM:DAT:MSG:HRSEVIRI) product contains provides 12 spectral channels. Level 1.5 image data corresponds to the geolocated and radiometrically pre-processed image data, ready for further processing, e.g. the extraction of meteorological products.\n", "\n", "### Highlights\n", "* Use `satpy` to load, visualise, and regrid SEVIRI level 1.5 data.\n", "* Fetch a fire database containing some 8080 fires from September 1st, 2020.\n", "* Visualisation of fire pixels from the database.\n", "* Visualisation of the fire pixels alongside bands from the SEVIRI satellite data.\n", "* Demonstration of how to write a custom `intake` driver for `satpy`.\n", "\n", ":::{note}\n", "The author acknowledges [EUMETSAT](https://www.eumetsat.int).\n", ":::" ] }, { "cell_type": "markdown", "id": "18738ce3-a86f-4618-a300-69c159d4cb8f", "metadata": {}, "source": [ "## Install and load libraries" ] }, { "cell_type": "code", "id": "2c1778e3-4131-4339-8afd-f730dfb1c6ae", "metadata": { "pycharm": { "is_executing": true }, "tags": [ "hide-input" ] }, "source": [ "!pip -q install pyspectral\n", "!pip -q install 'satpy==0.26.0'\n", "!pip -q install pyorbital\n", "!pip -q install geopandas\n", "!pip -q install geoviews" ], "outputs": [], "execution_count": null }, { "cell_type": "code", "id": "5c6e9348-a1db-4463-8d81-53f6b5206e29", "metadata": { "tags": [ "hide-input" ] }, "source": [ "import pandas as pd\n", "import numpy as np\n", "import geopandas\n", "import intake\n", "import fsspec, aiohttp\n", "import hvplot.xarray\n", "import hvplot.pandas\n", "import holoviews as hv\n", "import panel as pn\n", "import satpy\n", "import xarray as xr\n", "import tempfile\n", "from scipy.spatial import cKDTree\n", "from satpy.writers import get_enhanced_image\n", "from getpass import getpass\n", "from pathlib import Path\n", "from pyresample import geometry\n", "from pyresample import create_area_def\n", "import datetime\n", "import urllib.request\n", "import os.path\n", "import requests\n", "from pathlib import Path\n", "from dotenv import load_dotenv\n", "import warnings\n", "warnings.filterwarnings(action='ignore')" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "id": "47657442-aadd-49a8-a182-20d32b18a111", "metadata": {}, "source": [ "## Set project structure" ] }, { "cell_type": "code", "id": "8ceaf536-eff6-4d5a-a69e-682896bc991e", "metadata": {}, "source": [ "notebook_folder = Path('./notebook')\n", "if not notebook_folder.exists():\n", " notebook_folder.mkdir(exist_ok=True)" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "id": "7f7e01b0-e91d-4b6e-8732-eaf8322d0da9", "metadata": {}, "source": [ "## Fetch Data\n", "\n", "\n", ":::{note}\n", "To download data from the [EUMETSAT's Data site](https://data.eumetsat.int/data/map/EO:EUM:DAT:MSG:HRSEVIRI) you must have a valid account. Please register with EUMETSAT's data sevices if you do not already have an account. Then provide your consumer key and consumer secret when prompted in the cell below. Your consumer key and consumer secret can be found at the following url: https://api.eumetsat.int/api-key/\n", " \n", "Now you should successfully be able to download SEVIRI data.\n", ":::" ] }, { "cell_type": "code", "id": "c3b896d7-71e2-4534-bca0-41da8c661a36", "metadata": {}, "source": [ "file_path = notebook_folder / 'MSG4-SEVI-MSG15-0100-NA-20200901104243.230000000Z-NA.nat'\n", "\n", "if not file_path.is_file() or file_path.stat().st_size == 0:\n", " load_dotenv()\n", " consumer_key = os.environ.get(\"EUMESAT_API_KEY\") #replace for your EUMESAT consumer key if run local or in Binder\n", " consumer_secret = os.environ.get(\"EUMESAT_API_SECRET\") #replace for your EUMESAT consumer secret if run local or in Binder\n", "\n", " token_url = 'https://api.eumetsat.int/token'\n", " response = requests.post(\n", " token_url,\n", " auth=requests.auth.HTTPBasicAuth(consumer_key, consumer_secret),\n", " data = {'grant_type': 'client_credentials'},\n", " headers = {\"Content-Type\" : \"application/x-www-form-urlencoded\"}\n", " )\n", "\n", " access_token = response.json()['access_token']\n", "\n", " filename = 'MSG4-SEVI-MSG15-0100-NA-20200901104243.230000000Z-NA.nat'\n", "\n", " product_url = 'https://api.eumetsat.int/data/download/collections/EO%3AEUM%3ADAT%3AMSG%3AHRSEVIRI/products/MSG4-SEVI-MSG15-0100-NA-20200901104243.230000000Z-NA/entry'\n", " product_url += f'?name={filename}'\n", " product_url += f'&access_token={access_token}'\n", "\n", " urllib.request.urlretrieve(product_url, str(file_path))" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "id": "66331370-493e-426a-b730-245913150d8f", "metadata": {}, "source": [ "Download the fire pixel data for this day from Zenodo. This data is not directly downloadable from the internet, so we share a subset of fires for this imagery here." ] }, { "cell_type": "code", "id": "fc8b87ad-f365-45b0-8ef4-5c052420c5c6", "metadata": {}, "source": [ "filename = 'HDF5_LSASAF_MSG_FRP-PIXEL-ListProduct_MSG-Disk_202009011045'\n", "file_path = notebook_folder / filename\n", "url = f'https://zenodo.org/record/5717106/files/{filename}?download=1'\n", "\n", "if not file_path.is_file() or file_path.stat().st_size == 0:\n", " urllib.request.urlretrieve(url, file_path)" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "id": "7b2e4579-288e-4ae9-b9e3-e66f17d5ffa2", "metadata": {}, "source": [ "Load an intake catalog for the downloaded data" ] }, { "cell_type": "code", "id": "30d99b41-a3d4-4c0c-af52-292be4a0a2a1", "metadata": {}, "source": [ "catalog_file = notebook_folder / 'catalog.yaml'\n", "\n", "with catalog_file.open('w') as f:\n", " f.write('''\n", "sources:\n", " seviri_l1b:\n", " args:\n", " urlpath: 'notebook/MSG4-SEVI-MSG15-0100-NA-20200901104243.230000000Z-NA.nat'\n", " reader: 'seviri_l1b_native'\n", " description: \"SEVIRI Level 1.5 Products\"\n", " driver: SatpySource\n", " metadata: {}\n", " seviri_fires:\n", " args:\n", " urlpath: 'notebook/HDF5_LSASAF_MSG_FRP-PIXEL-ListProduct_MSG-Disk_202009011045'\n", " description: \"SEVIRI Level 2 Fires\"\n", " driver: netcdf\n", " metadata: {}\n", "''')" ], "outputs": [], "execution_count": null }, { "cell_type": "code", "id": "b4ce368f-5bdc-4200-9119-00943fdfd83e", "metadata": {}, "source": [ "from intake.source.base import PatternMixin\n", "from intake.source.base import DataSource, Schema\n", "\n", "import glob\n", "\n", "class SatpySource(DataSource, PatternMixin):\n", " \"\"\"Intake driver for data supported by satpy.Scene\"\"\"\n", " \n", " container = 'python'\n", " name = 'satpy'\n", " version = '0.0.1'\n", " partition_access = False\n", "\n", " def __init__(self, urlpath, reader=None, metadata=None, path_as_pattern=True):\n", " \"\"\"\n", " Parameters\n", " ----------\n", " path: str, location of model pkl file\n", " Either the absolute or relative path to the file\n", " opened. Some examples:\n", " - ``{{ CATALOG_DIR }}/file.nat``\n", " - ``{{ CATALOG_DIR }}/*.nc``\n", " reader: str, optional\n", " Name of the satpy reader to use when loading data (ie. ``modis_l1b``)\n", " metadata: dict, optional\n", " Additional metadata to pass to the intake catalog\n", " path_as_pattern: bool or str, optional\n", " Whether to treat the path as a pattern (ie. ``data_{field}.nc``)\n", " and create new coodinates in the output corresponding to pattern\n", " fields. If str, is treated as pattern to match on. Default is True.\n", " \"\"\"\n", "\n", " self.path_as_pattern = path_as_pattern\n", " self.urlpath = urlpath\n", " self._reader = reader\n", " super(SatpySource, self).__init__(metadata=metadata)\n", "\n", " def _load(self):\n", " files = self.urlpath\n", " files = list(glob.glob(files))\n", " return satpy.Scene(files, reader=self._reader)\n", " \n", " def _get_schema(self):\n", " self._schema = Schema(\n", " npartitions=1,\n", " extra_metadata={}\n", " )\n", " return self._schema\n", "\n", " def read(self):\n", " self._load_metadata()\n", " return self._load()\n", "\n", "intake.register_driver('SatpySource', SatpySource)\n", "cat = intake.open_catalog(catalog_file)" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "id": "b1f4605d-8ae1-4390-afc4-7f79a5e575a5", "metadata": {}, "source": [ "## Load SEVIRI Satellite Data\n", "\n", "Here we use `satpy` to load the SEVIRI level 1b data into memory. As well as loading the image data, `satpy` provides a easy way to regrid the data to different coordinate systems." ] }, { "cell_type": "code", "id": "6863149d-1cd5-4713-8057-0bc6a50db1fc", "metadata": {}, "source": [ "scn = cat['seviri_l1b'].read()\n", "scn" ], "outputs": [], "execution_count": null }, { "cell_type": "code", "id": "ac05742d-6d30-4b8f-9e22-c578d61339f8", "metadata": { "tags": [ "hide-input" ] }, "source": [ "scn.load(['natural_color', 'IR_039', 'IR_108'])\n", "plot_seviri_raw = scn.show('natural_color')\n", "plot_seviri_raw" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "id": "0b782471", "metadata": {}, "source": [ "### Resample to a different projection\n", "\n", "In the plot above you can see that the raw SEVIRI data has distortion towards edge of the image. By regridding the data we can avoid some of this distortion." ] }, { "cell_type": "code", "id": "bfa5c2e1", "metadata": {}, "source": [ "area_id = \"Southern Africa\"\n", "description = \"Southern Africa in Mercator-Projection\"\n", "proj_id = \"Southern Africa\"\n", "proj_dict = {\"proj\": \"eqc\"}\n", "\n", "width = 1000 # width of the result domain in pixels\n", "height = 1000 # height of the result domain in pixels\n", "\n", "llx = 0 # projection x coordinate of lower left corner of lower left pixel\n", "lly = -30e5 # projection y coordinate of lower left corner of lower left pixel\n", "urx = 55e5 # projection x coordinate of upper right corner of upper right pixel\n", "ury = 10e5 # projection y coordinate of upper right corner of upper right pixel\n", "\n", "area_extent = (llx,lly,urx,ury)\n", "\n", "resolution=3000\n", "center =(0,0)\n", "area_def = create_area_def(area_id, proj_dict, resolution=resolution, area_extent=area_extent)\n", "\n", "seviri_scn = scn.resample(area_def, radius_of_influence=10000)" ], "outputs": [], "execution_count": null }, { "cell_type": "code", "id": "ee8a3dc0", "metadata": { "tags": [ "hide-input" ] }, "source": [ "plot_seviri_scn = seviri_scn.show('natural_color')\n", "plot_seviri_scn" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "id": "efe2162d", "metadata": {}, "source": [ "## Load SEVIRI Fire Database\n", "\n", "Now we're going to load the SEVIRI fire database from HDF file. This contains the longitude, latitude, and time of where fires have been detected to occur. It also includes an estimate of the Fire Radiative Power (FRP) {cite:p}`wooster2015lsa, wooster_he_xu_lattanzio`, a measure of the intensity of the fire, for each fire detected." ] }, { "cell_type": "code", "id": "fd402579", "metadata": { "tags": [ "hide-input" ] }, "source": [ "# Scale factors for the SEVIRI fire product from:\n", "# https://nextcloud.lsasvcs.ipma.pt/s/CZa8RwDxjGqYezS?dir=undefined&openfile=105793\n", "\n", "SCALE_FACTORS = dict(\n", " FRP=10,\n", " FRP_UNCERTAINTY=100,\n", " ERR_FRP_COEFF=10000,\n", " ERR_BACKGROUND=10000,\n", " ERR_ATM_TRANS=10000,\n", " ERR_VERT_COMP=10000,\n", " ERR_RADIOMETRIC=10000,\n", " LATITUDE=100,\n", " LONGITUDE=100,\n", " FIRE_CONFIDENCE=100,\n", " BT_MIR=10,\n", " BT_TIR=10,\n", " BW_BT_MIR=10,\n", " BW_BTD=10,\n", " PIXEL_SIZE=100,\n", " PIXEL_VZA=100,\n", " PIXEL_ATM_TRANS=10000,\n", " RAD_PIX=10000,\n", " STD_BCK=10000\n", ")\n", "\n", "def process_fire_product(fire_pixels):\n", " # Hack: xarray randomly loads some columns as coordinates, which don't get converted to a dataframe\n", " # Correct here by removing them as coords and renaming them back to their original name\n", " coords = list(fire_pixels.coords.keys())\n", " fire_pixels = fire_pixels.reset_index(coords).reset_coords()\n", " fire_pixels = fire_pixels.swap_dims({key: f'phony_dim_0' for key in list(fire_pixels.dims.keys())})\n", " fire_pixels = fire_pixels.rename({f\"{name}_\": name for name in coords})\n", " fire_pixels = fire_pixels.to_dataframe()\n", "\n", " for c in fire_pixels.columns:\n", " if c in SCALE_FACTORS:\n", " fire_pixels[c] = fire_pixels[c] / SCALE_FACTORS[c]\n", " \n", " fire_pixels['ABS_LINE_1KM'] = fire_pixels.ABS_LINE // 2\n", " fire_pixels['ABS_PIXEL_1KM'] = fire_pixels.ABS_PIXEL // 2\n", " fire_pixels.index.name = 'index'\n", " return fire_pixels\n", "\n", "def parse_l2_timestamp(product_name):\n", " \"\"\"Parse the timestamp from the SEVIRI L2 product name\"\"\"\n", " timestamp = product_name.split('_')[-1]\n", " timestamp = pd.to_datetime(timestamp, format='%Y%m%d%H%M')\n", " return timestamp\n", "\n", "# Read in fires, scale and rename dimensions\n", "fires = cat['seviri_fires'].read()\n", "fires = process_fire_product(fires)\n", "fires = fires.rename({'LONGITUDE': 'longitude', 'LATITUDE': 'latitude', 'FRP': 'frp'}, axis=1)\n", "\n", "# Grab the timestamp of the product\n", "urlpath = cat['seviri_fires'].describe()['args']['urlpath']\n", "fires['timestamp'] = parse_l2_timestamp(urlpath)\n", "\n", "# Convert to geopandas\n", "fires = geopandas.GeoDataFrame(\n", " fires, geometry=geopandas.points_from_xy(fires.longitude, fires.latitude))\n", "\n", "# We're only interested in data from Southern Africa for now\n", "llx = 0 # projection x coordinate of lower left corner of lower left pixel\n", "lly = -30 # projection y coordinate of lower left corner of lower left pixel\n", "urx = 55 # projection x coordinate of upper right corner of upper right pixel\n", "ury = 10 # projection y coordinate of upper right corner of upper right pixel\n", "\n", "fires = fires.cx[llx:urx, lly:ury]\n", "fires = fires.sort_index()\n", "fires" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "id": "17cbf1b2", "metadata": {}, "source": [ "### Geographical distribution of Fires\n", "\n", "Here we plot the geographical distribution of fires detected by SEVIRI over Southern Africa." ] }, { "cell_type": "code", "id": "699458dc", "metadata": { "tags": [ "hide-input" ] }, "source": [ "fires.hvplot.points('longitude', 'latitude', c='frp', geo=True, alpha=0.2,\n", " tiles='ESRI', xlim=(llx, urx), ylim=(lly, ury), cmap='autumn', logz=True,\n", " dynamic=False)" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "id": "1be8f654", "metadata": {}, "source": [ "## Visualising Fire Pixels with Satellite Imagery\n", "\n", "Visualise a color image of the SEVIRI region using `hvplot`." ] }, { "cell_type": "code", "id": "071602fc", "metadata": { "tags": [ "hide-input" ] }, "source": [ "seviri_dataset = seviri_scn.to_xarray_dataset()\n", "\n", "img = get_enhanced_image(seviri_scn['natural_color'].squeeze())\n", "img = img.data\n", "img = img.clip(0, 1)\n", "img = (img * 255).astype(np.uint8)\n", "\n", "seviri_dataset['natural_color'] = img\n", "\n", "grid = seviri_scn.min_area().get_lonlats()\n", "lons, lats = grid[0][0], grid[1][:, 0]\n", "seviri_dataset = seviri_dataset.assign_coords(dict(x=lons, y=lats))\n", "seviri_dataset = seviri_dataset.rename(dict(x='longitude', y='latitude'))\n", "\n", "plot_SEVIRI_rgb = seviri_dataset['natural_color'].hvplot.rgb(x='longitude', y='latitude', bands='bands', rasterize=True, data_aspect=1)\n", "plot_SEVIRI_rgb" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "id": "70d298be", "metadata": {}, "source": [ "Now overlay the fire pixels ontop of the SEVIRI image, along with the FRP for each fire pixel. Try zooming in with `rasterize=False`, you should be able to see clear smoke trails at the locations of some of the fires!" ] }, { "cell_type": "code", "id": "c7db11ba", "metadata": { "tags": [ "hide-input" ] }, "source": [ "rgb = seviri_dataset['natural_color'].hvplot.rgb(x='longitude', y='latitude', bands='bands', data_aspect=1, hover=False, title='True Colour', rasterize=True)\n", "points = fires.hvplot.points('longitude', 'latitude', c='frp', cmap='autumn', logz=True, alpha=0.4)\n", "plot_fires_SEVIRI_rgb = rgb*points\n", "plot_fires_SEVIRI_rgb" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "id": "0497cc2d", "metadata": {}, "source": [ "We can also overlay the fire pixels ontop of the band difference between SEVIRI 3.9 and 10.8 micron wavelength bands. The 3.9 and 10.8 bands are thermal bands. Fires will appear as very bright pixels in the difference between these two bands. Try zooming in with `rasterize=False`, you should be able to clearly see bright spots at the location of each fire pixel." ] }, { "cell_type": "code", "id": "be16d973", "metadata": { "tags": [ "hide-input" ] }, "source": [ "band_difference = seviri_dataset['IR_039'] - seviri_dataset['IR_108']\n", "thermal = band_difference.hvplot.image(x='longitude', y='latitude', cmap='viridis', data_aspect=1, hover=False, title='Band 20: 3.7 micron', rasterize=True)\n", "points = fires.hvplot.points('longitude', 'latitude', c='frp', cmap='autumn', logz=True, alpha=0.4)\n", "plot_fires_SEVIRI_thermal = thermal*points\n", "plot_fires_SEVIRI_thermal" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "id": "8f0dc221", "metadata": {}, "source": [ "## Summary\n", "\n", "This notebook has demonstrated the use of:\n", " - The `satpy` package to easily load, plot and regrid satellite data from the SEVIRI sensor.\n", " - `hvplot` to interatively visualise wildfire pixels detected from the SEVIRI sensor.\n", " - `geopandas` to load, filter, and manipulate historical wildfire pixel data.\n", " - `intake` the load data, including using a custom plugin for obsecure file formats." ] }, { "cell_type": "markdown", "id": "46443ea9", "metadata": {}, "source": [ "## Citing this Notebook" ] }, { "cell_type": "markdown", "id": "7435c139", "metadata": {}, "source": "Please see [CITATION.cff](https://github.com/eds-book/bc30df18-fce2-42fa-aade-1ce5b7f3ca3c/blob/main/CITATION.cff) for the full citation information. The citation file can be exported to APA or BibTex formats (learn more [here](https://docs.github.com/en/repositories/managing-your-repositorys-settings-and-features/customizing-your-repository/about-citation-files))." }, { "cell_type": "markdown", "id": "3b2dd27b", "metadata": {}, "source": [ "## Additional information\n", "\n", "**Review**: This notebook has been reviewed by one or more members of the Environmental Data Science book community. The open review is available [here](https://github.com/alan-turing-institute/environmental-ds-book/pull/12).\n", "\n", "**Dataset originator/creator**:\n", "* SEVIRI Level 1.5 Image Data - MSG - 0 degree, European Organisation for the Exploitation of Meteorological Satellites (EUMETSAT)\n", "* FRPPIXEL, Land Surface Analysis, Satellite Application Facility on Land Surface Analysis (LSA SAF)\n", "\n", "**License**: The code in this notebook is licensed under the MIT License. The Environmental Data Science book is licensed under the Creative Commons by Attribution 4.0 license. See further details [here](https://github.com/alan-turing-institute/environmental-ds-book/blob/main/LICENSE).\n", "\n", "**Contact**: If you have any suggestion or report an issue with this notebook, feel free to [create an issue](https://github.com/alan-turing-institute/environmental-ds-book/issues/new/choose) or send a direct message to [environmental.ds.book@gmail.com](mailto:environmental.ds.book@gmail.com)." ] }, { "cell_type": "code", "id": "67fe9e8b", "metadata": { "tags": [ "remove-input" ] }, "source": [ "from datetime import date\n", "\n", "print('Notebook repository version: v2025.6.0')\n", "print(f'Last tested: {date.today()}')" ], "outputs": [], "execution_count": null } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.10" } }, "nbformat": 4, "nbformat_minor": 5 }