MTG Fire-Temperature RGB visualization using SATPy with compressed MTG data from EUMETCast and the FRP product¶

Notebook prepared by Ben Maathuis, ITC-University of Twente, Enschede. The Netherlands

  • For MTG Fire-Temperature information, see also: https://user.eumetsat.int/catalogue/EO:EUM:DAT:1046
  • For additional information on the MTG Fire Radiative Power product see also: https://user.eumetsat.int/catalogue/EO:EUM:DAT:1156

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.

In [1]:
#uncomment the line below to install the HDF5plugin (if not already installed), note: the MTG data delivered through EUMETCast is compressed
#and requires the fcidecomp - see in your python folder: \Lib\site-packages\hdf5plugin\plugins\libh5fcidecomp.dll
#!pip install hdf5plugin
In [2]:
#set the environment so python knows the hdf5plugin location
import os
import site

folder = [
    os.path.join(f, "hdf5plugin", "plugins")
    for f in site.getsitepackages()
    if os.path.exists(os.path.join(f, "hdf5plugin", "plugins"))
]

if folder:
    os.environ["HDF5_PLUGIN_PATH"] = folder[0]
    print(os.environ["HDF5_PLUGIN_PATH"])
else:
    raise RuntimeError("hdf5plugin not found")

C:\python313\Lib\site-packages\hdf5plugin\plugins
In [3]:
#import the required resources
import ilwis
import shutil
import sys
import glob
from satpy.scene import Scene
from satpy.composites.core import GenericCompositor
from satpy import find_files_and_readers
from satpy.enhancements.enhancer import get_enhanced_image
from datetime import datetime
import hdf5plugin
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
import matplotlib.ticker as mticker
import netCDF4 as nc
import h5py
import geopandas as gpd
from shapely.geometry import Point
import cartopy
from cartopy import crs as ccrs, feature as cfeature
from cartopy.io.shapereader import Reader

import warnings
warnings.filterwarnings('ignore')
In [4]:
#set ypour folders
MTG_dir = os.getcwd()+'/Fire_data'

print("current dir is: %s" % (os.getcwd()))
print("current input data directory is: ",MTG_dir) 

#local output folder
dst_dir = os.getcwd()+'/Fire_result'

print("current dir is: %s" % (os.getcwd()))
print("current working directory is:",dst_dir) 

if os.path.isdir(dst_dir):
    print("Folder exists")
else:
    print("Folder doesn't exists")
    os.mkdir(dst_dir)

print(dst_dir)

current dir is: d:\jupyter\notebook_scripts\Special\NBI\notebook_final\gnc_notebooks
current input data directory is:  d:\jupyter\notebook_scripts\Special\NBI\notebook_final\gnc_notebooks/Fire_data
current dir is: d:\jupyter\notebook_scripts\Special\NBI\notebook_final\gnc_notebooks
current working directory is: d:\jupyter\notebook_scripts\Special\NBI\notebook_final\gnc_notebooks/Fire_result
Folder exists
d:\jupyter\notebook_scripts\Special\NBI\notebook_final\gnc_notebooks/Fire_result
In [5]:
#set the working directory for ILWISPy
ilwis.setWorkingCatalog(dst_dir)
print(dst_dir)

d:\jupyter\notebook_scripts\Special\NBI\notebook_final\gnc_notebooks/Fire_result
In [6]:
#the timestamp of the sample data
year = 2026
month = 7
day = 29
hour = 15 #only the hour timestamp is required
In [7]:
date = datetime(year, month, day, hour)

date_str = date.strftime("%Y%m%d%H")
print(date_str)

2026072915
In [8]:
#copy the data from the source data folder to your destination folder
for seg in range(0, 41):   # 32-38 (segments covering the Mediterrenean region only)
    pattern = os.path.join(MTG_dir, f"*_{date_str}*_{seg:04d}.nc")

    for file in glob.glob(pattern):
        shutil.copy2(file, dst_dir)
        print(f"Copied: {os.path.basename(file)}")

