Мобильная версияКод: Выделить всё
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) данные с устройства формирования изображения молний Meteosat третьего поколения (MTG LI) в координаты широты и долготы?
Это информация 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
Буду признателен за любые предложения.
Подробнее здесь: https://stackoverflow.com/questions/793 ... -data-plot
