{ "cells": [ { "cell_type": "markdown", "id": "fdde40cf-4514-4432-9a82-7c7c0e4585ba", "metadata": {}, "source": [ "### MTG Fire-Temperature RGB visualization using SATPy with compressed MTG data from EUMETCast and the FRP product " ] }, { "cell_type": "markdown", "id": "5338013e-5289-4072-8758-350f0374c5b2", "metadata": {}, "source": [ "Notebook prepared by Ben Maathuis, ITC-University of Twente, Enschede. The Netherlands" ] }, { "cell_type": "markdown", "id": "17cd056a-5bc5-4ec8-8442-e26b9b03cf89", "metadata": {}, "source": [ "+ For MTG Fire-Temperature information, see also: https://user.eumetsat.int/catalogue/EO:EUM:DAT:1046\n", "+ For additional information on the MTG Fire Radiative Power product see also: https://user.eumetsat.int/catalogue/EO:EUM:DAT:1156" ] }, { "cell_type": "markdown", "id": "57cfe1f3-b807-4b57-b8f8-18d5c6ff4d0e", "metadata": {}, "source": [ "Sample data is available at: https://filetransfer.itc.nl/pub/52n/ilwis_py/sample_data_V2/Fire_data.zip. Unzip the file: note the content, 40 compressed segment files of a MTG time step (202607291500) in NetCDF format and a Fire Radiative Power file of the correponding time step. Here it is assumed that the unzipped data is situated in a folder '/Fire_data' which should be situated within this notebook folder! It is furthermore assumed that you have locally installed ILWIS386 for data visualization when data is written to disk." ] }, { "cell_type": "code", "execution_count": null, "id": "87a2fecd-3808-4d2a-969d-db1a281fa708", "metadata": {}, "outputs": [], "source": [ "#uncomment the line below to install the HDF5plugin (if not already installed), note: the MTG data delivered through EUMETCast is compressed\n", "#and requires the fcidecomp - see in your python folder: \\Lib\\site-packages\\hdf5plugin\\plugins\\libh5fcidecomp.dll\n", "#!pip install hdf5plugin " ] }, { "cell_type": "code", "execution_count": null, "id": "86dea368-8f28-44a5-ac74-ccbf8d6901b1", "metadata": {}, "outputs": [], "source": [ "#set the environment so python knows the hdf5plugin location\n", "import os\n", "import site\n", "\n", "folder = [\n", " os.path.join(f, \"hdf5plugin\", \"plugins\")\n", " for f in site.getsitepackages()\n", " if os.path.exists(os.path.join(f, \"hdf5plugin\", \"plugins\"))\n", "]\n", "\n", "if folder:\n", " os.environ[\"HDF5_PLUGIN_PATH\"] = folder[0]\n", " print(os.environ[\"HDF5_PLUGIN_PATH\"])\n", "else:\n", " raise RuntimeError(\"hdf5plugin not found\")" ] }, { "cell_type": "code", "execution_count": null, "id": "71d0888b-d22b-4649-8ac8-b7776ed54ffa", "metadata": {}, "outputs": [], "source": [ "#import the required resources\n", "import ilwis\n", "import shutil\n", "import sys\n", "import glob\n", "from satpy.scene import Scene\n", "from satpy.composites.core import GenericCompositor\n", "from satpy import find_files_and_readers\n", "from satpy.enhancements.enhancer import get_enhanced_image\n", "from datetime import datetime\n", "import hdf5plugin\n", "import numpy as np\n", "import pandas as pd\n", "import matplotlib.pyplot as plt\n", "from matplotlib.colors import LogNorm\n", "import matplotlib.ticker as mticker\n", "import netCDF4 as nc\n", "import h5py\n", "import geopandas as gpd\n", "from shapely.geometry import Point\n", "import cartopy\n", "from cartopy import crs as ccrs, feature as cfeature\n", "from cartopy.io.shapereader import Reader\n", "\n", "import warnings\n", "warnings.filterwarnings('ignore')" ] }, { "cell_type": "code", "execution_count": null, "id": "8c5b86eb-52f8-4a6f-8aae-7d17a1640246", "metadata": {}, "outputs": [], "source": [ "#set ypour folders\n", "MTG_dir = os.getcwd()+'/Fire_data'\n", "\n", "print(\"current dir is: %s\" % (os.getcwd()))\n", "print(\"current input data directory is: \",MTG_dir) \n", "\n", "#local output folder\n", "dst_dir = os.getcwd()+'/Fire_result'\n", "\n", "print(\"current dir is: %s\" % (os.getcwd()))\n", "print(\"current working directory is:\",dst_dir) \n", "\n", "if os.path.isdir(dst_dir):\n", " print(\"Folder exists\")\n", "else:\n", " print(\"Folder doesn't exists\")\n", " os.mkdir(dst_dir)\n", "\n", "print(dst_dir)" ] }, { "cell_type": "code", "execution_count": null, "id": "e1c31d50-7300-4bd4-b159-1177c30d0b49", "metadata": {}, "outputs": [], "source": [ "#set the working directory for ILWISPy\n", "ilwis.setWorkingCatalog(dst_dir)\n", "print(dst_dir)" ] }, { "cell_type": "code", "execution_count": null, "id": "956aa085-38d2-4f24-94ee-5a6216b4c370", "metadata": {}, "outputs": [], "source": [ "#the timestamp of the sample data\n", "year = 2026\n", "month = 7\n", "day = 29\n", "hour = 15 #only the hour timestamp is required" ] }, { "cell_type": "code", "execution_count": null, "id": "c8a7519b-a04d-4593-a721-8decd5fb7aca", "metadata": {}, "outputs": [], "source": [ "date = datetime(year, month, day, hour)\n", "\n", "date_str = date.strftime(\"%Y%m%d%H\")\n", "print(date_str)" ] }, { "cell_type": "code", "execution_count": null, "id": "fe600d17-bff5-4c75-9972-81f0a2854125", "metadata": {}, "outputs": [], "source": [ "#copy the data from the source data folder to your destination folder\n", "for seg in range(0, 41): # 32-38 (segments covering the Mediterrenean region only)\n", " pattern = os.path.join(MTG_dir, f\"*_{date_str}*_{seg:04d}.nc\")\n", "\n", " for file in glob.glob(pattern):\n", " shutil.copy2(file, dst_dir)\n", " print(f\"Copied: {os.path.basename(file)}\")" ] }, { "cell_type": "code", "execution_count": null, "id": "fcb1cddd-07ac-413f-a939-9f90a3fe4da3", "metadata": {}, "outputs": [], "source": [ "#set the environment of satpy\n", "files = find_files_and_readers(base_dir=dst_dir, reader='fci_l1c_nc')" ] }, { "cell_type": "code", "execution_count": null, "id": "fbdc8694-2a74-4e72-a6b4-0fe1c0d030c7", "metadata": {}, "outputs": [], "source": [ "#add the MTG image in satpy memory\n", "scn_tc = Scene(filenames=files)" ] }, { "cell_type": "code", "execution_count": null, "id": "afbae158-0b4e-4814-9bce-7eb543edf1aa", "metadata": {}, "outputs": [], "source": [ "#print(scn.available_dataset_names()) #uncomment if you want to see the content of the imported image in satpy" ] }, { "cell_type": "markdown", "id": "cec286a9-76ac-4717-ae1b-afa7b2ad1ca8", "metadata": {}, "source": [ "### First create a natural color composite\n", "Get an idea of the location of clouds and land features, also to see if smoke plumes (resulting from fires) can be identified. Note the MTG NR data files for the Visible channels are having 11136 lines / columns, so be a bit patient while executing the code fields below!" ] }, { "cell_type": "code", "execution_count": null, "id": "a15819f8-261f-4b32-86dc-8ae37a724508", "metadata": {}, "outputs": [], "source": [ "# scn.load(['true_color'], upper_right_corner='NE')\n", "# #scn.show('natural_color') #uncomment if you want to see the full MTG (daytime) image\n", "# nc = scn['true_color'].data" ] }, { "cell_type": "code", "execution_count": null, "id": "dd7f038c-745f-4a44-b911-0a7cc2330e63", "metadata": {}, "outputs": [], "source": [ "#load data for selected color visualization\n", "image = 'true_color'\n", "scn_tc.load([image], upper_right_corner='NE')" ] }, { "cell_type": "code", "execution_count": null, "id": "92713fde-168a-4ba3-9f1b-2812a37618c5", "metadata": {}, "outputs": [], "source": [ "plt.figure(figsize=(8,8))\n", "img = get_enhanced_image(scn_tc['true_color'])\n", "# get DataArray out of `XRImage` object\n", "img_data = img.data\n", "#comment the line below if you don't want to see the full MTG (daytime) image\n", "img_data.plot.imshow(rgb='bands', vmin=0, vmax=1); #uncomment if you want to see the full MTG (daytime) image" ] }, { "cell_type": "code", "execution_count": null, "id": "b88a8694-d3cc-4c22-a437-004bd6882ac8", "metadata": {}, "outputs": [], "source": [ "#rotate the array and create a 1D data array\n", "data1_rot0 = np.rot90(img_data[0], k=4, axes=(1, 0))\n", "data1_rot1 = np.rot90(img_data[1], k=4, axes=(1, 0))\n", "data1_rot2 = np.rot90(img_data[2], k=4, axes=(1, 0))\n", "\n", "data1_rot = np.array([data1_rot2, data1_rot1, data1_rot0]).flatten()" ] }, { "cell_type": "code", "execution_count": null, "id": "0ddb9017-6b15-4132-9ce2-2c172d4a2b9c", "metadata": {}, "outputs": [], "source": [ "#create an empty raster\n", "grf = ilwis.GeoReference('code=georef:type=corners,csy=proj4:+proj=geos +h=35786400 +a=6378137 +rf=298.257223563,envelope=-5567999.9986 5567999.9986 5567999.9986 -5567999.9986,gridsize=11136 11136,cornerofcorners=yes')\n", "dfNum = ilwis.DataDefinition(ilwis.NumericDomain('code=value'), ilwis.NumericRange(-100.0, 500.0, 0.001))\n", "rcNew = ilwis.RasterCoverage()\n", "rcNew.setSize(ilwis.Size(11136,11136,3))\n", "rcNew.setGeoReference(grf)\n", "rcNew.setDataDef(dfNum)" ] }, { "cell_type": "code", "execution_count": null, "id": "b521c31c-491a-4a15-acf6-b90912500d71", "metadata": {}, "outputs": [], "source": [ "#add the data to the ilwis raster\n", "rcNew.array2raster(data1_rot)" ] }, { "cell_type": "code", "execution_count": null, "id": "e153e80e-a8c8-4431-a506-f5092a1a0626", "metadata": {}, "outputs": [], "source": [ "#remove if (satpy) nodata\n", "rcNew1 = ilwis.do('mapcalc', 'iff(@1>=0,@1,0)', rcNew) " ] }, { "cell_type": "code", "execution_count": null, "id": "ed01304c-ccba-4141-b959-f6762b53d6c6", "metadata": {}, "outputs": [], "source": [ "#store the results as an ilwis maplist - display the map using the ilwis386 desktop software\n", "rcNew1.store('finalTC.mpl')" ] }, { "cell_type": "code", "execution_count": null, "id": "cef8b729-a12b-4e39-a8cb-ce914c79d56a", "metadata": {}, "outputs": [], "source": [ "#define AoI - here Mediterranean region\n", "#note pixelsize Eastings total = 50 degree (-11 to 39) * 110 km per degree / No columns (9001) = roughly 600 meters/pixel\n", "grf_LL= ilwis.GeoReference('code=georef:type=corners, csy=epsg:4326, envelope=-11 55 39 34, gridsize=9001 3781, cornerofcorners=yes')" ] }, { "cell_type": "code", "execution_count": null, "id": "bf489ffb-38c7-40a5-afec-24d0575c39dd", "metadata": {}, "outputs": [], "source": [ "# Resample the map list and store the results\n", "rcNew_med = ilwis.do(\"resample\", rcNew1, grf_LL, \"nearestneighbour\")" ] }, { "cell_type": "code", "execution_count": null, "id": "d7a844e5-90dd-4d3a-b547-7c35c6a87d22", "metadata": {}, "outputs": [], "source": [ "#stretch the rasterbands contained in rcNew_med using a loop\n", "multiple_stretch = []\n", "multiple_bands = ilwis.do('selection',rcNew_med,\"rasterbands(0..2)\") \n", "ls = ilwis.do('linearstretch',multiple_bands, 1) #using an upper and lower data limit defined by the cumulative 1 and 99 % thresholds\n", "mb_stretch = ilwis.do('setvaluerange', ls, 0, 255, 1) #set ouput to byte range" ] }, { "cell_type": "code", "execution_count": null, "id": "bb7f1f76-2078-44c4-890e-922984d757c3", "metadata": {}, "outputs": [], "source": [ "#store the results as an ILWIS maplist - display the map using the ilwis386 desktop software\n", "mb_stretch.store('med_'+str(date_str)+'00.mpl')" ] }, { "cell_type": "code", "execution_count": null, "id": "2916b942-5cce-4062-9d92-467a61bc862b", "metadata": {}, "outputs": [], "source": [ "#load the individual spectral channels\n", "Blues = ilwis.do('selection',mb_stretch,\"rasterbands(0)\")\n", "Greens = ilwis.do('selection',mb_stretch,\"rasterbands(1)\")\n", "Reds = ilwis.do('selection',mb_stretch,\"rasterbands(2)\")" ] }, { "cell_type": "code", "execution_count": null, "id": "d5f1a839-8606-4662-b084-de2b66a165d7", "metadata": {}, "outputs": [], "source": [ "#transform the spectral channels from iwlsi format to a numpy array using the iterator\n", "Blues_2np = np.fromiter(iter(Blues), np.ubyte, Blues.size().linearSize()) \n", "Blues_2np = Blues_2np.reshape((Blues.size().ysize, Blues.size().xsize))\n", "\n", "Greens_2np = np.fromiter(iter(Greens), np.ubyte, Blues.size().linearSize()) \n", "Greens_2np = Greens_2np.reshape((Blues.size().ysize, Blues.size().xsize))\n", "\n", "Reds_2np = np.fromiter(iter(Reds), np.ubyte, Blues.size().linearSize()) \n", "Reds_2np = Reds_2np.reshape((Blues.size().ysize, Blues.size().xsize))" ] }, { "cell_type": "code", "execution_count": null, "id": "69f51657-e5c3-42f6-8065-5336b43b338e", "metadata": {}, "outputs": [], "source": [ "#create a numpy 3D data stack\n", "ncol = np.dstack((Reds_2np, Greens_2np, Blues_2np))" ] }, { "cell_type": "code", "execution_count": null, "id": "803c8803-e111-4265-8212-e3a8ad4442c4", "metadata": {}, "outputs": [], "source": [ "#plot the result using matplotlib / cartopy\n", "#note the features, like clouds and smoke plumes (Iberian Peninsula, Greece - Crete and along the south-west coast of Türkiye)\n", "img_extent = (-11, 39, 34, 55) \n", "fig = plt.figure(figsize=(12, 10))\n", "ax = plt.axes(projection=ccrs.PlateCarree())\n", "plt.title('Meteosat Third Generation composite of '+ (date_str)+'00')\n", "\n", "#add Natural earth shape files\n", "ax.add_feature(cfeature.BORDERS, linewidth=0.5, color='yellow', zorder=3)\n", "ax.add_feature(cfeature.COASTLINE, linewidth=0.5, color='blue', zorder=3)\n", "\n", "#data raster \n", "ax.imshow(ncol, origin='upper', extent=img_extent, transform=ccrs.PlateCarree())\n", "\n", "gl = ax.gridlines(draw_labels=True,color = 'grey', linestyle = '--', linewidth = 0.5)\n", "gl.top_labels = False\n", "gl.right_labels = False\n", "gl.ylocator = mticker.FixedLocator([35, 40, 45, 50])\n", "\n", "plt.show()\n", "#fig.savefig(dst_dir+'/mtg_tc'+(date_str)+'00.jpg')" ] }, { "cell_type": "markdown", "id": "46d341d0-dd57-4f17-815a-bb5eae8b593c", "metadata": {}, "source": [ "### Create the MTG Fire Temperature visualization" ] }, { "cell_type": "markdown", "id": "9819349d-ee7b-4707-8f63-9cfdc061a3c4", "metadata": {}, "source": [ "#### Process the ir_38 channel (Red)" ] }, { "cell_type": "code", "execution_count": null, "id": "64eef86f-1f4b-409d-b3aa-9c78db60b064", "metadata": {}, "outputs": [], "source": [ "scn_tc.load(['ir_38'])\n", "#uncomment the line below to show the spectral channel\n", "#scn_tc.show('ir_38')" ] }, { "cell_type": "code", "execution_count": null, "id": "943b6080-83cf-4fd6-8092-3678139e2702", "metadata": {}, "outputs": [], "source": [ "ir_38_values = scn_tc['ir_38'].values\n", "# uncomment the line below to retrieve the value for a specific location in the array \n", "#print(ir_38_values[3000,3000])" ] }, { "cell_type": "code", "execution_count": null, "id": "08eccd38-1a89-4f82-9cc3-d49a24e3c1e8", "metadata": {}, "outputs": [], "source": [ "#Rotate / flip the image\n", "ir_38_rot = np.flipud(ir_38_values)" ] }, { "cell_type": "code", "execution_count": null, "id": "d89f84ae-7459-449d-813b-2619aa6b155a", "metadata": {}, "outputs": [], "source": [ "#Create a 1D data list\n", "data38 = np.array([ir_38_rot]).flatten()\n", "data38.shape" ] }, { "cell_type": "code", "execution_count": null, "id": "ed56f699-6dba-4cab-b046-0e4474a0c3e4", "metadata": {}, "outputs": [], "source": [ "#create an empty raster - note TIR = 5568 lines and columns - 2 km SSP\n", "grf = ilwis.GeoReference('code=georef:type=corners,csy=proj4:+proj=geos +h=35786400 +a=6378137 +rf=298.257223563,envelope=-5567999.9942 5567999.9942 5567999.9942 -5567999.9942,gridsize=5568 5568,cornerofcorners=yes')\n", "dfNum = ilwis.DataDefinition(ilwis.NumericDomain('code=value'), ilwis.NumericRange(0.0, 500.0, 0.001))\n", "rcNew = ilwis.RasterCoverage()\n", "rcNew.setSize(ilwis.Size(5568,5568,1))\n", "rcNew.setGeoReference(grf)\n", "rcNew.setDataDef(dfNum)" ] }, { "cell_type": "code", "execution_count": null, "id": "2aecca49-dae7-4028-a41e-daf127f2f840", "metadata": {}, "outputs": [], "source": [ "rcNew.array2raster(data38)" ] }, { "cell_type": "code", "execution_count": null, "id": "c03179f0-e9a5-4b8a-a07e-1e1f8ac5b335", "metadata": {}, "outputs": [], "source": [ "#resmaple to AoI - here Mediterranean region\n", "rc38_res = ilwis.do('resample', rcNew, grf_LL, 'nearestneighbour')\n", "rc38_res = ilwis.do('mapcalc', 'iff(@1>=0,@1,0)', rc38_res) # remove (satpy) nodata\n", "print(rc38_res.size())" ] }, { "cell_type": "markdown", "id": "a876f44b-5d8e-49a2-8a16-4087afd4d0b3", "metadata": {}, "source": [ "#### Process the nir_22 channel (Green)" ] }, { "cell_type": "code", "execution_count": null, "id": "3bb3e118-58c5-4cbe-9022-d8e3d06c0eb2", "metadata": {}, "outputs": [], "source": [ "scn_tc.load([\"nir_22\"])\n", "# you can access the values of a dataset as a Numpy array with\n", "nir_22_values = scn_tc['nir_22'].values\n", "# uncomment to retrieve the value for a specific location in the array \n", "#print(nir_22_values[3000,3000]) #uncomment to retrieve the value for the specified location\n", "#uncomment to show the image\n", "#scn_tc.show('nir_22') " ] }, { "cell_type": "code", "execution_count": null, "id": "6a3a20ad-047b-45bc-a6e1-1d9b519cb943", "metadata": {}, "outputs": [], "source": [ "nir_22_rot = np.flipud(nir_22_values)" ] }, { "cell_type": "code", "execution_count": null, "id": "18c17cb7-577f-4007-a7cd-422ef4d6e025", "metadata": {}, "outputs": [], "source": [ "data22 = np.array([nir_22_rot]).flatten()" ] }, { "cell_type": "code", "execution_count": null, "id": "79ea52f2-c94b-49d2-8180-7aa86756d4bc", "metadata": {}, "outputs": [], "source": [ "#create an empty raster - note NIR = 11136 lines and columns - 1 km SSP\n", "grf = ilwis.GeoReference('code=georef:type=corners,csy=proj4:+proj=geos +h=35786400 +a=6378137 +rf=298.257223563,envelope=-5567999.9986 5567999.9986 5567999.9986 -5567999.9986,gridsize=11136 11136,cornerofcorners=yes')\n", "dfNum = ilwis.DataDefinition(ilwis.NumericDomain('code=value'), ilwis.NumericRange(-100.0, 500.0, 0.001))\n", "rcNew = ilwis.RasterCoverage()\n", "rcNew.setSize(ilwis.Size(11136,11136,1))\n", "rcNew.setGeoReference(grf)\n", "rcNew.setDataDef(dfNum)" ] }, { "cell_type": "code", "execution_count": null, "id": "92a9f686-5912-487b-b646-9dab1b78f6f9", "metadata": {}, "outputs": [], "source": [ "rcNew.array2raster(data22)" ] }, { "cell_type": "code", "execution_count": null, "id": "0076527b-7c45-4475-893b-82aa30003c5c", "metadata": {}, "outputs": [], "source": [ "rc22_res = ilwis.do('resample', rcNew, grf_LL, 'nearestneighbour')\n", "rc22_res = ilwis.do('mapcalc', 'iff(@1>=0,@1,0)', rc22_res) # remove (satpy) nodata\n", "print(rc22_res.size())" ] }, { "cell_type": "markdown", "id": "fcb2f10a-9fc1-4314-ac72-e0d97625a211", "metadata": {}, "source": [ "#### Process the nir_16 channel (Blue)" ] }, { "cell_type": "code", "execution_count": null, "id": "f6dce7c8-8bd3-4f39-842a-ab62fcddeb49", "metadata": {}, "outputs": [], "source": [ "scn_tc.load([\"nir_16\"])\n", "# you can access the values of a dataset as a Numpy array with\n", "nir_16_values = scn_tc['nir_16'].values\n", "# retrieve the value for a specific location in the array \n", "#uncomment to print the pixel value\n", "#(nir_16_values[3000,3000])\n", "#uncomment to shwo the image\n", "#scn_tc.show('nir_16')" ] }, { "cell_type": "code", "execution_count": null, "id": "cf0d6826-31c0-4e06-b43f-fdbf5bd5d399", "metadata": {}, "outputs": [], "source": [ "nir_16_rot = np.flipud(nir_16_values)" ] }, { "cell_type": "code", "execution_count": null, "id": "6dba5908-b0d4-4906-9448-574cc9f95d91", "metadata": {}, "outputs": [], "source": [ "data16 = np.array([nir_16_rot]).flatten()" ] }, { "cell_type": "code", "execution_count": null, "id": "a80fd5ac-d201-4e0b-b37b-c548cea01be2", "metadata": {}, "outputs": [], "source": [ "rcNew.array2raster(data16)" ] }, { "cell_type": "code", "execution_count": null, "id": "cf660139-1f98-4c71-86e1-b0f16f05ae11", "metadata": {}, "outputs": [], "source": [ "rc16_res = ilwis.do('resample', rcNew, grf_LL, 'nearestneighbour')\n", "rc16_res = ilwis.do('mapcalc', 'iff(@1>=0,@1,0)', rc16_res) # remove (satpy) nodata\n", "print(rc16_res.size())" ] }, { "cell_type": "markdown", "id": "6c181f5b-11ab-4964-bd48-82ac0ab51a33", "metadata": {}, "source": [ "### Image enhancement and creating a color composite\n", "\n", "Take the resampled images (ir_38, nri_22 and nir 16), conduct a linear constrast stretch, using a 1 % cut-off threshold and stretch the data to a byte range (0-255), subsequently transform the result into a dictionary containing the 3 numpy arrays, which are used later to display the Fire-Temperature images as an RGB using maplotlib" ] }, { "cell_type": "code", "execution_count": null, "id": "7b26c932-4494-4e45-88d2-1bcc35337727", "metadata": {}, "outputs": [], "source": [ "# Dictionary containing the input raster maps\n", "rasters = {\n", " \"rc38s_2np\": rc38_res,\n", " \"rc22s_2np\": rc22_res,\n", " \"rc16s_2np\": rc16_res,\n", "}\n", "\n", "# Output dictionary\n", "rasters_np = {}\n", "\n", "for name, raster in rasters.items():\n", "\n", " # Create histogram\n", " hist = raster.statistics(ilwis.PropertySets.pHISTOGRAM, 65535)\n", "\n", " # Calculate 1% stretch\n", " minPerc, maxPerc = hist.calcStretchRange(1)\n", "\n", " # Stretch and set value range\n", " stretched = ilwis.do('linearstretch', raster, minPerc, maxPerc)\n", " stretched = ilwis.do('setvaluerange', stretched, 0, 255, 1)\n", "\n", " # Convert to NumPy array\n", " arr = np.fromiter(iter(stretched), np.ubyte, stretched.size().linearSize())\n", " arr = arr.reshape((stretched.size().ysize, stretched.size().xsize))\n", "\n", " # Store result\n", " rasters_np[name] = arr" ] }, { "cell_type": "code", "execution_count": null, "id": "6395182a-16b5-4aab-8c00-3cbf65fc542c", "metadata": {}, "outputs": [], "source": [ "#create a 3D stack for visualization using matplotlib / cartopy\n", "rgb = np.dstack((\n", " rasters_np[\"rc38s_2np\"], # R\n", " rasters_np[\"rc22s_2np\"], # G\n", " rasters_np[\"rc16s_2np\"] # B\n", "))" ] }, { "cell_type": "code", "execution_count": null, "id": "527bdfd6-c158-4563-b1b4-6eb160879bb7", "metadata": {}, "outputs": [], "source": [ "img_extent = (-11, 39, 34, 55) \n", "fig = plt.figure(figsize=(12, 10))\n", "ax = plt.axes(projection=ccrs.PlateCarree())\n", "plt.title('Meteosat Third Generation Fire Temperature composite of '+ (date_str)+'00')\n", "\n", "#add Natural earth shape files\n", "ax.add_feature(cfeature.BORDERS, linewidth=0.5, color='yellow', zorder=3)\n", "ax.add_feature(cfeature.COASTLINE, linewidth=0.5, color='blue', zorder=3)\n", "\n", "#data raster \n", "ax.imshow(rgb, origin='upper', extent=img_extent, transform=ccrs.PlateCarree())\n", "\n", "gl = ax.gridlines(draw_labels=True,color = 'grey', linestyle = '--', linewidth = 0.5)\n", "gl.top_labels = False\n", "gl.right_labels = False\n", "gl.ylocator = mticker.FixedLocator([35, 40, 45, 50])\n", "\n", "plt.show()\n", "#uncomment line below to save as a JPG image\n", "#fig.savefig(dst_dir+'/mtg_FT'+(date_str)+'00.jpg')" ] }, { "cell_type": "code", "execution_count": null, "id": "85aac68e-1ff4-4ae3-b5e9-8707410a142f", "metadata": {}, "outputs": [], "source": [ "#create a 1D list\n", "R = rasters_np[\"rc38s_2np\"]\n", "G = rasters_np[\"rc22s_2np\"]\n", "B = rasters_np[\"rc16s_2np\"]\n", "\n", "data = np.array([B, G, R]).flatten()\n", "data.shape" ] }, { "cell_type": "code", "execution_count": null, "id": "9855c40b-3f14-424e-a029-94532358b9eb", "metadata": {}, "outputs": [], "source": [ "#create empty raster \n", "rcAll_ilw = ilwis.RasterCoverage()\n", "dfNumrc = ilwis.DataDefinition(ilwis.NumericDomain('code=value'), ilwis.NumericRange(0, 255, 0))\n", "rcAll_ilw.setDataDef(dfNumrc)\n", "rcAll_ilw.setSize(ilwis.Size(9001, 3781, 3))\n", "rcAll_ilw.setGeoReference(rc38_res.geoReference())" ] }, { "cell_type": "code", "execution_count": null, "id": "f358368b-e93a-4d2f-aa92-c15d828854cf", "metadata": {}, "outputs": [], "source": [ "#add the data to the raster\n", "rcAll_ilw.array2raster(data) \n", "print(rcAll_ilw.size())" ] }, { "cell_type": "code", "execution_count": null, "id": "74a812f6-13ef-407b-90a9-acb349d9d201", "metadata": {}, "outputs": [], "source": [ "#store the results as an ILWIS maplist - display the map using the ilwis386 desktop software\n", "rcAll_ilw.store('FT_comp_'+date_str+'.mpl')" ] }, { "cell_type": "markdown", "id": "864cca0c-7bdb-4c72-ba5e-cccaacf5d420", "metadata": {}, "source": [ "### Display the results for some selected areas" ] }, { "cell_type": "code", "execution_count": null, "id": "e2f6b587-6d91-4381-b75b-6ddb4398ddc6", "metadata": {}, "outputs": [], "source": [ "# Full image extent (same for both images)\n", "img_extent = (-11, 39, 34, 55)\n", "\n", "# Submap extent: lon_min, lon_max, lat_min, lat_max\n", "sub_extent = (-9.55, -4.10, 40.50, 43.79)\n", "\n", "fig, axes = plt.subplots(\n", " 1, 2,\n", " figsize=(16, 8),\n", " subplot_kw={\"projection\": ccrs.PlateCarree()}\n", ")\n", "\n", "titles = [\n", " f\"MTG Fire Temperature composite {date_str}00\",\n", " f\"MTG True Color composite {date_str}00\"\n", "]\n", "\n", "images = [rgb, ncol]\n", "\n", "for ax, img, title in zip(axes, images, titles):\n", "\n", " # Zoom to submap\n", " ax.set_extent(sub_extent, crs=ccrs.PlateCarree())\n", "\n", " # Add borders and coastlines\n", " ax.add_feature(cfeature.BORDERS, linewidth=0.5,\n", " color='yellow', zorder=3)\n", " ax.add_feature(cfeature.COASTLINE, linewidth=0.5,\n", " color='blue', zorder=3)\n", "\n", " # Plot image\n", " ax.imshow(\n", " img,\n", " origin='upper',\n", " extent=img_extent,\n", " transform=ccrs.PlateCarree()\n", " )\n", "\n", " # Gridlines\n", " gl = ax.gridlines(\n", " draw_labels=True,\n", " color='grey',\n", " linestyle='--',\n", " linewidth=0.5\n", " )\n", "\n", " gl.top_labels = False\n", " gl.right_labels = False\n", "\n", " gl.xlocator = mticker.FixedLocator([-8,-7, -6, - 5, -4])\n", " gl.ylocator = mticker.FixedLocator([41, 42, 43])\n", "\n", " ax.set_title(title)\n", "\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "id": "721ec7b4-df25-449e-9628-eb1fa50bb6d0", "metadata": {}, "outputs": [], "source": [ "# Full image extent (same for both images)\n", "img_extent = (-11, 39, 34, 55)\n", "\n", "# Submap extent: lon_min, lon_max, lat_min, lat_max\n", "sub_extent = (23.34, 30.12, 34.13, 40.32)\n", "\n", "fig, axes = plt.subplots(\n", " 1, 2,\n", " figsize=(16, 8),\n", " subplot_kw={\"projection\": ccrs.PlateCarree()}\n", ")\n", "\n", "titles = [\n", " f\"MTG Fire Temperature composite {date_str}00\",\n", " f\"MTG True Color composite {date_str}00\"\n", "]\n", "\n", "images = [rgb, ncol]\n", "\n", "for ax, img, title in zip(axes, images, titles):\n", "\n", " # Zoom to submap\n", " ax.set_extent(sub_extent, crs=ccrs.PlateCarree())\n", "\n", " # Add borders and coastlines\n", " ax.add_feature(cfeature.BORDERS, linewidth=0.5,\n", " color='yellow', zorder=3)\n", " ax.add_feature(cfeature.COASTLINE, linewidth=0.5,\n", " color='blue', zorder=3)\n", "\n", " # Plot image\n", " ax.imshow(\n", " img,\n", " origin='upper',\n", " extent=img_extent,\n", " transform=ccrs.PlateCarree()\n", " )\n", "\n", " # Gridlines\n", " gl = ax.gridlines(\n", " draw_labels=True,\n", " color='grey',\n", " linestyle='--',\n", " linewidth=0.5\n", " )\n", "\n", " gl.top_labels = False\n", " gl.right_labels = False\n", "\n", " gl.xlocator = mticker.FixedLocator([24, 26, 28, 30])\n", " gl.ylocator = mticker.FixedLocator([36, 38, 40])\n", "\n", " ax.set_title(title)\n", "\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "8d5e06c9-7427-478b-8ecb-32dc3b340c02", "metadata": {}, "source": [ "#### Process Fire Radiative Power\n", "\n", "See also: https://data.eumetsat.int/product/EO:EUM:DAT:1156\n", "\n", "file name convention = W_PT-LSASAF-LISBON,SATELLITE,LSA-509_MTG_MTFRPPIXEL_MTG-FD_C_LPMG_20260729150000.nc" ] }, { "cell_type": "code", "execution_count": null, "id": "b8f5bb53-1593-4c21-8031-fa9042f1d13b", "metadata": {}, "outputs": [], "source": [ "#check timestamp\n", "date_str" ] }, { "cell_type": "code", "execution_count": null, "id": "8b6c1ade-9807-4aa9-8e38-0c36520efca3", "metadata": {}, "outputs": [], "source": [ "for filename in glob.glob(os.path.join(MTG_dir+'/W_PT-LSASAF-LISBON,SATELLITE,LSA-509_MTG_MTFRPPIXEL_MTG-FD_C_LPMG_'+(date_str)+'0000.nc')):\n", " shutil.copy(filename, dst_dir)" ] }, { "cell_type": "code", "execution_count": null, "id": "071dc983-3b18-42e6-9fbd-e75cba038467", "metadata": {}, "outputs": [], "source": [ "fn = dst_dir+'/W_PT-LSASAF-LISBON,SATELLITE,LSA-509_MTG_MTFRPPIXEL_MTG-FD_C_LPMG_'+(date_str)+'0000.nc'\n", "ds = nc.Dataset(fn)\n", "\n", "print(ds.groups.keys())\n", "print(ds.groups[\"ListProduct\"].variables.keys())" ] }, { "cell_type": "code", "execution_count": null, "id": "097eeef8-3065-427f-bc23-bf1b23c873ec", "metadata": {}, "outputs": [], "source": [ "# Open the ListProduct group\n", "grp = ds.groups[\"ListProduct\"]\n", "\n", "# Create DataFrame\n", "df = pd.DataFrame({\n", " \"latitude\": grp.variables[\"LATITUDE\"][:],\n", " \"longitude\": grp.variables[\"LONGITUDE\"][:],\n", " \"fire_radiative_power\": grp.variables[\"FRP\"][:],\n", " \"fire_confidence\": grp.variables[\"FIRE_CONFIDENCE\"][:],\n", " \"frp_uncertainty\": grp.variables[\"FRP_UNCERTAINTY\"][:]\n", "})\n", "\n", "print(df.head())\n", "#print(df.info())" ] }, { "cell_type": "code", "execution_count": null, "id": "77af5a19-61dc-4013-9950-bf463e1b4239", "metadata": {}, "outputs": [], "source": [ "len(df)" ] }, { "cell_type": "code", "execution_count": null, "id": "6122ca93-4d05-4582-8d92-8e8d271ca8a9", "metadata": {}, "outputs": [], "source": [ "# Create point geometry from longitude and latitude\n", "geometry = gpd.points_from_xy(\n", " df[\"longitude\"],\n", " df[\"latitude\"]\n", ")\n", "\n", "# Create GeoDataFrame\n", "gdf = gpd.GeoDataFrame(\n", " df,\n", " geometry=geometry,\n", " crs=\"EPSG:4326\" # WGS84 latitude/longitude\n", ")\n", "\n", "# Save as ESRI Shapefile\n", "gdf.to_file(dst_dir+\"/MTG_FRP.shp\", driver=\"ESRI Shapefile\")\n", "\n", "print(f\"Created shapefile with {len(gdf)} fire detections.\")" ] }, { "cell_type": "code", "execution_count": null, "id": "fcc28b4d-979d-487b-97bc-38acd72ea33b", "metadata": {}, "outputs": [], "source": [ "# Filter to area of interest\n", "df_aoi = df[\n", " (df[\"latitude\"] >= 34.0) &\n", " (df[\"latitude\"] <= 55.0) &\n", " (df[\"longitude\"] >= -11.0) &\n", " (df[\"longitude\"] <= 39.0)\n", "].copy()\n", "\n", "print(f\"Number of fire pixels in AOI: {len(df_aoi)}\")\n", "print(df_aoi.head())" ] }, { "cell_type": "code", "execution_count": null, "id": "71170a64-409a-4b75-92c4-6e2f668a32ff", "metadata": {}, "outputs": [], "source": [ "# Create point geometry from longitude and latitude\n", "geometry = gpd.points_from_xy(\n", " df_aoi[\"longitude\"],\n", " df_aoi[\"latitude\"]\n", ")\n", "\n", "# Create GeoDataFrame\n", "gdf = gpd.GeoDataFrame(\n", " df_aoi,\n", " geometry=geometry,\n", " crs=\"EPSG:4326\" # WGS84 latitude/longitude\n", ")\n", "\n", "# Save as ESRI Shapefile\n", "gdf.to_file(dst_dir+\"/MTG_FRP_AOI.shp\", driver=\"ESRI Shapefile\")\n", "\n", "print(f\"Created shapefile with {len(gdf)} fire detections.\")" ] }, { "cell_type": "markdown", "id": "75bb4aee-623f-4ad8-877f-9a7d2f9d1907", "metadata": {}, "source": [ "#### Plot final results" ] }, { "cell_type": "code", "execution_count": null, "id": "e259d6e7-7981-4404-80d1-fde74572a251", "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(\n", " figsize=(12, 10),\n", " subplot_kw={\"projection\": ccrs.PlateCarree()}\n", ")\n", "\n", "extent = (-11, 39, 34, 55)\n", "\n", "# RGB image\n", "ax.imshow(\n", " rgb,\n", " extent=extent,\n", " origin=\"upper\",\n", " transform=ccrs.PlateCarree()\n", ")\n", "\n", "# Fire detections\n", "norm = LogNorm(\n", " vmin=max(1, df_aoi[\"fire_radiative_power\"].min()), # avoid zero for LogNorm\n", " vmax=df_aoi[\"fire_radiative_power\"].max()\n", ")\n", "\n", "sc = ax.scatter(\n", " df_aoi[\"longitude\"],\n", " df_aoi[\"latitude\"],\n", " c=df_aoi[\"fire_radiative_power\"],\n", " cmap=\"hot\",\n", " norm=norm,\n", " s=40,\n", " edgecolors=\"black\",\n", " linewidth=0.3,\n", " transform=ccrs.PlateCarree(),\n", " zorder=10\n", ")\n", "\n", "# Coastlines and borders\n", "ax.add_feature(cfeature.COASTLINE, linewidth=0.5)\n", "ax.add_feature(cfeature.BORDERS, linewidth=0.5)\n", "\n", "# Colorbar\n", "cbar = plt.colorbar(sc, ax=ax, shrink=0.3, extend=\"max\")\n", "cbar.set_label(\"Fire Radiative Power (MW)\")\n", "\n", "# Map extent\n", "ax.set_extent(extent, crs=ccrs.PlateCarree())\n", "\n", "# Gridlines\n", "gl = ax.gridlines(\n", " draw_labels=True,\n", " color=\"grey\",\n", " linestyle=\"--\",\n", " linewidth=0.5\n", ")\n", "\n", "gl.top_labels = False\n", "gl.right_labels = False\n", "\n", "gl.xlocator = mticker.FixedLocator([-10, -5, 0, 5, 10, 15, 20, 25, 30, 35])\n", "gl.ylocator = mticker.FixedLocator([35, 40, 45, 50, 55])\n", "\n", "ax.set_title(\"MTG Fire Temperature RGB with Fire Radiative Power Detections\")\n", "\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "67041743-c158-4f33-a6f4-c5a61f0ba308", "metadata": {}, "source": [ "### Plot some selected areas for better visualization" ] }, { "cell_type": "code", "execution_count": null, "id": "bd6d72fe-36be-4b0f-9da0-7e3f9cf9c690", "metadata": {}, "outputs": [], "source": [ "# ------------------------------------------------------------------\n", "# Image extent (full resampled image)\n", "# ------------------------------------------------------------------\n", "img_extent = (-11, 39, 34, 55)\n", "\n", "# ------------------------------------------------------------------\n", "# Create figure with two subplots\n", "# ------------------------------------------------------------------\n", "fig, axes = plt.subplots(\n", " 1, 2,\n", " figsize=(16, 8),\n", " subplot_kw={\"projection\": ccrs.PlateCarree()}\n", ")\n", "\n", "# Overall title\n", "fig.suptitle(\n", " f\"NW Iberian Peninsula MTG Fire Temperature RGB with Fire Radiative Power Detections\\nMTG FCI - {date_str}00 UTC\",\n", " fontsize=18,\n", " fontweight=\"bold\"\n", ")\n", "\n", "# Plot settings\n", "titles = [\n", " \"Fire Temperature RGB\",\n", " \"Natural Color RGB\"\n", "]\n", "\n", "images = [\n", " rgb,\n", " ncol\n", "]\n", "\n", "# Keep only positive FRP values for logarithmic colouring\n", "df_plot = df_aoi[df_aoi[\"fire_radiative_power\"] > 0]\n", "\n", "# ------------------------------------------------------------------\n", "# Plot both maps\n", "# ------------------------------------------------------------------\n", "for ax, image, title in zip(axes, images, titles):\n", "\n", " # Zoom to Iberian Peninsula\n", " ax.set_extent(\n", " [-9.00, -5.20, 40.50, 43.79],\n", " crs=ccrs.PlateCarree()\n", " )\n", " \n", " # Background RGB image\n", " ax.imshow(\n", " image,\n", " origin=\"upper\",\n", " extent=img_extent,\n", " transform=ccrs.PlateCarree()\n", " )\n", "\n", " # Coastlines and borders\n", " ax.add_feature(\n", " cfeature.COASTLINE,\n", " linewidth=0.5,\n", " color=\"blue\",\n", " zorder=3\n", " )\n", "\n", " ax.add_feature(\n", " cfeature.BORDERS,\n", " linewidth=0.5,\n", " color=\"yellow\",\n", " zorder=3\n", " )\n", "\n", " # Fire detections\n", " sc = ax.scatter(\n", " df_plot[\"longitude\"],\n", " df_plot[\"latitude\"],\n", " c=df_plot[\"fire_radiative_power\"],\n", " cmap=\"hot\",\n", " norm=LogNorm(\n", " vmin=df_plot[\"fire_radiative_power\"].min(),\n", " vmax=df_plot[\"fire_radiative_power\"].max()\n", " ),\n", " s=35,\n", " edgecolors=\"black\",\n", " linewidth=0.3,\n", " transform=ccrs.PlateCarree(),\n", " zorder=10\n", " )\n", "\n", " # Gridlines\n", " gl = ax.gridlines(\n", " draw_labels=True,\n", " color=\"grey\",\n", " linestyle=\"--\",\n", " linewidth=0.5\n", " )\n", "\n", " gl.top_labels = False\n", " gl.right_labels = False\n", "\n", " gl.xlocator = mticker.FixedLocator(\n", " [-9, -8, -7, -6, -5, -4, -2, 0]\n", " )\n", "\n", " gl.ylocator = mticker.FixedLocator(\n", " [37, 38, 39, 40, 41, 42, 43]\n", " )\n", "\n", " ax.set_title(title, fontsize=13)\n", "\n", "# -----------------------------------------------------------\n", "# Create a shared logarithmic colorbar\n", "# -----------------------------------------------------------\n", "\n", "# Set logarithmic scaling limits\n", "vmin = 1 # MW (must be > 0)\n", "vmax = df_plot[\"fire_radiative_power\"].max()\n", "\n", "# Update the scatter plots to use the same normalization\n", "norm = LogNorm(vmin=vmin, vmax=vmax)\n", "\n", "# In the plotting loop, use:\n", "sc = ax.scatter(\n", " df_plot[\"longitude\"],\n", " df_plot[\"latitude\"],\n", " c=df_plot[\"fire_radiative_power\"],\n", " cmap=\"hot\",\n", " norm=norm,\n", " s=35,\n", " edgecolors=\"black\",\n", " linewidth=0.3,\n", " transform=ccrs.PlateCarree(),\n", " zorder=10\n", ")\n", "\n", "# -----------------------------------------------------------\n", "# Create colorbar axis on the right\n", "# -----------------------------------------------------------\n", "cbar = fig.colorbar(\n", " sc,\n", " ax=axes,\n", " location=\"right\",\n", " shrink=0.55,\n", " pad=0.03,\n", " extend = 'max'\n", ")\n", "\n", "# Get current position\n", "pos = cbar.ax.get_position()\n", "\n", "# Move it upwards\n", "cbar.ax.set_position([\n", " pos.x0 + 0.14, # left\n", " pos.y0 + 0.04, # move up\n", " pos.width, # width\n", " pos.height # height\n", "])\n", "cbar.set_label(\"Fire Radiative Power (MW)\", fontsize=11)\n", "\n", "# Nice logarithmic tick locations\n", "cbar.set_ticks([1, 3, 10, 30, 100, 300, 1000])\n", "\n", "fig.subplots_adjust(\n", " top=0.96, # controls top margin\n", " wspace=0.10 # space between the two maps\n", ")\n", "\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "id": "d93faec6-a474-4dc6-a61a-588beb9dc84a", "metadata": {}, "outputs": [], "source": [ "# ------------------------------------------------------------------\n", "# Image extent (full resampled image)\n", "# ------------------------------------------------------------------\n", "img_extent = (-11, 39, 34, 55)\n", "\n", "# ------------------------------------------------------------------\n", "# Create figure with two subplots\n", "# ------------------------------------------------------------------\n", "fig, axes = plt.subplots(\n", " 1, 2,\n", " figsize=(16, 8),\n", " subplot_kw={\"projection\": ccrs.PlateCarree()}\n", ")\n", "\n", "# Overall title\n", "fig.suptitle(\n", " f\"Eastern Mediterranean MTG Fire Temperature RGB with Fire Radiative Power Detections\\nMTG FCI - {date_str}00 UTC\",\n", " fontsize=18,\n", " fontweight=\"bold\"\n", ")\n", "\n", "# Plot settings\n", "titles = [\n", " \"Fire Temperature RGB\",\n", " \"Natural Color RGB\"\n", "]\n", "\n", "images = [\n", " rgb,\n", " ncol\n", "]\n", "\n", "# Keep only positive FRP values for logarithmic colouring\n", "df_plot = df_aoi[df_aoi[\"fire_radiative_power\"] > 0]\n", "\n", "# ------------------------------------------------------------------\n", "# Plot both maps\n", "# ------------------------------------------------------------------\n", "for ax, image, title in zip(axes, images, titles):\n", "\n", " # Zoom to Eastern Mediterranean\n", " ax.set_extent(\n", " [23.34, 30.12, 34.13, 40.32],\n", " crs=ccrs.PlateCarree()\n", " )\n", " \n", " # Background RGB image\n", " ax.imshow(\n", " image,\n", " origin=\"upper\",\n", " extent=img_extent,\n", " transform=ccrs.PlateCarree()\n", " )\n", "\n", " # Coastlines and borders\n", " ax.add_feature(\n", " cfeature.COASTLINE,\n", " linewidth=0.5,\n", " color=\"blue\",\n", " zorder=3\n", " )\n", "\n", " ax.add_feature(\n", " cfeature.BORDERS,\n", " linewidth=0.5,\n", " color=\"yellow\",\n", " zorder=3\n", " )\n", "\n", " # Fire detections\n", " sc = ax.scatter(\n", " df_plot[\"longitude\"],\n", " df_plot[\"latitude\"],\n", " c=df_plot[\"fire_radiative_power\"],\n", " cmap=\"hot\",\n", " norm=LogNorm(\n", " vmin=df_plot[\"fire_radiative_power\"].min(),\n", " vmax=df_plot[\"fire_radiative_power\"].max()\n", " ),\n", " s=35,\n", " edgecolors=\"black\",\n", " linewidth=0.3,\n", " transform=ccrs.PlateCarree(),\n", " zorder=10\n", " )\n", "\n", " # Gridlines\n", " gl = ax.gridlines(\n", " draw_labels=True,\n", " color=\"grey\",\n", " linestyle=\"--\",\n", " linewidth=0.5\n", " )\n", "\n", " gl.top_labels = False\n", " gl.right_labels = False\n", "\n", " gl.xlocator = mticker.FixedLocator(\n", " [24, 25, 26, 27, 28, 29, 30]\n", " )\n", "\n", " gl.ylocator = mticker.FixedLocator(\n", " [35, 36, 37, 38, 39, 40]\n", " )\n", "\n", " ax.set_title(title, fontsize=13)\n", "\n", "# -----------------------------------------------------------\n", "# Create a shared logarithmic colorbar\n", "# -----------------------------------------------------------\n", "\n", "# Set logarithmic scaling limits\n", "vmin = 1 # MW (must be > 0)\n", "vmax = df_plot[\"fire_radiative_power\"].max()\n", "\n", "# Update the scatter plots to use the same normalization\n", "norm = LogNorm(vmin=vmin, vmax=vmax)\n", "\n", "# In the plotting loop, use:\n", "sc = ax.scatter(\n", " df_plot[\"longitude\"],\n", " df_plot[\"latitude\"],\n", " c=df_plot[\"fire_radiative_power\"],\n", " cmap=\"hot\",\n", " norm=norm,\n", " s=35,\n", " edgecolors=\"black\",\n", " linewidth=0.3,\n", " transform=ccrs.PlateCarree(),\n", " zorder=10\n", ")\n", "\n", "# -----------------------------------------------------------\n", "# Create colorbar axis on the right\n", "# -----------------------------------------------------------\n", "cbar = fig.colorbar(\n", " sc,\n", " ax=axes,\n", " location=\"right\",\n", " shrink=0.55,\n", " pad=0.03,\n", " extend = 'max'\n", ")\n", "\n", "# Get current position\n", "pos = cbar.ax.get_position()\n", "\n", "# Move it upwards\n", "cbar.ax.set_position([\n", " pos.x0 + 0.14, # left\n", " pos.y0 + 0.04, # move up\n", " pos.width, # width\n", " pos.height # height\n", "])\n", "cbar.set_label(\"Fire Radiative Power (MW)\", fontsize=11)\n", "\n", "# Nice logarithmic tick locations\n", "cbar.set_ticks([1, 3, 10, 30, 100, 300, 1000])\n", "\n", "\n", "fig.subplots_adjust(\n", " top=0.96, # controls top margin\n", " wspace=0.10 # space between the two maps\n", ")\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "37c3b254-0c40-43c8-95b7-2fe025c5ae8f", "metadata": {}, "source": [ "#### Remove the EUMETCast MTG segments and RFP from the destination folder\n", "\n", "New instance of kernel is started as files are still occupied by a previous instance of python" ] }, { "cell_type": "code", "execution_count": null, "id": "6b9982ca-b53a-4402-a748-39e1d09557e2", "metadata": {}, "outputs": [], "source": [ "import IPython\n", "IPython.Application.instance().kernel.do_shutdown(True) #automatically restarts kernel" ] }, { "cell_type": "code", "execution_count": null, "id": "6640fad2-918f-4ad7-b876-65520d82edad", "metadata": {}, "outputs": [], "source": [ "import sys, os \n", "import glob" ] }, { "cell_type": "code", "execution_count": null, "id": "7c83991f-dd1b-4dec-be10-855d4141b6d1", "metadata": {}, "outputs": [], "source": [ "#local output folder\n", "dst_dir = os.getcwd()+'/Fire_result'" ] }, { "cell_type": "code", "execution_count": null, "id": "83ad95cb-ff15-4da2-94af-8a34b9842a6c", "metadata": {}, "outputs": [], "source": [ "list = glob.glob(dst_dir+'/W_XX-EUMETSAT-Darmstadt,IMG+SAT,MTI1+FCI-1C-RRAD-FDHSI-FD--CHK-BODY--DIS-NC4E_C_EUMT_*.nc')\n", "list = [elem[len(dst_dir)+1:] for elem in list]\n", "for elem in list:\n", " os.remove(dst_dir+'/'+elem)" ] }, { "cell_type": "code", "execution_count": null, "id": "63ddafdf-6bc6-48ab-b43c-491dafb8cdfb", "metadata": {}, "outputs": [], "source": [ "list = glob.glob(dst_dir+'/W_PT-LSASAF-LISBON,SATELLITE,LSA-509_MTG_MTFRPPIXEL_MTG-FD_C_LPMG_*0000.nc')\n", "list = [elem[len(dst_dir)+1:] for elem in list]\n", "for elem in list:\n", " os.remove(dst_dir+'/'+elem)" ] }, { "cell_type": "code", "execution_count": null, "id": "b2d03d51-2d09-4b47-a5db-a85be5fc2ecb", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "id": "5fe0e73f-9b35-4e88-b7a8-7b73c0b527ab", "metadata": {}, "outputs": [], "source": [] } ], "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.13.3" } }, "nbformat": 4, "nbformat_minor": 5 }