Copied: W_XX-EUMETSAT-Darmstadt,IMG+SAT,MTI1+FCI-1C-RRAD-FDHSI-FD--CHK-BODY--DIS-NC4E_C_EUMT_20260729150431_IDPFI_OPE_20260729150007_20260729150024_N_JLS_O_0091_0001.nc
Copied: W_XX-EUMETSAT-Darmstadt,IMG+SAT,MTI1+FCI-1C-RRAD-FDHSI-FD--CHK-BODY--DIS-NC4E_C_EUMT_20260729150438_IDPFI_OPE_20260729150007_20260729150031_N_JLS_O_0091_0002.nc
Copied: W_XX-EUMETSAT-Darmstadt,IMG+SAT,MTI1+FCI-1C-RRAD-FDHSI-FD--CHK-BODY--DIS-NC4E_C_EUMT_20260729150453_IDPFI_OPE_20260729150018_20260729150045_N_JLS_O_0091_0003.nc
Copied: W_XX-EUMETSAT-Darmstadt,IMG+SAT,MTI1+FCI-1C-RRAD-FDHSI-FD--CHK-BODY--DIS-NC4E_C_EUMT_20260729150503_IDPFI_OPE_20260729150024_20260729150100_N_JLS_O_0091_0004.nc
Copied: W_XX-EUMETSAT-Darmstadt,IMG+SAT,MTI1+FCI-1C-RRAD-FDHSI-FD--CHK-BODY--DIS-NC4E_C_EUMT_20260729150517_IDPFI_OPE_20260729150031_20260729150115_N_JLS_O_0091_0005.nc
Copied: W_XX-EUMETSAT-Darmstadt,IMG+SAT,MTI1+FCI-1C-RRAD-FDHSI-FD--CHK-BODY--DIS-NC4E_C_EUMT_20260729150520_IDPFI_OPE_20260729150045_20260729150124_N_JLS_O_0091_0006.nc
Copied: W_XX-EUMETSAT-Darmstadt,IMG+SAT,MTI1+FCI-1C-RRAD-FDHSI-FD--CHK-BODY--DIS-NC4E_C_EUMT_20260729150529_IDPFI_OPE_20260729150100_20260729150141_N_JLS_O_0091_0007.nc
Copied: W_XX-EUMETSAT-Darmstadt,IMG+SAT,MTI1+FCI-1C-RRAD-FDHSI-FD--CHK-BODY--DIS-NC4E_C_EUMT_20260729150536_IDPFI_OPE_20260729150107_20260729150149_N_JLS_O_0091_0008.nc
Copied: W_XX-EUMETSAT-Darmstadt,IMG+SAT,MTI1+FCI-1C-RRAD-FDHSI-FD--CHK-BODY--DIS-NC4E_C_EUMT_20260729150544_IDPFI_OPE_20260729150124_20260729150207_N_JLS_O_0091_0009.nc
Copied: W_XX-EUMETSAT-Darmstadt,IMG+SAT,MTI1+FCI-1C-RRAD-FDHSI-FD--CHK-BODY--DIS-NC4E_C_EUMT_20260729150557_IDPFI_OPE_20260729150141_20260729150225_N_JLS_O_0091_0010.nc
Copied: W_XX-EUMETSAT-Darmstadt,IMG+SAT,MTI1+FCI-1C-RRAD-FDHSI-FD--CHK-BODY--DIS-NC4E_C_EUMT_20260729150602_IDPFI_OPE_20260729150158_20260729150235_N_JLS_O_0091_0011.nc
Copied: W_XX-EUMETSAT-Darmstadt,IMG+SAT,MTI1+FCI-1C-RRAD-FDHSI-FD--CHK-BODY--DIS-NC4E_C_EUMT_20260729150612_IDPFI_OPE_20260729150207_20260729150254_N_JLS_O_0091_0012.nc
Copied: W_XX-EUMETSAT-Darmstadt,IMG+SAT,MTI1+FCI-1C-RRAD-FDHSI-FD--CHK-BODY--DIS-NC4E_C_EUMT_20260729150615_IDPFI_OPE_20260729150225_20260729150303_N_JLS_O_0091_0013.nc
Copied: W_XX-EUMETSAT-Darmstadt,IMG+SAT,MTI1+FCI-1C-RRAD-FDHSI-FD--CHK-BODY--DIS-NC4E_C_EUMT_20260729150623_IDPFI_OPE_20260729150244_20260729150322_N_JLS_O_0091_0014.nc
Copied: W_XX-EUMETSAT-Darmstadt,IMG+SAT,MTI1+FCI-1C-RRAD-FDHSI-FD--CHK-BODY--DIS-NC4E_C_EUMT_20260729150627_IDPFI_OPE_20260729150303_20260729150332_N_JLS_O_0091_0015.nc
Copied: W_XX-EUMETSAT-Darmstadt,IMG+SAT,MTI1+FCI-1C-RRAD-FDHSI-FD--CHK-BODY--DIS-NC4E_C_EUMT_20260729150635_IDPFI_OPE_20260729150322_20260729150352_N_JLS_O_0091_0016.nc
Copied: W_XX-EUMETSAT-Darmstadt,IMG+SAT,MTI1+FCI-1C-RRAD-FDHSI-FD--CHK-BODY--DIS-NC4E_C_EUMT_20260729150643_IDPFI_OPE_20260729150332_20260729150411_N_JLS_O_0091_0017.nc
Copied: W_XX-EUMETSAT-Darmstadt,IMG+SAT,MTI1+FCI-1C-RRAD-FDHSI-FD--CHK-BODY--DIS-NC4E_C_EUMT_20260729150648_IDPFI_OPE_20260729150352_20260729150421_N_JLS_O_0091_0018.nc
Copied: W_XX-EUMETSAT-Darmstadt,IMG+SAT,MTI1+FCI-1C-RRAD-FDHSI-FD--CHK-BODY--DIS-NC4E_C_EUMT_20260729150659_IDPFI_OPE_20260729150411_20260729150441_N_JLS_O_0091_0019.nc
Copied: W_XX-EUMETSAT-Darmstadt,IMG+SAT,MTI1+FCI-1C-RRAD-FDHSI-FD--CHK-BODY--DIS-NC4E_C_EUMT_20260729150708_IDPFI_OPE_20260729150421_20260729150501_N_JLS_O_0091_0020.nc
Copied: W_XX-EUMETSAT-Darmstadt,IMG+SAT,MTI1+FCI-1C-RRAD-FDHSI-FD--CHK-BODY--DIS-NC4E_C_EUMT_20260729150729_IDPFI_OPE_20260729150441_20260729150521_N_JLS_O_0091_0021.nc
Copied: W_XX-EUMETSAT-Darmstadt,IMG+SAT,MTI1+FCI-1C-RRAD-FDHSI-FD--CHK-BODY--DIS-NC4E_C_EUMT_20260729150747_IDPFI_OPE_20260729150501_20260729150540_N_JLS_O_0091_0022.nc
Copied: W_XX-EUMETSAT-Darmstadt,IMG+SAT,MTI1+FCI-1C-RRAD-FDHSI-FD--CHK-BODY--DIS-NC4E_C_EUMT_20260729150805_IDPFI_OPE_20260729150511_20260729150600_N_JLS_O_0091_0023.nc
Copied: W_XX-EUMETSAT-Darmstadt,IMG+SAT,MTI1+FCI-1C-RRAD-FDHSI-FD--CHK-BODY--DIS-NC4E_C_EUMT_20260729150823_IDPFI_OPE_20260729150530_20260729150619_N_JLS_O_0091_0024.nc
Copied: W_XX-EUMETSAT-Darmstadt,IMG+SAT,MTI1+FCI-1C-RRAD-FDHSI-FD--CHK-BODY--DIS-NC4E_C_EUMT_20260729150834_IDPFI_OPE_20260729150540_20260729150629_N_JLS_O_0091_0025.nc
Copied: W_XX-EUMETSAT-Darmstadt,IMG+SAT,MTI1+FCI-1C-RRAD-FDHSI-FD--CHK-BODY--DIS-NC4E_C_EUMT_20260729150853_IDPFI_OPE_20260729150600_20260729150648_N_JLS_O_0091_0026.nc
Copied: W_XX-EUMETSAT-Darmstadt,IMG+SAT,MTI1+FCI-1C-RRAD-FDHSI-FD--CHK-BODY--DIS-NC4E_C_EUMT_20260729150859_IDPFI_OPE_20260729150619_20260729150707_N_JLS_O_0091_0027.nc
Copied: W_XX-EUMETSAT-Darmstadt,IMG+SAT,MTI1+FCI-1C-RRAD-FDHSI-FD--CHK-BODY--DIS-NC4E_C_EUMT_20260729150915_IDPFI_OPE_20260729150629_20260729150726_N_JLS_O_0091_0028.nc
Copied: W_XX-EUMETSAT-Darmstadt,IMG+SAT,MTI1+FCI-1C-RRAD-FDHSI-FD--CHK-BODY--DIS-NC4E_C_EUMT_20260729150932_IDPFI_OPE_20260729150648_20260729150744_N_JLS_O_0091_0029.nc
Copied: W_XX-EUMETSAT-Darmstadt,IMG+SAT,MTI1+FCI-1C-RRAD-FDHSI-FD--CHK-BODY--DIS-NC4E_C_EUMT_20260729150949_IDPFI_OPE_20260729150707_20260729150801_N_JLS_O_0091_0030.nc
Copied: W_XX-EUMETSAT-Darmstadt,IMG+SAT,MTI1+FCI-1C-RRAD-FDHSI-FD--CHK-BODY--DIS-NC4E_C_EUMT_20260729151005_IDPFI_OPE_20260729150716_20260729150818_N_JLS_O_0091_0031.nc
Copied: W_XX-EUMETSAT-Darmstadt,IMG+SAT,MTI1+FCI-1C-RRAD-FDHSI-FD--CHK-BODY--DIS-NC4E_C_EUMT_20260729151007_IDPFI_OPE_20260729150726_20260729150826_N_JLS_O_0091_0032.nc
Copied: W_XX-EUMETSAT-Darmstadt,IMG+SAT,MTI1+FCI-1C-RRAD-FDHSI-FD--CHK-BODY--DIS-NC4E_C_EUMT_20260729151018_IDPFI_OPE_20260729150744_20260729150834_N_JLS_O_0091_0033.nc
Copied: W_XX-EUMETSAT-Darmstadt,IMG+SAT,MTI1+FCI-1C-RRAD-FDHSI-FD--CHK-BODY--DIS-NC4E_C_EUMT_20260729151031_IDPFI_OPE_20260729150801_20260729150850_N_JLS_O_0091_0034.nc
Copied: W_XX-EUMETSAT-Darmstadt,IMG+SAT,MTI1+FCI-1C-RRAD-FDHSI-FD--CHK-BODY--DIS-NC4E_C_EUMT_20260729151037_IDPFI_OPE_20260729150810_20260729150904_N_JLS_O_0091_0035.nc
Copied: W_XX-EUMETSAT-Darmstadt,IMG+SAT,MTI1+FCI-1C-RRAD-FDHSI-FD--CHK-BODY--DIS-NC4E_C_EUMT_20260729151046_IDPFI_OPE_20260729150826_20260729150911_N_JLS_O_0091_0036.nc
Copied: W_XX-EUMETSAT-Darmstadt,IMG+SAT,MTI1+FCI-1C-RRAD-FDHSI-FD--CHK-BODY--DIS-NC4E_C_EUMT_20260729151048_IDPFI_OPE_20260729150842_20260729150917_N_JLS_O_0091_0037.nc
Copied: W_XX-EUMETSAT-Darmstadt,IMG+SAT,MTI1+FCI-1C-RRAD-FDHSI-FD--CHK-BODY--DIS-NC4E_C_EUMT_20260729151058_IDPFI_OPE_20260729150850_20260729150929_N_JLS_O_0091_0038.nc
Copied: W_XX-EUMETSAT-Darmstadt,IMG+SAT,MTI1+FCI-1C-RRAD-FDHSI-FD--CHK-BODY--DIS-NC4E_C_EUMT_20260729151059_IDPFI_OPE_20260729150904_20260729150935_N_JLS_O_0091_0039.nc
Copied: W_XX-EUMETSAT-Darmstadt,IMG+SAT,MTI1+FCI-1C-RRAD-FDHSI-FD--CHK-BODY--DIS-NC4E_C_EUMT_20260729151057_IDPFI_OPE_20260729150924_20260729150935_N_JLS_O_0091_0040.nc
In [9]:
#set the environment of satpy
files = find_files_and_readers(base_dir=dst_dir, reader='fci_l1c_nc')
In [10]:
#add the MTG image in satpy memory
scn_tc = Scene(filenames=files)
In [11]:
#print(scn.available_dataset_names()) #uncomment if you want to see the content of the imported image in satpy

