Add support for MTG Satellite imagery

[jacobbieker]

Detailed Description

The MTG satellite is going live soonish ((https://www.eumetsat.int/getting-ready-data-first-mtg-missions) and will slowly replace the current MSG full-disk imager over the course of 2024.

Context

The full-disk imaging is useful as the backup to the RSS imaging that stops every 28 days for 48 hours, and 1 month each year. So this is important for keeping the forecasts running correctly. The MTG data is also much higher resolution, and every 10 minutes, so could be worth switching to it full time.

Possible Implementation

Satip has support for the FCI imager https://satpy.readthedocs.io/en/stable/api/satpy.readers.fci_l1c_nc.html with some notes of what needs to be installed for it. A decoder needs to be installed as well.

[sanjayk0508]

Hi @jacobbieker, Could you please tell me the approach for adding support for MTG Satellite imagery? EUMETSAT has released the early November MTG data. I want to help out by working on this good first issue. Hope you don't mind me asking a beginner question.

[jacobbieker]

Yeah, no problem! Happy to have you help! The best place to put the processing code would probably be in satip/utils.py. The approach for adding support would be to download one of the raw native files from the MTG data, and adapt the convert_scene_to_dataarray or load_native_to_dataarray functions to support all the channels that MTG has. This would probably include adding either an argument to the function to tell whether its MTG or MSG data, or use the metadata in the native files to pick that automatically itself.

The other part of the support would be to enable Satip to download MTG data from EUMETSAT directly, but not sure if that is viable right now since I don't think they are disseminating it until the new year on an API. Feel free to ask any more questions, happy to help!

[bikramb98]

@jacobbieker I tried looking for MTG data in the EUMETSAT User Portal, but couldn't find any. Can you point me to where I can find a sample data from MTG, and I can get started on developing this feature?

[jacobbieker]

Yes, here is some example data: https://sftp.eumetsat.int/public/folder/UsCVknVOOkSyCdgpMimJNQ/User-Materials/Test-Data/MTG/MTG_FCI_L1C_RC72-20240113_TD-505_Feb2024/ which hopefully should help! The timeline for the full release is https://user.eumetsat.int/resources/user-guides/mtg-in-operations

[bikramb98]

@jacobbieker Using scene.available_dataset_names(), I get a list of 208 available bands for a file from the link above. Is there anyway I can find which bands out of those will be useful for OCF?
Also, I was succesful in reading a .nc file downloaded using the link above, and used the convert_scene_to_dataarray function in utils.py, and specifying a band I obtained from the list using scene.available_dataset_names(), I was able to output dataarray without any errors.

[jacobbieker]

Could you send a list of all the names? In general, we want to make a public archive of the data, so we would want all of the bands. At a minimum though, we want the 16 spectral channels that the imager takes. And that's great! Glad the code works without much changes then. Could you post the output from print(dataarray)? As I am curious what it looks like.

[bikramb98]

@jacobbieker Yes sure, here are the 208 availables bands:

mtg_scene = Scene(reader="fci_l1c_nc", filenames=[file_path])
available_bands = mtg_scene.available_dataset_names()
print(available_bands)

['ir_105', 'ir_105_earth_sun_distance', 'ir_105_index_map', 'ir_105_pixel_quality', 'ir_105_platform_altitude', 'ir_105_subsatellite_latitude', 'ir_105_subsatellite_longitude', 'ir_105_subsolar_latitude', 'ir_105_subsolar_longitude', 'ir_105_sun_satellite_distance', 'ir_105_swath_direction', 'ir_105_swath_number', 'ir_105_time', 'ir_123', 'ir_123_earth_sun_distance', 'ir_123_index_map', 'ir_123_pixel_quality', 'ir_123_platform_altitude', 'ir_123_subsatellite_latitude', 'ir_123_subsatellite_longitude', 'ir_123_subsolar_latitude', 'ir_123_subsolar_longitude', 'ir_123_sun_satellite_distance', 'ir_123_swath_direction', 'ir_123_swath_number', 'ir_123_time', 'ir_133', 'ir_133_earth_sun_distance', 'ir_133_index_map', 'ir_133_pixel_quality', 'ir_133_platform_altitude', 'ir_133_subsatellite_latitude', 'ir_133_subsatellite_longitude', 'ir_133_subsolar_latitude', 'ir_133_subsolar_longitude', 'ir_133_sun_satellite_distance', 'ir_133_swath_direction', 'ir_133_swath_number', 'ir_133_time', 'ir_38', 'ir_38_earth_sun_distance', 'ir_38_index_map', 'ir_38_pixel_quality', 'ir_38_platform_altitude', 'ir_38_subsatellite_latitude', 'ir_38_subsatellite_longitude', 'ir_38_subsolar_latitude', 'ir_38_subsolar_longitude', 'ir_38_sun_satellite_distance', 'ir_38_swath_direction', 'ir_38_swath_number', 'ir_38_time', 'ir_87', 'ir_87_earth_sun_distance', 'ir_87_index_map', 'ir_87_pixel_quality', 'ir_87_platform_altitude', 'ir_87_subsatellite_latitude', 'ir_87_subsatellite_longitude', 'ir_87_subsolar_latitude', 'ir_87_subsolar_longitude', 'ir_87_sun_satellite_distance', 'ir_87_swath_direction', 'ir_87_swath_number', 'ir_87_time', 'ir_97', 'ir_97_earth_sun_distance', 'ir_97_index_map', 'ir_97_pixel_quality', 'ir_97_platform_altitude', 'ir_97_subsatellite_latitude', 'ir_97_subsatellite_longitude', 'ir_97_subsolar_latitude', 'ir_97_subsolar_longitude', 'ir_97_sun_satellite_distance', 'ir_97_swath_direction', 'ir_97_swath_number', 'ir_97_time'

'nir_13', 'nir_13_earth_sun_distance', 'nir_13_index_map', 'nir_13_pixel_quality', 'nir_13_platform_altitude', 'nir_13_subsatellite_latitude', 'nir_13_subsatellite_longitude', 'nir_13_subsolar_latitude', 'nir_13_subsolar_longitude', 'nir_13_sun_satellite_distance', 'nir_13_swath_direction', 'nir_13_swath_number', 'nir_13_time', 'nir_16', 'nir_16_earth_sun_distance', 'nir_16_index_map', 'nir_16_pixel_quality', 'nir_16_platform_altitude', 'nir_16_subsatellite_latitude', 'nir_16_subsatellite_longitude', 'nir_16_subsolar_latitude', 'nir_16_subsolar_longitude', 'nir_16_sun_satellite_distance', 'nir_16_swath_direction', 'nir_16_swath_number', 'nir_16_time', 'nir_22', 'nir_22_earth_sun_distance', 'nir_22_index_map', 'nir_22_pixel_quality', 'nir_22_platform_altitude', 'nir_22_subsatellite_latitude', 'nir_22_subsatellite_longitude', 'nir_22_subsolar_latitude', 'nir_22_subsolar_longitude', 'nir_22_sun_satellite_distance', 'nir_22_swath_direction', 'nir_22_swath_number', 'nir_22_time'

'vis_04', 'vis_04_earth_sun_distance', 'vis_04_index_map', 'vis_04_pixel_quality', 'vis_04_platform_altitude', 'vis_04_subsatellite_latitude', 'vis_04_subsatellite_longitude', 'vis_04_subsolar_latitude', 'vis_04_subsolar_longitude', 'vis_04_sun_satellite_distance', 'vis_04_swath_direction', 'vis_04_swath_number', 'vis_04_time', 'vis_05', 'vis_05_earth_sun_distance', 'vis_05_index_map', 'vis_05_pixel_quality', 'vis_05_platform_altitude', 'vis_05_subsatellite_latitude', 'vis_05_subsatellite_longitude', 'vis_05_subsolar_latitude', 'vis_05_subsolar_longitude', 'vis_05_sun_satellite_distance', 'vis_05_swath_direction', 'vis_05_swath_number', 'vis_05_time', 'vis_06', 'vis_06_earth_sun_distance', 'vis_06_index_map', 'vis_06_pixel_quality', 'vis_06_platform_altitude', 'vis_06_subsatellite_latitude', 'vis_06_subsatellite_longitude', 'vis_06_subsolar_latitude', 'vis_06_subsolar_longitude', 'vis_06_sun_satellite_distance', 'vis_06_swath_direction', 'vis_06_swath_number', 'vis_06_time', 'vis_08', 'vis_08_earth_sun_distance', 'vis_08_index_map', 'vis_08_pixel_quality', 'vis_08_platform_altitude', 'vis_08_subsatellite_latitude', 'vis_08_subsatellite_longitude', 'vis_08_subsolar_latitude', 'vis_08_subsolar_longitude', 'vis_08_sun_satellite_distance', 'vis_08_swath_direction', 'vis_08_swath_number', 'vis_08_time', 'vis_09', 'vis_09_earth_sun_distance', 'vis_09_index_map', 'vis_09_pixel_quality', 'vis_09_platform_altitude', 'vis_09_subsatellite_latitude', 'vis_09_subsatellite_longitude', 'vis_09_subsolar_latitude', 'vis_09_subsolar_longitude', 'vis_09_sun_satellite_distance', 'vis_09_swath_direction', 'vis_09_swath_number', 'vis_09_time'

'wv_63', 'wv_63_earth_sun_distance', 'wv_63_index_map', 'wv_63_pixel_quality', 'wv_63_platform_altitude', 'wv_63_subsatellite_latitude', 'wv_63_subsatellite_longitude', 'wv_63_subsolar_latitude', 'wv_63_subsolar_longitude', 'wv_63_sun_satellite_distance', 'wv_63_swath_direction', 'wv_63_swath_number', 'wv_63_time', 'wv_73', 'wv_73_earth_sun_distance', 'wv_73_index_map', 'wv_73_pixel_quality', 'wv_73_platform_altitude', 'wv_73_subsatellite_latitude', 'wv_73_subsatellite_longitude', 'wv_73_subsolar_latitude', 'wv_73_subsolar_longitude', 'wv_73_sun_satellite_distance', 'wv_73_swath_direction', 'wv_73_swath_number', 'wv_73_time']

The output dataarray looks like this:

<xarray.DataArray (time: 1, variable: 1, y_geostationary: 853,
                   x_geostationary: 1901)>
dask.array<broadcast_to, shape=(1, 1, 853, 1901), dtype=uint16, chunksize=(1, 1, 278, 1901), chunktype=numpy.ndarray>
Coordinates:
  * y_geostationary  (y_geostationary) float64 4.229e+06 4.23e+06 ... 5.081e+06
  * x_geostationary  (x_geostationary) float64 -1.166e+06 ... 7.335e+05
  * variable         (variable) object 'vis_08'
    x_osgb           (y_geostationary, x_geostationary) float32 -6.912e+05 .....
    y_osgb           (y_geostationary, x_geostationary) float32 -4.409e+05 .....
  * time             (time) datetime64[ns] 2024-01-13T11:50:00
Attributes: (12/52)
    orbital_parameters:          {'satellite_actual_longitude': -3.5532407760...
    _FillValue:                  65535
    warm_add_offset:             0.0
    wavelength:                  0.865 µm (0.815-0.915 µm)
    sensor:                      fci
    platform_name:               MTG-I1
    ...                          ...
    vis_08_ancillary_variables:  [<xarray.DataArray 'pixel_quality' (y: 853, ...
    vis_08_start_time:           2024-01-13 11:51:37
    vis_08_end_time:             2024-01-13 11:52:30
    vis_08_reader:               fci_l1c_nc
    vis_08_area:                 Area ID: fill\nDescription: fill\nProjection...
    vis_08__satpy_id:            DataID(name='vis_08', wavelength=WavelengthR...

[bikramb98]

Do you know if .nc in the filename unique to MTG satellites? If so, I can make the Scene reader to be fci_l1c_nc in case of a .nc file

[jacobbieker]

The .nc filename is unique for EUMETSAT data for the MTG satellite so far. If we do #222 other satellites use .nc as their filetype as well, so if there is a different way to determine its the MTG file, that would be ideal.

For all the bands, I would go with just the ones that are the band name vis_04 for example, and just keep the attributes when we save the file to Zarr. As I see there is the vis_08_ancillary_variables for example that seems to hold the pixel quality mask, and other ones.

[bikramb98]

@jacobbieker I can determine if its EUMETSAT data by looking at the filename. From what I've noticed, the filename such as 'W_XX-EUMETSAT-Darmstadt,IMG+SAT,MTI1+FCI-1C-RRAD-FDHSI-FD--CHK-BODY--DIS-NC4E_C_EUMT_20240113115555_IDPFI_OPE_20240113115137_20240113115230_N_JLS_C_0072_0011.nc' contains the keyword EUMETSAT.

For the bands, I will use the ones such as vis_04, wv_63, nir_13 etc, and ignore anything with any added bits such as vis_06_sun_satellite_distance. Is that the right thing to do?

Can you explain what you mean by keeping the attributes when saving the file to Zarr? I've not worked with satellite data before so apologies for the multiple questions!

[jacobbieker]

Cool! For the bands, yes, that is what I meant!

For keeping the attributes, I mean keep these when saving out to Zarr

Attributes: (12/52)
    orbital_parameters:          {'satellite_actual_longitude': -3.5532407760...
    _FillValue:                  65535
    warm_add_offset:             0.0
    wavelength:                  0.865 µm (0.815-0.915 µm)
    sensor:                      fci
    platform_name:               MTG-I1
    ...                          ...
    vis_08_ancillary_variables:  [<xarray.DataArray 'pixel_quality' (y: 853, ...
    vis_08_start_time:           2024-01-13 11:51:37
    vis_08_end_time:             2024-01-13 11:52:30
    vis_08_reader:               fci_l1c_nc
    vis_08_area:                 Area ID: fill\nDescription: fill\nProjection...
    vis_08__satpy_id:            DataID(name='vis_08', wavelength=WavelengthR...

In general, with Xarray, attributes might disappear when doing operations. like rescaling data or combining timesteps, so we just want to make sure to keep them, as they include useful information, especially the ancillar_variables, which seems to include other datasets that give the useful information.

And no worries on the questions! Happy to answer any more, this is super helpful work!