{ "cells": [ { "cell_type": "markdown", "id": "5fa9b6de", "metadata": {}, "source": "# MODIS MOD021KM and FIRMS" }, { "cell_type": "markdown", "id": "4e468e0f-d850-420e-9d78-872164e7e8a0", "metadata": {}, "source": [ "## Context\n", "### Purpose\n", "Explore [MODIS](https://cmr.earthdata.nasa.gov/search/concepts/C1378227407-LAADS.html) satellite imagery and wildfire data that is open and free for scientific use.\n", "\n", "### Sensor description\n", "The [MOD021KM](https://cmr.earthdata.nasa.gov/search/concepts/C1378227407-LAADS.html) {cite:p}`firms` product contains calibrated and geolocated at-aperture radiances for 36 discrete bands located in the 0.4 to 14.4 micron region of the electromagnetic spectrum.\n", "\n", "### Highlights\n", "* Use `satpy` to load, visualise, and regrid MODIS level 1B data.\n", "* Fetch a fire database containing some 497364 fires from 2020.\n", "* Visualisation of fire pixels from the database.\n", "* Visualisation of the fire pixels alongside bands from the MODIS satellite data.\n", "\n", ":::{note}\n", "The author acknowledges [MODIS Science Team](https://doi.org/10.5067/modis/mod021km.061) and the use of data and/or imagery from NASA's Fire Information for Resource Management System (FIRMS) (https://earthdata.nasa.gov/firms), part of NASA's Earth Observing System Data and Information System (EOSDIS).\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 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 pathlib import Path\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", "import datetime\n", "import urllib.request\n", "import os.path\n", "from dotenv import load_dotenv\n", "import warnings\n", "warnings.filterwarnings(action='ignore')" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "id": "8d7f7241", "metadata": {}, "source": [ "## Set project structure" ] }, { "cell_type": "code", "id": "88741dea", "metadata": {}, "source": [ "notebook_folder = './notebook'\n", "if not os.path.exists(notebook_folder):\n", " os.makedirs(notebook_folder)" ], "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 [NASA's Earth Data site](https://urs.earthdata.nasa.gov/) you must have a valid Earth Data account. Please register with Earth Data is you do not already have an account and then provide your username and password when prompted in the cell below.\n", ":::\n", "\n", ":::{important}\n", "In order to download MODIS level1 data you must have authorized access to `LAADS Web` on your account. To do this navigate to your earthdata profile page\n", "\n", " - Click \"Applications\"\n", " - Click \"Authorized Apps\"\n", " - Click \"Autorize More Applications\"\n", " - Search for \"LAADS Web\"\n", " - Click \"Authorize\"\n", " \n", "Now you should successfully be able to download MODIS data.\n", ":::" ] }, { "cell_type": "code", "id": "c8972bbc-980d-45e2-a19f-b46f37268c0d", "metadata": {}, "source": [ "fname = 'MOD021KM.A2020245.0840.061.2020245193036.hdf'\n", "\n", "if not os.path.isfile(os.path.join(notebook_folder, fname)) or os.path.getsize(os.path.join(notebook_folder, fname)) == 0:\n", " load_dotenv()\n", " username = os.environ.get(\"NASA_EARTHDATA_USER\") #replace for your EarthData username if run local or in Binder\n", " password = os.environ.get(\"NASA_EARTHDATA_PWD\") #replace for your EarthData password if run local or in Binder\n", "\n", " fsspec.config.conf['https'] = dict(client_kwargs={'auth': aiohttp.BasicAuth(username, password)})\n", "\n", " url = f'https://ladsweb.modaps.eosdis.nasa.gov/archive/allData/61/MOD021KM/2020/245/{fname}'\n", " filename = url.split('/')[-1]\n", " with fsspec.open(url) as f:\n", " with Path(os.path.join(notebook_folder, filename)).open('wb') as handle:\n", " data = f.read()\n", " try:\n", " data.decode('utf-8')\n", " raise RuntimeError('Could not download MODIS data! Have you authorized LAADS Web in your Eathdata account above?')\n", " except UnicodeDecodeError:\n", " handle.write(data)" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "id": "bfc54c32-0e44-4d49-97d3-b3aea7d373c3", "metadata": {}, "source": [ "Download the database of MODIS wildfires from a publically accessible location." ] }, { "cell_type": "code", "id": "1ac532b6-1ea5-4b0e-8159-acbb75855fcd", "metadata": { "tags": [ "scroll-output" ] }, "source": [ "filepath = 'https://firms2.modaps.eosdis.nasa.gov/data/country/zips'\n", "filename = 'modis_2020_all_countries.zip'\n", "\n", "if not os.path.isfile(filename) or os.path.getsize(filename) == 0:\n", " urllib.request.urlretrieve(os.path.join(filepath, filename), os.path.join(notebook_folder, filename))\n", " !unzip -o ./notebook/modis_2020_all_countries.zip -d ./notebook" ], "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": [ "# set catalogue location\n", "catalog_file = os.path.join(notebook_folder, 'catalog.yaml')\n", "\n", "with open(catalog_file, 'w') as f:\n", " f.write('''\n", "sources:\n", " modis_l1b:\n", " args:\n", " urlpath: 'notebook/MOD021KM.A2020245.0840.061.2020245193036.hdf'\n", " reader: 'modis_l1b'\n", " description: \"MODIS Level 1B Products\"\n", " driver: SatpySource\n", " metadata: {}\n", " modis_fires:\n", " args:\n", " urlpath: 'notebook/modis/2020/modis_*.csv'\n", " description: \"MODIS Level 2 Fires\"\n", " driver: csv\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 MODIS Satellite Data\n", "\n", "Here we use `satpy` to load the MODIS 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['modis_l1b'].read()\n", "scn" ], "outputs": [], "execution_count": null }, { "cell_type": "code", "id": "ac05742d-6d30-4b8f-9e22-c578d61339f8", "metadata": { "tags": [ "hide-input", "scroll-output" ] }, "source": [ "scn.load(['true_color', '20'], resolution=1000)\n", "plot_modis_raw = scn.show('true_color')\n", "plot_modis_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 MODIS 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 = 23E5 # projection x coordinate of lower left corner of lower left pixel\n", "lly = -22.9E5 # projection y coordinate of lower left corner of lower left pixel\n", "urx = 48E5 # projection x coordinate of upper right corner of upper right pixel\n", "ury = -1.9E5 # projection y coordinate of upper right corner of upper right pixel\n", "\n", "area_extent = (llx,lly,urx,ury)\n", "from pyresample import create_area_def\n", "resolution=1000\n", "center =(0,0)\n", "area_def = create_area_def(area_id, proj_dict, center=center, resolution=resolution)\n", "\n", "modis_scn = scn.resample(area_def, radius_of_influence=10000)" ], "outputs": [], "execution_count": null }, { "cell_type": "code", "id": "ee8a3dc0", "metadata": { "tags": [ "hide-input" ] }, "source": [ "plot_modis_scn = modis_scn.show('true_color')\n", "plot_modis_scn" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "id": "efe2162d", "metadata": {}, "source": [ "## Load MODIS Fire Database\n", "\n", "Now we're going to load the modis fire database from CSV files. 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), a measure of the intensity of the fire, for each fire detected." ] }, { "cell_type": "code", "id": "fd402579", "metadata": { "tags": [ "hide-input" ] }, "source": [ "fires = cat['modis_fires'].read()\n", "\n", "time = fires.acq_date.astype(str) + ' ' + fires.acq_time.astype(str).str.zfill(4)\n", "fires['timestamp'] = pd.to_datetime(time, format='%Y-%m-%d %H%M')\n", "fires = fires.loc[fires.satellite == 'Terra']\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.set_index('timestamp', drop=False, inplace=True)\n", "fires = fires.sort_index()\n", "fires = fires.loc['2020-09-01 00:00:00':'2020-10-01 00:00:00']\n", "fires" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "id": "1abfcc6e", "metadata": {}, "source": [ "### Daily Fire Radiative Power\n", "\n", "Here we plot the Fire Radiative Power (FRP) for each day of the month of September. FRP is a measure of the intensity of a fire in units of MegaWatts." ] }, { "cell_type": "code", "id": "856d45fd", "metadata": { "tags": [ "hide-input" ] }, "source": [ "from bokeh.models.formatters import DatetimeTickFormatter\n", "\n", "def plot_timeseries(data, y_variable):\n", " \"\"\"Timeseries plot showing the mean fire radiative power at each day as well as the 5% and 95%\"\"\"\n", "\n", " def percentile(n):\n", " def percentile_(x):\n", " return np.percentile(x, n)\n", " percentile_.__name__ = 'percentile_%s' % n\n", " return percentile_\n", "\n", " formatter = DatetimeTickFormatter(months='%b')\n", "\n", " bounds = data.groupby([data.index.day_of_year])[y_variable].agg([percentile(5), percentile(95)])\n", " avg = data.groupby([data.index.day_of_year])[y_variable].mean()\n", "\n", " bounds.index = data.groupby([data.index.day_of_year])['timestamp'].first()\n", " avg.index = bounds.index\n", "\n", " tseries = hv.Overlay([\n", " (bounds.hvplot.area('timestamp', 'percentile_5', 'percentile_95',\n", " alpha=0.2, color='red', xformatter=formatter)),\n", " avg.hvplot.line(x='timestamp', label=f'Average Daily FRP', color='red', xformatter=formatter)])\n", "\n", " return tseries.options(width=800, height=400, xlabel='Day',ylabel=f'FRP (MW)',legend_position='top_left')\n", "\n", "plot_FIRMS_ts = plot_timeseries(fires, 'frp')\n", "plot_FIRMS_ts" ], "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 MODIS over Southern Africa." ] }, { "cell_type": "code", "id": "699458dc", "metadata": { "tags": [ "hide-input" ] }, "source": [ "plot_FIRMS_geo = fires.hvplot.points('longitude', 'latitude', groupby='index.day_of_year', c='frp', geo=True, alpha=0.2,\n", " tiles='ESRI', xlim=(llx, urx), ylim=(lly, ury), cmap='autumn', logz=True,\n", " dynamic=False)\n", "plot_FIRMS_geo" ], "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 MODIS region using `hvplot`." ] }, { "cell_type": "code", "id": "071602fc", "metadata": { "tags": [ "hide-input" ] }, "source": [ "modis_dataset = modis_scn.to_xarray_dataset()\n", "\n", "img = get_enhanced_image(modis_scn['true_color'].squeeze())\n", "img = img.data\n", "img = img.clip(0, 1)\n", "img = (img * 255).astype(np.uint8)\n", "\n", "modis_dataset['true_color'] = img\n", "\n", "grid = modis_scn.min_area().get_lonlats()\n", "lons, lats = grid[0][0], grid[1][:, 0]\n", "modis_dataset = modis_dataset.assign_coords(dict(x=lons, y=lats))\n", "modis_dataset = modis_dataset.rename(dict(x='longitude', y='latitude'))\n", "\n", "plot_MODIS_rgb = modis_dataset['true_color'].hvplot.rgb(x='longitude', y='latitude', bands='bands', rasterize=True)\n", "plot_MODIS_rgb" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "id": "13aa66a0", "metadata": {}, "source": [ "Get the fire pixels that are visible in the the MODIS scene by filtering the time, longitude & latitude. We can also compute which row and column the fire pixel is in for this projection." ] }, { "cell_type": "code", "id": "4ca5cbfe", "metadata": { "tags": [ "hide-input" ] }, "source": [ "grid = modis_scn.min_area().get_lonlats()\n", "\n", "# Mask any fire pixels not in this scene\n", "time_mask = np.logical_and(fires.timestamp >= scn.start_time,\n", " fires.timestamp <= scn.end_time)\n", "fire_points = fires.loc[time_mask]\n", "fire_points = fire_points.cx[grid[0].min():grid[0].max(), grid[1].min():grid[1].max()]\n", "\n", "# extract the x and y coordinates as flat arrays\n", "x = np.ravel(grid[0])\n", "y = np.ravel(grid[1])\n", "points = np.dstack([x, y]).squeeze()\n", "\n", "# Find locations of fire pixels within satellite image\n", "tree = cKDTree(points)\n", "distance, idx = tree.query(fire_points[['longitude', 'latitude']], k=1)\n", "index = np.unravel_index(idx, grid[0].shape)\n", "index = np.vstack(index).T\n", "index\n", "\n", "# fire_points[['y', 'x']] = index\n", "fire_points = fire_points[img.values[0, index[:, 0], index[:, 1]] != 0]\n", "fire_points" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "id": "4b7e704a", "metadata": {}, "source": [ "Now overlay the fire pixels ontop of the MODIS image, along with the FRP for each fire pixel. Try zooming in, you should be able to see clear smoke trails at the locations of some of the fires!" ] }, { "cell_type": "code", "id": "56e3e46d", "metadata": { "tags": [ "hide-input" ] }, "source": [ "rgb = modis_dataset['true_color'].hvplot.rgb(x='longitude', y='latitude', bands='bands', data_aspect=1, hover=False, title='True Colour', rasterize=True)\n", "points = fire_points.hvplot.points('longitude', 'latitude', c='frp', cmap='autumn', logz=True, alpha=0.4)\n", "plot_FIRMS_MODIS_rgb = rgb*points\n", "plot_FIRMS_MODIS_rgb" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "id": "c834983c", "metadata": {}, "source": [ "We can also overlay the fire pixels ontop of the MODIS 3.7 micron wavelength band. The 3.7 band is a thermal band. Fires will appear as very bright pixels in image. Try zooming in, you should be able to clearly see bright spots at the location of each fire pixel." ] }, { "cell_type": "code", "id": "3e4a9f6a", "metadata": { "tags": [ "hide-input" ] }, "source": [ "thermal = modis_dataset['20'].hvplot.image(x='longitude', y='latitude', cmap='viridis', data_aspect=1, hover=False, title='Band 20: 3.7 micron', rasterize=True)\n", "points = fire_points.hvplot.points('longitude', 'latitude', c='frp', cmap='autumn', logz=True, alpha=0.4)\n", "plot_FIRMS_MODIS_thermal = thermal*points\n", "plot_FIRMS_MODIS_thermal" ], "outputs": [], "execution_count": null }, { "cell_type": "markdown", "id": "8f0dc221", "metadata": {}, "source": [ "## Summary\n", "\n", "This notebook has demonstrated the use of:\n", "\n", "* The `satpy` package to easily load, plot and regrid satellite data from the MODIS sensor.\n", "* `hvplot` to interactively visualise wildfire pixels detected from the MODIS sensor.\n", "* `geopandas` to load, filter, and manipulate historical wildfire pixel data." ] }, { "cell_type": "markdown", "id": "85668bd9", "metadata": {}, "source": [ "## Citing this Notebook" ] }, { "cell_type": "markdown", "id": "828485b4", "metadata": {}, "source": "Please see [CITATION.cff](https://github.com/eds-book/a2875cdc-ba6a-49dc-aab7-cdf7c4fc0fa8/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": "7d0be357", "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/2).\n", "\n", "**Dataset**: MOD021KM & FIRMS (Collection 61) {cite:p}`firms`.\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": "94cec629", "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 }