First create a natural color composite¶

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!

In [12]:
# scn.load(['true_color'], upper_right_corner='NE')
# #scn.show('natural_color') #uncomment if you want to see the full MTG (daytime) image
# nc = scn['true_color'].data
In [13]:
#load data for selected color visualization
image = 'true_color'
scn_tc.load([image], upper_right_corner='NE')
In [14]:
plt.figure(figsize=(8,8))
img = get_enhanced_image(scn_tc['true_color'])
# get DataArray out of `XRImage` object
img_data = img.data
#comment the line below if you don't want to see the full MTG (daytime) image
img_data.plot.imshow(rgb='bands', vmin=0, vmax=1); #uncomment if you want to see the full MTG (daytime) image
No description has been provided for this image
In [15]:
#rotate the array and create a 1D data array
data1_rot0 = np.rot90(img_data[0], k=4, axes=(1, 0))
data1_rot1 = np.rot90(img_data[1], k=4, axes=(1, 0))
data1_rot2 = np.rot90(img_data[2], k=4, axes=(1, 0))

data1_rot = np.array([data1_rot2, data1_rot1, data1_rot0]).flatten()
In [16]:
#create an empty raster
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')
dfNum = ilwis.DataDefinition(ilwis.NumericDomain('code=value'), ilwis.NumericRange(-100.0, 500.0, 0.001))
rcNew = ilwis.RasterCoverage()
rcNew.setSize(ilwis.Size(11136,11136,3))
rcNew.setGeoReference(grf)
rcNew.setDataDef(dfNum)
In [17]:
#add the data to the ilwis raster
rcNew.array2raster(data1_rot)
In [18]:
#remove if (satpy) nodata
rcNew1 = ilwis.do('mapcalc', 'iff(@1>=0,@1,0)', rcNew)
In [19]:
#store the results as an ilwis maplist - display the map using the ilwis386 desktop software
rcNew1.store('finalTC.mpl')
In [20]:
#define AoI - here Mediterranean region
#note pixelsize Eastings total = 50 degree (-11 to 39) * 110 km per degree / No columns (9001) = roughly 600 meters/pixel
grf_LL= ilwis.GeoReference('code=georef:type=corners, csy=epsg:4326, envelope=-11 55 39 34, gridsize=9001 3781, cornerofcorners=yes')
In [21]:
# Resample the map list and store the results
rcNew_med = ilwis.do("resample", rcNew1, grf_LL, "nearestneighbour")
In [22]:
#stretch the rasterbands contained in rcNew_med using a loop
multiple_stretch = []
multiple_bands = ilwis.do('selection',rcNew_med,"rasterbands(0..2)") 
ls = ilwis.do('linearstretch',multiple_bands, 1) #using an upper and lower data limit defined by the cumulative 1 and 99 % thresholds
mb_stretch = ilwis.do('setvaluerange', ls, 0, 255, 1) #set ouput to byte range
In [23]:
#store the results as an ILWIS maplist - display the map using the ilwis386 desktop software
mb_stretch.store('med_'+str(date_str)+'00.mpl')
In [24]:
#load the individual spectral channels
Blues = ilwis.do('selection',mb_stretch,"rasterbands(0)")
Greens = ilwis.do('selection',mb_stretch,"rasterbands(1)")
Reds = ilwis.do('selection',mb_stretch,"rasterbands(2)")
In [25]:
#transform the spectral channels from iwlsi format to a numpy array using the iterator
Blues_2np = np.fromiter(iter(Blues), np.ubyte, Blues.size().linearSize()) 
Blues_2np = Blues_2np.reshape((Blues.size().ysize, Blues.size().xsize))

