2
$\begingroup$

I am trying to read an MSG (Meteosat Second Generation) file given in the native(.nat) format using satpy. I tried using:

import numpy as np
import os
from satpy import Scene
import matplotlib.pyplot as plt
import pandas as pd
from pyproj import Proj
from pyproj import transform
file = '/data/swadhin/meteosat/l1.5_data/native/MSG1-SEVI-MSG15-0100-NA-20211201035743.145000000Z-NA.nat'

# define reader
reader = "seviri_l1b_native"
# read the file
scn = Scene(filenames = {reader:[file]})
# extract data set names
dataset_names = scn.all_dataset_names()
# print available datasets
print('\n'.join(map(str, dataset_names)))

which gave me

HRV
IR_016
IR_039
IR_087
IR_097
IR_108
IR_120
IR_134
VIS006
VIS008
WV_062
WV_073

I am now trying to read the 10.8 and 12 micron bands as follows:

import matplotlib.pyplot as plt

scn.load(['IR_108'])

# Access the dataset
ir108 = scn['IR_108']
    
# Print the metadata
print(ir108.attrs)

which again gives me:

{'orbital_parameters': {'projection_longitude': np.float32(41.5), 'projection_latitude': 0.0, 'projection_altitude': 35785831.0, 'satellite_nominal_longitude': np.float32(41.5), 'satellite_nominal_latitude': 0.0, 'satellite_actual_longitude': 41.639685531217864, 'satellite_actual_latitude': 6.875188281971676, 'satellite_actual_altitude': 35779932.76250567}, 'units': 'K', 'wavelength': WavelengthRange(min=9.8, central=10.8, max=11.8, unit='µm'), 'standard_name': 'toa_brightness_temperature', 'platform_name': 'Meteosat-8', 'sensor': 'seviri', 'georef_offset_corrected': np.True_, 'time_parameters': {'nominal_start_time': datetime.datetime(2021, 12, 1, 3, 45), 'nominal_end_time': datetime.datetime(2021, 12, 1, 4, 0), 'observation_start_time': datetime.datetime(2021, 12, 1, 3, 45, 11, 146000), 'observation_end_time': datetime.datetime(2021, 12, 1, 3, 57, 43, 145000)}, 'start_time': datetime.datetime(2021, 12, 1, 3, 45), 'end_time': datetime.datetime(2021, 12, 1, 4, 0), 'reader': 'seviri_l1b_native', 'area': Area ID: msg_seviri_iodc_3km
Description: MSG SEVIRI Indian Ocean Data Coverage service area definition with 3 km resolution
Projection: {'a': '6378169', 'h': '35785831', 'lon_0': '41.5', 'no_defs': 'None', 'proj': 'geos', 'rf': '295.488065897014', 'type': 'crs', 'units': 'm', 'x_0': '0', 'y_0': '0'}
Number of columns: 3712
Number of rows: 3712
Area extent: (np.float32(5567248.0), np.float32(5570248.5), np.float32(-5570248.5), np.float32(-5567248.0)), 'name': 'IR_108', 'resolution': 3000.403165817, 'calibration': 'brightness_temperature', 'modifiers': (), '_satpy_id': DataID(name='IR_108', wavelength=WavelengthRange(min=9.8, central=10.8, max=11.8, unit='µm'), resolution=3000.403165817, calibration=<2>, modifiers=()), 'ancillary_variables': []}

But when I try to get the brightness temperature for a particular lat and long, I tried constructing the grid as:

# Get the brightness temperature data
brightness_temp = ir108.data.compute() 

# Get the area definition (contains geolocation information)
area_def = ir108.attrs['area']

# Projection parameters from ir108.attrs
proj_params = area_def.proj_dict

# Define the projection using pyproj
proj = Proj(proj='geos', h=35785831, lon_0=41.5, a=6378169, rf=295.488065897014, units='m')

# Create a grid of pixel coordinates
x = np.linspace(area_def.area_extent[0], area_def.area_extent[2], area_def.shape[1])
y = np.linspace(area_def.area_extent[1], area_def.area_extent[3], area_def.shape[0])
x_coords, y_coords = np.meshgrid(x, y)

# Transform pixel coordinates to geographic coordinates
lons, lats = proj(x_coords, y_coords, inverse=True)

# Flatten the arrays for DataFrame creation
flat_lons = lons.flatten()
flat_lats = lats.flatten()
flat_brightness_temp = brightness_temp.flatten()

# Get the time information from the metadata
start_time = ir108.attrs['start_time']
end_time = ir108.attrs['end_time']

# Create a DataFrame
df = pd.DataFrame({
    'latitude': flat_lats,
    'longitude': flat_lons,
    'brightness_temperature': flat_brightness_temp,
    'start_time': start_time,
    'end_time': end_time
})

# Display the DataFrame
print(df.head())

which gave me

  latitude  longitude  brightness_temperature          start_time  \
0       0.0  41.500000                     NaN 2021-12-01 03:45:00   
1       0.0  41.500009                     NaN 2021-12-01 03:45:00   
2       0.0  41.500018                     NaN 2021-12-01 03:45:00   
3       0.0  41.500027                     NaN 2021-12-01 03:45:00   
4       0.0  41.500036                     NaN 2021-12-01 03:45:00   

             end_time  
0 2021-12-01 04:00:00  
1 2021-12-01 04:00:00  
2 2021-12-01 04:00:00  
3 2021-12-01 04:00:00  
4 2021-12-01 04:00:00  

Where am I doing it wrong? I am doing it for the first time and can't figure it out.

I am uploading the file here. File link - .nat file

$\endgroup$

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.