MTG Lightning AFA data plot

Question

I am attempting to plot the MTG LI level 2 data (AFA), but I am encountering difficulties. The code I have written is as follows:

import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import numpy as np
from pyproj import Proj
# Load the dataset
file_nc = "W_XX-EUMETSAT-Darmstadt,IMG+SAT,MTI1+LI-2-AFA--FD--CHK-BODY--ARC-NC4E_C_EUMT_20241217110019_L2PF_OPE_20241217105000_20241217110000_N__O_0066_0001.nc"
data = xr.open_dataset(file_nc)
# Extract variables
accumulated_flash_area = data['accumulated_flash_area'].values
scale_factor_x = -5.58871526031607E-5
add_offset_x = 0.155617776423501
scale_factor_y = 5.58871526031607E-5
add_offset_y = -0.155617776423501
x_raw = data['x'].values
y_raw = data['y'].values
# Decode x and y with scale and offset
x_coords = x_raw * scale_factor_x + add_offset_x
y_coords = y_raw * scale_factor_y + add_offset_y
# Define geostationary projection
proj_metadata = data['mtg_geos_projection'].attrs
geos_proj = ccrs.Geostationary(
 central_longitude=proj_metadata.get('longitude_of_projection_origin', 0.0),
 satellite_height=proj_metadata.get('perspective_point_height', 35786000.0)
)
# Initialize the plot
fig, ax = plt.subplots(subplot_kw={'projection': geos_proj}, figsize=(12, 8))
# Add map features
ax.add_feature(cfeature.LAND, facecolor='lightgrey')
ax.add_feature(cfeature.OCEAN, facecolor='lightblue')
# Plot the filtered data
sc = ax.scatter(
 x_coords, y_coords, c=accumulated_flash_area, 
 s=5, cmap='viridis', 
 transform=geos_proj, 
 alpha=0.8
)
# Add a colorbar
plt.colorbar(sc, ax=ax, orientation='vertical', label='Accumulated Flash Area')
# Add a title
ax.set_title("Accumulated Flash Area for Europe (Geostationary Projection)", fontsize=14)
# Show the plot
plt.show()

In this way, we will get plots like this, and the plots will be without land and ocean.

Accumulated Flash Area(AFA)

How to fix this?

Second question is how to convert Accumulated Flash Area (AFA) data from the Meteosat Third Generation Lightning Imager (MTG LI) into latitude and longitude coordinates ? This is nc info :

dimensions:
 enumtype_dim = 1;
 auxiliary_dataset = 1;
 pixels = 577609;
 accumulations = 20;
 variables:
 String auxiliary_dataset_identifier(auxiliary_dataset=1);
 :long_name = "Auxiliary dataset identifier";
 :title = "Identifier of auxiliary dataset or type of auxiliary dataset used to produce this product";
 enum auxiliary_dataset_status_type auxiliary_dataset_status(auxiliary_dataset=1);
 :long_name = "Status of auxiliary dataset";
 :meaning = "0 = OK, 1 = out_of_validity_time, 2 = not_available";
 int mtg_geos_projection;
 :long_name = "MTG geostationary projection";
 :grid_mapping_name = "geostationary";
 :perspective_point_height = 3.57864E7; // double
 :semi_major_axis = 6378137.0; // double
 :semi_minor_axis = 6356752.31424518; // double
 :inverse_flattening = 298.257223563; // double
 :latitude_of_projection_origin = 0.0; // double
 :longitude_of_projection_origin = 0.0; // double
 :sweep_angle_axis = "y";
 double accumulation_start_times(accumulations=20);
 :long_name = "Accumulation start time";
 :units = "seconds since 2000年01月01日 00:00:00.0";
 :precision = "0.001";
 uint accumulation_offsets(accumulations=20);
 :_Unsigned = "true";
 short x(pixels=577609);
 :long_name = "azimuth angle encoded as column";
 :axis = "X";
 :units = "radian";
 :standard_name = "projection_x_coordinate";
 :scale_factor = -5.58871526031607E-5; // double
 :add_offset = 0.155617776423501; // double
 :valid_range = 1S, 5568S; // short
 :_ChunkSizes = 131072U; // uint
 short y(pixels=577609);
 :long_name = "elevation angle encoded as row";
 :axis = "Y";
 :units = "radian";
 :standard_name = "projection_y_coordinate";
 :scale_factor = 5.58871526031607E-5; // double
 :add_offset = -0.155617776423501; // double
 :valid_range = 1S, 5568S; // short
 :_ChunkSizes = 131072U; // uint
 uint accumulated_flash_area(pixels=577609);
 :_Unsigned = "true";
 :long_name = "Number of contributing unique flashes to each pixel";
 :grid_mapping = "mtg_geos_projection";
 :coordinate = "sparse: x y";
 :_ChunkSizes = 65536U; // uint
 byte l1b_missing_warning(accumulations=20);
 :long_name = "Expected L1b inputs missing";
 byte l1b_geolocation_warning(accumulations=20);
 :long_name = "L1b event geolocation warning";
 byte l1b_radiometric_warning(accumulations=20);
 :long_name = "L1b event radiometric warning";
 ubyte average_flash_qa(accumulations=20);
 :_Unsigned = "true";
 :_FillValue = 255UB; // ubyte
 :long_name = "average flash confidence value";
 :add_offset = 0.0f; // float
 :scale_factor = 0.004f; // float