Greens_2np = np.fromiter(iter(Greens), np.ubyte, Blues.size().linearSize()) 
Greens_2np = Greens_2np.reshape((Blues.size().ysize, Blues.size().xsize))

Reds_2np = np.fromiter(iter(Reds), np.ubyte, Blues.size().linearSize()) 
Reds_2np = Reds_2np.reshape((Blues.size().ysize, Blues.size().xsize))
In [26]:
#create a numpy 3D data stack
ncol = np.dstack((Reds_2np, Greens_2np, Blues_2np))
In [27]:
#plot the result using matplotlib / cartopy
#note the features, like clouds and smoke plumes (Iberian Peninsula, Greece - Crete and along the south-west coast of Türkiye)
img_extent = (-11, 39, 34, 55) 
fig = plt.figure(figsize=(12, 10))
ax = plt.axes(projection=ccrs.PlateCarree())
plt.title('Meteosat Third Generation composite of '+ (date_str)+'00')

#add Natural earth shape files
ax.add_feature(cfeature.BORDERS, linewidth=0.5, color='yellow', zorder=3)
ax.add_feature(cfeature.COASTLINE, linewidth=0.5, color='blue', zorder=3)

#data raster 
ax.imshow(ncol, origin='upper', extent=img_extent, transform=ccrs.PlateCarree())

gl = ax.gridlines(draw_labels=True,color = 'grey', linestyle = '--', linewidth = 0.5)
gl.top_labels = False
gl.right_labels = False
gl.ylocator = mticker.FixedLocator([35, 40, 45, 50])

plt.show()
#fig.savefig(dst_dir+'/mtg_tc'+(date_str)+'00.jpg')
No description has been provided for this image

Create the MTG Fire Temperature visualization¶

Process the ir_38 channel (Red)¶

In [29]:
scn_tc.load(['ir_38'])
#uncomment the line below to show the spectral channel
#scn_tc.show('ir_38')
In [31]:
ir_38_values = scn_tc['ir_38'].values
# uncomment the line below to retrieve the value for a specific location in the array 
#print(ir_38_values[3000,3000])
In [32]:
#Rotate / flip the image
ir_38_rot = np.flipud(ir_38_values)
In [33]:
#Create a 1D data list
data38 = np.array([ir_38_rot]).flatten()
data38.shape
Out[33]:
(31002624,)
In [34]:
#create an empty raster - note TIR = 5568 lines and columns - 2 km SSP
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')
dfNum = ilwis.DataDefinition(ilwis.NumericDomain('code=value'), ilwis.NumericRange(0.0, 500.0, 0.001))
rcNew = ilwis.RasterCoverage()
rcNew.setSize(ilwis.Size(5568,5568,1))
rcNew.setGeoReference(grf)
rcNew.setDataDef(dfNum)
In [35]:
rcNew.array2raster(data38)
In [36]:
#resmaple to AoI  - here Mediterranean region
rc38_res = ilwis.do('resample', rcNew, grf_LL, 'nearestneighbour')
rc38_res = ilwis.do('mapcalc', 'iff(@1>=0,@1,0)', rc38_res) # remove (satpy) nodata
print(rc38_res.size())

Size(9001, 3781, 1)

Process the nir_22 channel (Green)¶

