python - MTG Lightning AFA data plot - Stack Overflow

admin管理员组

文章数量:1124566

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.

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: /resources/user-guides/mtg-li-level-2-data-guide and here .php for GOES Imager Projection (similar to EUMSAT)

I would appreciate any suggestions.

本文标签： pythonMTG Lightning AFA data plotStack Overflow