MTG闪电成像仪（LI）L2级累积闪区面积（AFA）绘图无海陆特征及经纬度转换错误的技术求助

问题一：绘制AFA数据时无法显示海陆特征 #

我尝试绘制MTG LI L2级的累积闪区面积（AFA）数据，但遇到了问题——生成的图像里没有显示陆地和海洋特征。我编写的代码如下：

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()

运行这段代码后得到的图像如下，完全看不到海陆特征：

请问该如何修复这个问题？

问题二：AFA数据经纬度转换结果不正确 #

我还需要将MTG LI的AFA数据转换为经纬度坐标，但自己编写的代码得到的结果不正确。以下是NC文件的关键信息：

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

我尝试了下面的转换代码，但得到的经纬度结果不正确：

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

我参考了MTG LI L2数据指南及类似的静止卫星成像仪投影相关资料，但还是没找到问题所在。希望能得到大家的建议！

备注：内容来源于stack exchange，提问作者meteorolog