In [38]:
scn_tc.load(["nir_22"])
# you can access the values of a dataset as a Numpy array with
nir_22_values = scn_tc['nir_22'].values
# uncomment to retrieve the value for a specific location in the array 
#print(nir_22_values[3000,3000]) #uncomment to retrieve the value for the specified location
#uncomment to show the image
#scn_tc.show('nir_22')
In [39]:
nir_22_rot = np.flipud(nir_22_values)
In [40]:
data22 = np.array([nir_22_rot]).flatten()
In [41]:
#create an empty raster - note NIR = 11136 lines and columns - 1 km SSP
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')
dfNum = ilwis.DataDefinition(ilwis.NumericDomain('code=value'), ilwis.NumericRange(-100.0, 500.0, 0.001))
rcNew = ilwis.RasterCoverage()
rcNew.setSize(ilwis.Size(11136,11136,1))
rcNew.setGeoReference(grf)
rcNew.setDataDef(dfNum)
In [42]:
rcNew.array2raster(data22)
In [43]:
rc22_res = ilwis.do('resample', rcNew, grf_LL, 'nearestneighbour')
rc22_res = ilwis.do('mapcalc', 'iff(@1>=0,@1,0)', rc22_res) # remove (satpy) nodata
print(rc22_res.size())

Size(9001, 3781, 1)

Process the nir_16 channel (Blue)¶

In [44]:
scn_tc.load(["nir_16"])
# you can access the values of a dataset as a Numpy array with
nir_16_values = scn_tc['nir_16'].values
# retrieve the value for a specific location in the array 
#uncomment to print the pixel value
#(nir_16_values[3000,3000])
#uncomment to shwo the image
#scn_tc.show('nir_16')
In [45]:
nir_16_rot = np.flipud(nir_16_values)
In [46]:
data16 = np.array([nir_16_rot]).flatten()
In [47]:
rcNew.array2raster(data16)
In [48]:
rc16_res = ilwis.do('resample', rcNew, grf_LL, 'nearestneighbour')
rc16_res = ilwis.do('mapcalc', 'iff(@1>=0,@1,0)', rc16_res) # remove (satpy) nodata
print(rc16_res.size())

Size(9001, 3781, 1)

Image enhancement and creating a color composite¶

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

In [49]:
# Dictionary containing the input raster maps
rasters = {
    "rc38s_2np": rc38_res,
    "rc22s_2np": rc22_res,
    "rc16s_2np": rc16_res,
}

# Output dictionary
rasters_np = {}

for name, raster in rasters.items():

    # Create histogram
    hist = raster.statistics(ilwis.PropertySets.pHISTOGRAM, 65535)

    # Calculate 1% stretch
    minPerc, maxPerc = hist.calcStretchRange(1)

    # Stretch and set value range
    stretched = ilwis.do('linearstretch', raster, minPerc, maxPerc)
    stretched = ilwis.do('setvaluerange', stretched, 0, 255, 1)

    # Convert to NumPy array
    arr = np.fromiter(iter(stretched), np.ubyte, stretched.size().linearSize())
    arr = arr.reshape((stretched.size().ysize, stretched.size().xsize))

    # Store result
    rasters_np[name] = arr
In [50]:
#create a 3D stack for visualization using matplotlib / cartopy
rgb = np.dstack((
    rasters_np["rc38s_2np"],   # R
    rasters_np["rc22s_2np"],   # G
    rasters_np["rc16s_2np"]    # B
))
In [51]:
img_extent = (-11, 39, 34, 55) 
fig = plt.figure(figsize=(12, 10))
ax = plt.axes(projection=ccrs.PlateCarree())
plt.title('Meteosat Third Generation Fire Temperature composite of '+ (date_str)+'00')

#add Natural earth shape files
ax.add_feature(cfeature.BORDERS, linewidth=0.5, color='yellow', zorder=3)
ax.add_feature(cfeature.COASTLINE, linewidth=0.5, color='blue', zorder=3)

#data raster 
ax.imshow(rgb, origin='upper', extent=img_extent, transform=ccrs.PlateCarree())

gl = ax.gridlines(draw_labels=True,color = 'grey', linestyle = '--', linewidth = 0.5)
gl.top_labels = False
gl.right_labels = False
gl.ylocator = mticker.FixedLocator([35, 40, 45, 50])

plt.show()
#uncomment line below to save as a JPG image
#fig.savefig(dst_dir+'/mtg_FT'+(date_str)+'00.jpg')
No description has been provided for this image
In [52]:
#create a 1D list
R = rasters_np["rc38s_2np"]
G = rasters_np["rc22s_2np"]
B = rasters_np["rc16s_2np"]

data = np.array([B, G, R]).flatten()
data.shape
Out[52]:
(102098343,)
In [53]:
#create empty raster 
rcAll_ilw = ilwis.RasterCoverage()
dfNumrc = ilwis.DataDefinition(ilwis.NumericDomain('code=value'), ilwis.NumericRange(0, 255, 0))
rcAll_ilw.setDataDef(dfNumrc)
rcAll_ilw.setSize(ilwis.Size(9001, 3781, 3))
rcAll_ilw.setGeoReference(rc38_res.geoReference())
In [54]:
#add the data to the raster
rcAll_ilw.array2raster(data) 
print(rcAll_ilw.size())

Size(9001, 3781, 3)
In [55]:
#store the results as an ILWIS maplist - display the map using the ilwis386 desktop software
rcAll_ilw.store('FT_comp_'+date_str+'.mpl')

Display the results for some selected areas¶

In [56]:
# Full image extent (same for both images)
img_extent = (-11, 39, 34, 55)

# Submap extent: lon_min, lon_max, lat_min, lat_max
sub_extent = (-9.55, -4.10, 40.50, 43.79)

fig, axes = plt.subplots(
    1, 2,
    figsize=(16, 8),
    subplot_kw={"projection": ccrs.PlateCarree()}
)

titles = [
    f"MTG Fire Temperature composite {date_str}00",
    f"MTG True Color composite {date_str}00"
]

images = [rgb, ncol]

for ax, img, title in zip(axes, images, titles):

    # Zoom to submap
    ax.set_extent(sub_extent, crs=ccrs.PlateCarree())

    # Add borders and coastlines
    ax.add_feature(cfeature.BORDERS, linewidth=0.5,
                   color='yellow', zorder=3)
    ax.add_feature(cfeature.COASTLINE, linewidth=0.5,
                   color='blue', zorder=3)

    # Plot image
    ax.imshow(
        img,
        origin='upper',
        extent=img_extent,
        transform=ccrs.PlateCarree()
    )

    # Gridlines
    gl = ax.gridlines(
        draw_labels=True,
        color='grey',
        linestyle='--',
        linewidth=0.5
    )

    gl.top_labels = False
    gl.right_labels = False

    gl.xlocator = mticker.FixedLocator([-8,-7, -6, - 5, -4])
    gl.ylocator = mticker.FixedLocator([41, 42, 43])

    ax.set_title(title)