I did this:

import numpy as np
from netCDF4 import Dataset
def calculate_lat_lon(file_path):
 # Open the MTG LI data file
 with Dataset(file_path, 'r') as nc_file:
 # Read projection variables and constants
 x = nc_file.variables['x'][:] # E/W scanning angle in radians
 y = nc_file.variables['y'][:] # N/S elevation angle in radians
 proj_info = nc_file.variables['goes_imager_projection']
 lon_origin = proj_info.longitude_of_projection_origin
 H = proj_info.perspective_point_height + proj_info.semi_major_axis
 r_eq = proj_info.semi_major_axis
 r_pol = proj_info.semi_minor_axis
 # Create 2D coordinate matrices from 1D coordinate vectors
 x_2d, y_2d = np.meshgrid(x, y)
 # Calculate intermediate variables
 lambda_0 = np.deg2rad(lon_origin)
 a = np.sin(x_2d)**2 + (np.cos(x_2d)**2) * (np.cos(y_2d)**2 + ((r_eq**2 / r_pol**2) * np.sin(y_2d)**2))
 b = -2 * H * np.cos(x_2d) * np.cos(y_2d)
 c = H**2 - r_eq**2
 # Calculate the distance from the satellite to each point
 r_s = (-b - np.sqrt(b**2 - 4 * a * c)) / (2 * a)
 # Calculate the geocentric coordinates
 s_x = r_s * np.cos(x_2d) * np.cos(y_2d)
 s_y = -r_s * np.sin(x_2d)
 s_z = r_s * np.cos(x_2d) * np.sin(y_2d)
 # Calculate latitude and longitude
 lat = np.rad2deg(np.arctan((r_eq**2 / r_pol**2) * (s_z / np.sqrt((H - s_x)**2 + s_y**2))))
 lon = np.rad2deg(lambda_0 - np.arctan(s_y / (H - s_x)))
 return lat, lon

However, the lat and long are incorrect. Can you tell me what I'm doing wrong? Here is more information about projection: https://user.eumetsat.int/resources/user-guides/mtg-li-level-2-data-guide and here https://www.star.nesdis.noaa.gov/atmospheric-composition-training/python_abi_lat_lon.php for GOES Imager Projection (similar to EUMSAT)

I would appreciate any suggestions.

• ThanasisN: Do you had any progress? I have a similar problem with MTLST

Answer 1

I am also trying to work with this data.

As far as I managed to progress, to get x,y (pixel) coordinates that range from [0,5568) I used the following mapping:

 lambda_s = np.array(nc.variables['x'][:]) 
 lambda_0 = nc.variables['x'].add_offset
 ags = nc.variables['x'].scale_factor
 f_s = np.array(nc.variables['y'][:])
 f_0 = nc.variables['y'].add_offset
 egs = nc.variables['y'].scale_factor
 c = (lambda_s - lambda_0)/ags #+ 1
 r = -(f_s - f_0)/egs #+ 1
 c = np.rint(c)
 r = np.rint(r)
 
 li_map = np.zeros((5568,5568))
 index_map = np.zeros_like(li_map)
 # print(self.type,type(self.type))
 
 try:
 accumulated_data = nc.variables[self.type ][:] * nc.variables[self.type].scale_factor
 except:
 accumulated_data = nc.variables[self.type ][:]
 indecies = np.array(list(zip(np.array(c,dtype=int),np.array(r,dtype=int))))
 for k,(az,el) in enumerate(indecies):
 li_map [el,az] += accumulated_data[k]

I am also stuck in transforming(projecting) the 5568x5568 image(matrix), which is in Geostationary coordinates, into Plate Carree (Geographical). My final goal is to align it with the Geographical images from FCI (2km grid) images.

It would be interesting if someone managed this. I Hope this helps.

PS, maybe this also will be helpful: Projecting GOES-16 Geostationary data into Plate Carree Cartopy