plt.tight_layout()
plt.show()
No description has been provided for this image
In [57]:
# Full image extent (same for both images)
img_extent = (-11, 39, 34, 55)

# Submap extent: lon_min, lon_max, lat_min, lat_max
sub_extent = (23.34, 30.12, 34.13, 40.32)

fig, axes = plt.subplots(
    1, 2,
    figsize=(16, 8),
    subplot_kw={"projection": ccrs.PlateCarree()}
)

titles = [
    f"MTG Fire Temperature composite {date_str}00",
    f"MTG True Color composite {date_str}00"
]

images = [rgb, ncol]

for ax, img, title in zip(axes, images, titles):

    # Zoom to submap
    ax.set_extent(sub_extent, crs=ccrs.PlateCarree())

    # Add borders and coastlines
    ax.add_feature(cfeature.BORDERS, linewidth=0.5,
                   color='yellow', zorder=3)
    ax.add_feature(cfeature.COASTLINE, linewidth=0.5,
                   color='blue', zorder=3)

    # Plot image
    ax.imshow(
        img,
        origin='upper',
        extent=img_extent,
        transform=ccrs.PlateCarree()
    )

    # Gridlines
    gl = ax.gridlines(
        draw_labels=True,
        color='grey',
        linestyle='--',
        linewidth=0.5
    )

    gl.top_labels = False
    gl.right_labels = False

    gl.xlocator = mticker.FixedLocator([24, 26, 28, 30])
    gl.ylocator = mticker.FixedLocator([36, 38, 40])

    ax.set_title(title)

plt.tight_layout()
plt.show()
No description has been provided for this image

Process Fire Radiative Power¶

See also: https://data.eumetsat.int/product/EO:EUM:DAT:1156

file name convention = W_PT-LSASAF-LISBON,SATELLITE,LSA-509_MTG_MTFRPPIXEL_MTG-FD_C_LPMG_20260729150000.nc

In [58]:
#check timestamp
date_str
Out[58]:
'2026072915'
In [59]:
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')):
        shutil.copy(filename, dst_dir)
In [60]:
fn = dst_dir+'/W_PT-LSASAF-LISBON,SATELLITE,LSA-509_MTG_MTFRPPIXEL_MTG-FD_C_LPMG_'+(date_str)+'0000.nc'
ds = nc.Dataset(fn)

print(ds.groups.keys())
print(ds.groups["ListProduct"].variables.keys())

dict_keys(['ListProduct', 'QualityProduct'])
dict_keys(['BW_SIZE', 'BW_NUMPIX', 'BW_BT_MIR', 'BW_BTD', 'RAD_BCK', 'STD_BCK', 'FIRE_CONFIDENCE', 'BT_MIR', 'BT_TIR', 'RAD_PIX', 'PIXEL_VZA', 'PIXEL_SZA', 'PIXEL_SIZE', 'ACQTIME', 'ABS_LINE', 'ABS_SAMP', 'LONGITUDE', 'LATITUDE', 'LONGITUDE_PARALLAX', 'LATITUDE_PARALLAX', 'FRP', 'ERR_FRP_COEFF', 'ERR_ATM_TRANS', 'ERR_RADIOMETRIC', 'ERR_BACKGROUND', 'ERR_VERT_COMP', 'FRP_UNCERTAINTY', 'PIXEL_ATM_TRANS', 'index'])
In [61]:
# Open the ListProduct group
grp = ds.groups["ListProduct"]

# Create DataFrame
df = pd.DataFrame({
    "latitude": grp.variables["LATITUDE"][:],
    "longitude": grp.variables["LONGITUDE"][:],
    "fire_radiative_power": grp.variables["FRP"][:],
    "fire_confidence": grp.variables["FIRE_CONFIDENCE"][:],
    "frp_uncertainty": grp.variables["FRP_UNCERTAINTY"][:]
})

print(df.head())
#print(df.info())

    latitude  longitude  fire_radiative_power  fire_confidence  \
0  51.174240   3.807763             17.378065         0.902039   
1  49.691208  18.666277             22.724171         0.622751   
2  49.692051  18.682489             15.127853         0.216106   
3  49.673248  18.658327             15.058973         0.752949   
4  49.674091  18.674532             14.384563         0.751470   

   frp_uncertainty  
0         5.782651  
1         5.070397  
2         4.323100  
3         4.150695  
4         4.143226
In [62]:
len(df)
Out[62]:
7954
In [63]:
# Create point geometry from longitude and latitude
geometry = gpd.points_from_xy(
    df["longitude"],
    df["latitude"]
)

# Create GeoDataFrame
gdf = gpd.GeoDataFrame(
    df,
    geometry=geometry,
    crs="EPSG:4326"     # WGS84 latitude/longitude
)

# Save as ESRI Shapefile
gdf.to_file(dst_dir+"/MTG_FRP.shp", driver="ESRI Shapefile")

print(f"Created shapefile with {len(gdf)} fire detections.")

Created shapefile with 7954 fire detections.
In [64]:
# Filter to area of interest
df_aoi = df[
    (df["latitude"] >= 34.0) &
    (df["latitude"] <= 55.0) &
    (df["longitude"] >= -11.0) &
    (df["longitude"] <= 39.0)
].copy()

print(f"Number of fire pixels in AOI: {len(df_aoi)}")
print(df_aoi.head())

Number of fire pixels in AOI: 272
    latitude  longitude  fire_radiative_power  fire_confidence  \
0  51.174240   3.807763             17.378065         0.902039   
1  49.691208  18.666277             22.724171         0.622751   
2  49.692051  18.682489             15.127853         0.216106   
3  49.673248  18.658327             15.058973         0.752949   
4  49.674091  18.674532             14.384563         0.751470   

   frp_uncertainty  
0         5.782651  
1         5.070397  
2         4.323100  
3         4.150695  
4         4.143226
In [65]:
# Create point geometry from longitude and latitude
geometry = gpd.points_from_xy(
    df_aoi["longitude"],
    df_aoi["latitude"]
)

# Create GeoDataFrame
gdf = gpd.GeoDataFrame(
    df_aoi,
    geometry=geometry,
    crs="EPSG:4326"     # WGS84 latitude/longitude
)

# Save as ESRI Shapefile
gdf.to_file(dst_dir+"/MTG_FRP_AOI.shp", driver="ESRI Shapefile")

print(f"Created shapefile with {len(gdf)} fire detections.")

Created shapefile with 272 fire detections.

Plot final results¶

In [66]:
fig, ax = plt.subplots(
    figsize=(12, 10),
    subplot_kw={"projection": ccrs.PlateCarree()}
)

extent = (-11, 39, 34, 55)

# RGB image
ax.imshow(
    rgb,
    extent=extent,
    origin="upper",
    transform=ccrs.PlateCarree()
)

# Fire detections
norm = LogNorm(
    vmin=max(1, df_aoi["fire_radiative_power"].min()),  # avoid zero for LogNorm
    vmax=df_aoi["fire_radiative_power"].max()
)

sc = ax.scatter(
    df_aoi["longitude"],
    df_aoi["latitude"],
    c=df_aoi["fire_radiative_power"],
    cmap="hot",
    norm=norm,
    s=40,
    edgecolors="black",
    linewidth=0.3,
    transform=ccrs.PlateCarree(),
    zorder=10
)

# Coastlines and borders
ax.add_feature(cfeature.COASTLINE, linewidth=0.5)
ax.add_feature(cfeature.BORDERS, linewidth=0.5)

# Colorbar
cbar = plt.colorbar(sc, ax=ax, shrink=0.3, extend="max")
cbar.set_label("Fire Radiative Power (MW)")

# Map extent
ax.set_extent(extent, crs=ccrs.PlateCarree())

# Gridlines
gl = ax.gridlines(
    draw_labels=True,
    color="grey",
    linestyle="--",
    linewidth=0.5
)

gl.top_labels = False
gl.right_labels = False

gl.xlocator = mticker.FixedLocator([-10, -5, 0, 5, 10, 15, 20, 25, 30, 35])
gl.ylocator = mticker.FixedLocator([35, 40, 45, 50, 55])

ax.set_title("MTG Fire Temperature RGB with Fire Radiative Power Detections")

plt.tight_layout()
plt.show()
No description has been provided for this image

Plot some selected areas for better visualization¶

In [67]:
# ------------------------------------------------------------------
# Image extent (full resampled image)
# ------------------------------------------------------------------
img_extent = (-11, 39, 34, 55)

# ------------------------------------------------------------------
# Create figure with two subplots
# ------------------------------------------------------------------
fig, axes = plt.subplots(
    1, 2,
    figsize=(16, 8),
    subplot_kw={"projection": ccrs.PlateCarree()}
)

# Overall title
fig.suptitle(
    f"NW Iberian Peninsula MTG Fire Temperature RGB with Fire Radiative Power Detections\nMTG FCI - {date_str}00 UTC",
    fontsize=18,
    fontweight="bold"
)

# Plot settings
titles = [
    "Fire Temperature RGB",
    "Natural Color RGB"
]

images = [
    rgb,
    ncol
]

# Keep only positive FRP values for logarithmic colouring
df_plot = df_aoi[df_aoi["fire_radiative_power"] > 0]

# ------------------------------------------------------------------
# Plot both maps
# ------------------------------------------------------------------
for ax, image, title in zip(axes, images, titles):

    # Zoom to Iberian Peninsula
    ax.set_extent(
        [-9.00, -5.20, 40.50, 43.79],
        crs=ccrs.PlateCarree()
    )
    
    # Background RGB image
    ax.imshow(
        image,
        origin="upper",
        extent=img_extent,
        transform=ccrs.PlateCarree()
    )

    # Coastlines and borders
    ax.add_feature(
        cfeature.COASTLINE,
        linewidth=0.5,
        color="blue",
        zorder=3
    )

    ax.add_feature(
        cfeature.BORDERS,
        linewidth=0.5,
        color="yellow",
        zorder=3
    )

    # Fire detections
    sc = ax.scatter(
        df_plot["longitude"],
        df_plot["latitude"],
        c=df_plot["fire_radiative_power"],
        cmap="hot",
        norm=LogNorm(
            vmin=df_plot["fire_radiative_power"].min(),
            vmax=df_plot["fire_radiative_power"].max()
        ),
        s=35,
        edgecolors="black",
        linewidth=0.3,
        transform=ccrs.PlateCarree(),
        zorder=10
    )

    # Gridlines
    gl = ax.gridlines(
        draw_labels=True,
        color="grey",
        linestyle="--",
        linewidth=0.5
    )

    gl.top_labels = False
    gl.right_labels = False

    gl.xlocator = mticker.FixedLocator(
        [-9, -8, -7, -6, -5, -4, -2, 0]
    )

    gl.ylocator = mticker.FixedLocator(
        [37, 38, 39, 40, 41, 42, 43]
    )

    ax.set_title(title, fontsize=13)

# -----------------------------------------------------------
# Create a shared logarithmic colorbar
# -----------------------------------------------------------

# Set logarithmic scaling limits
vmin = 1      # MW (must be > 0)
vmax = df_plot["fire_radiative_power"].max()

# Update the scatter plots to use the same normalization
norm = LogNorm(vmin=vmin, vmax=vmax)

# In the plotting loop, use:
sc = ax.scatter(
    df_plot["longitude"],
    df_plot["latitude"],
    c=df_plot["fire_radiative_power"],
    cmap="hot",
    norm=norm,
    s=35,
    edgecolors="black",
    linewidth=0.3,
    transform=ccrs.PlateCarree(),
    zorder=10
)

# -----------------------------------------------------------
# Create colorbar axis on the right
# -----------------------------------------------------------
cbar = fig.colorbar(
    sc,
    ax=axes,
    location="right",
    shrink=0.55,
    pad=0.03,
    extend = 'max'
)

# Get current position
pos = cbar.ax.get_position()

# Move it upwards
cbar.ax.set_position([
    pos.x0 + 0.14,    # left
    pos.y0 + 0.04,    # move up
    pos.width,        # width
    pos.height        # height
])
cbar.set_label("Fire Radiative Power (MW)", fontsize=11)

# Nice logarithmic tick locations
cbar.set_ticks([1, 3, 10, 30, 100, 300, 1000])

fig.subplots_adjust(
    top=0.96,      # controls top margin
    wspace=0.10   # space between the two maps
)

plt.show()
No description has been provided for this image
In [68]:
# ------------------------------------------------------------------
# Image extent (full resampled image)
# ------------------------------------------------------------------
img_extent = (-11, 39, 34, 55)

# ------------------------------------------------------------------
# Create figure with two subplots
# ------------------------------------------------------------------
fig, axes = plt.subplots(
    1, 2,
    figsize=(16, 8),
    subplot_kw={"projection": ccrs.PlateCarree()}
)

# Overall title
fig.suptitle(
    f"Eastern Mediterranean MTG Fire Temperature RGB with Fire Radiative Power Detections\nMTG FCI - {date_str}00 UTC",
    fontsize=18,
    fontweight="bold"
)

# Plot settings
titles = [
    "Fire Temperature RGB",
    "Natural Color RGB"
]

images = [
    rgb,
    ncol
]

# Keep only positive FRP values for logarithmic colouring
df_plot = df_aoi[df_aoi["fire_radiative_power"] > 0]

# ------------------------------------------------------------------
# Plot both maps
# ------------------------------------------------------------------
for ax, image, title in zip(axes, images, titles):

    # Zoom to Eastern Mediterranean
    ax.set_extent(
        [23.34, 30.12, 34.13, 40.32],
        crs=ccrs.PlateCarree()
    )
    
    # Background RGB image
    ax.imshow(
        image,
        origin="upper",
        extent=img_extent,
        transform=ccrs.PlateCarree()
    )

    # Coastlines and borders
    ax.add_feature(
        cfeature.COASTLINE,
        linewidth=0.5,
        color="blue",
        zorder=3
    )

    ax.add_feature(
        cfeature.BORDERS,
        linewidth=0.5,
        color="yellow",
        zorder=3
    )

    # Fire detections
    sc = ax.scatter(
        df_plot["longitude"],
        df_plot["latitude"],
        c=df_plot["fire_radiative_power"],
        cmap="hot",
        norm=LogNorm(
            vmin=df_plot["fire_radiative_power"].min(),
            vmax=df_plot["fire_radiative_power"].max()
        ),
        s=35,
        edgecolors="black",
        linewidth=0.3,
        transform=ccrs.PlateCarree(),
        zorder=10
    )

    # Gridlines
    gl = ax.gridlines(
        draw_labels=True,
        color="grey",
        linestyle="--",
        linewidth=0.5
    )

    gl.top_labels = False
    gl.right_labels = False

    gl.xlocator = mticker.FixedLocator(
        [24, 25, 26, 27, 28, 29, 30]
    )

    gl.ylocator = mticker.FixedLocator(
        [35, 36, 37, 38, 39, 40]
    )

    ax.set_title(title, fontsize=13)

# -----------------------------------------------------------
# Create a shared logarithmic colorbar
# -----------------------------------------------------------

# Set logarithmic scaling limits
vmin = 1      # MW (must be > 0)
vmax = df_plot["fire_radiative_power"].max()

# Update the scatter plots to use the same normalization
norm = LogNorm(vmin=vmin, vmax=vmax)

# In the plotting loop, use:
sc = ax.scatter(
    df_plot["longitude"],
    df_plot["latitude"],
    c=df_plot["fire_radiative_power"],
    cmap="hot",
    norm=norm,
    s=35,
    edgecolors="black",
    linewidth=0.3,
    transform=ccrs.PlateCarree(),
    zorder=10
)

# -----------------------------------------------------------
# Create colorbar axis on the right
# -----------------------------------------------------------
cbar = fig.colorbar(
    sc,
    ax=axes,
    location="right",
    shrink=0.55,
    pad=0.03,
    extend = 'max'
)

# Get current position
pos = cbar.ax.get_position()

# Move it upwards
cbar.ax.set_position([
    pos.x0 + 0.14,    # left
    pos.y0 + 0.04,    # move up
    pos.width,        # width
    pos.height        # height
])
cbar.set_label("Fire Radiative Power (MW)", fontsize=11)

# Nice logarithmic tick locations
cbar.set_ticks([1, 3, 10, 30, 100, 300, 1000])


fig.subplots_adjust(
    top=0.96,      # controls top margin
    wspace=0.10   # space between the two maps
)

plt.show()
No description has been provided for this image

Remove the EUMETCast MTG segments and RFP from the destination folder¶

New instance of kernel is started as files are still occupied by a previous instance of python

In [69]:
import IPython
IPython.Application.instance().kernel.do_shutdown(True) #automatically restarts kernel
Out[69]:
{'status': 'ok', 'restart': True}
In [1]:
import sys, os 
import glob
In [2]:
#local output folder
dst_dir = os.getcwd()+'/Fire_result'
In [3]:
list = glob.glob(dst_dir+'/W_XX-EUMETSAT-Darmstadt,IMG+SAT,MTI1+FCI-1C-RRAD-FDHSI-FD--CHK-BODY--DIS-NC4E_C_EUMT_*.nc')
list = [elem[len(dst_dir)+1:] for elem in list]
for elem in list:
    os.remove(dst_dir+'/'+elem)
In [4]:
list = glob.glob(dst_dir+'/W_PT-LSASAF-LISBON,SATELLITE,LSA-509_MTG_MTFRPPIXEL_MTG-FD_C_LPMG_*0000.nc')
list = [elem[len(dst_dir)+1:] for elem in list]
for elem in list:
    os.remove(dst_dir+'/'+elem)
In [ ]:
In [ ]: