Quickstart with MSG SEVIRI

For this tutorial, we will use Meteosat Second Generation (MSG) data in uncompressed EUMETSAT HRIT format, read it with satpy, resample it with pyresample and process it a bit. Be sure to have the necessary python packages installed, using eg:

pip install satpy

Software to uncompress HRIT data is callled Public Wavelet Transform Decompression Library Software and can be obtained from EUMETSAT on their software page. Compressed HRIT data is easily recognizable as the files end with C_, while uncompressed data files end with __ (two underscores).

Loading data and generating our first image

Loading the data is done on-the-fly with satpy when a composite is to be generated. So first we create a scene instance, pointing the base_dir directory to the place containing the uncompressed HRIT data files.

In [4]:
from satpy.scene import Scene
from satpy.resample import get_area_def
from satpy import find_files_and_readers
from datetime import datetime

files = find_files_and_readers(base_dir='../hrit_seviri/test_data', 
                               start_time=datetime(2017, 1, 19, 9, 30),
                               end_time=datetime(2017, 1, 19, 9, 45),
                               reader='hrit_msg')
scn = Scene(filenames=files)

We then define a composite to load and show it:

In [5]:
composite = 'natural'
scn.load([composite])
scn.show(composite)
Out[5]:
<trollimage.xrimage.XRImage at 0x4a44c90>

As you can see, earth is displayed with North facing down: remember that this is raw data from the hrit files, unprojected, so the pixels are displayed in scanning order.

To get a list of available composites to choose from:

In [4]:
scn.available_composite_ids()
Out[4]:
[DatasetID(name='airmass', wavelength=None, resolution=None, polarization=None, calibration=None, modifiers=None),
 DatasetID(name='airmass_corr', wavelength=None, resolution=None, polarization=None, calibration=None, modifiers=None),
 DatasetID(name='ash', wavelength=None, resolution=None, polarization=None, calibration=None, modifiers=None),
 DatasetID(name='cloudtop', wavelength=None, resolution=None, polarization=None, calibration=None, modifiers=None),
 DatasetID(name='convection', wavelength=None, resolution=None, polarization=None, calibration=None, modifiers=None),
 DatasetID(name='day_microphysics', wavelength=None, resolution=None, polarization=None, calibration=None, modifiers=None),
 DatasetID(name='dust', wavelength=None, resolution=None, polarization=None, calibration=None, modifiers=None),
 DatasetID(name='fog', wavelength=None, resolution=None, polarization=None, calibration=None, modifiers=None),
 DatasetID(name='green_snow', wavelength=None, resolution=None, polarization=None, calibration=None, modifiers=None),
 DatasetID(name='ir_overview', wavelength=None, resolution=None, polarization=None, calibration=None, modifiers=None),
 DatasetID(name='natural', wavelength=None, resolution=None, polarization=None, calibration=None, modifiers=None),
 DatasetID(name='natural_sun', wavelength=None, resolution=None, polarization=None, calibration=None, modifiers=None),
 DatasetID(name='night_fog', wavelength=None, resolution=None, polarization=None, calibration=None, modifiers=None),
 DatasetID(name='night_microphysics', wavelength=None, resolution=None, polarization=None, calibration=None, modifiers=None),
 DatasetID(name='overview', wavelength=None, resolution=None, polarization=None, calibration=None, modifiers=None),
 DatasetID(name='overview_sun', wavelength=None, resolution=None, polarization=None, calibration=None, modifiers=None),
 DatasetID(name='realistic_colors', wavelength=None, resolution=None, polarization=None, calibration=None, modifiers=None)]

Loading channel data and working with it

In order to load channel data to work with, the procedure is identical to what we presented above, except the name of a channel or its wavelength (along with the full-featured DatasetID) has to be passed. For example:

In [5]:
scn.load(['VIS006', 0.8])
print(scn)
<xarray.DataArray 'astype-7a343ed713478b648e6dadb526d4c114' (y: 3712, x: 3712)>
dask.array<concatenate, shape=(3712, 3712), dtype=float64, chunksize=(53, 3712)>
Coordinates:
  * x        (x) float64 -5.569e+06 -5.569e+06 -5.569e+06 -5.569e+06 ...
  * y        (y) float64 5.566e+06 5.566e+06 5.566e+06 5.566e+06 5.566e+06 ...
Attributes:
    units:                %
    wavelength:           (0.74, 0.81, 0.88)
    standard_name:        toa_bidirectional_reflectance
    platform_name:        Meteosat-10
    sensor:               seviri
    satellite_longitude:  0.0
    satellite_latitude:   0.0
    satellite_altitude:   35785831.0
    start_time:           2017-01-19 09:30:10.553000
    end_time:             2017-01-19 09:42:41.403000
    area:                 Area ID: some_area_name\nDescription: On-the-fly ar...
    name:                 VIS008
    resolution:           3000.40316582
    calibration:          reflectance
    polarization:         None
    modifiers:            ()
    ancillary_variables:  []
<xarray.DataArray 'astype-5c596270ef4d3597373e262146cb262e' (y: 3712, x: 3712)>
dask.array<concatenate, shape=(3712, 3712), dtype=float64, chunksize=(53, 3712)>
Coordinates:
  * x        (x) float64 -5.569e+06 -5.569e+06 -5.569e+06 -5.569e+06 ...
  * y        (y) float64 5.566e+06 5.566e+06 5.566e+06 5.566e+06 5.566e+06 ...
Attributes:
    units:                %
    wavelength:           (0.56, 0.635, 0.71)
    standard_name:        toa_bidirectional_reflectance
    platform_name:        Meteosat-10
    sensor:               seviri
    satellite_longitude:  0.0
    satellite_latitude:   0.0
    satellite_altitude:   35785831.0
    start_time:           2017-01-19 09:30:10.553000
    end_time:             2017-01-19 09:42:41.403000
    area:                 Area ID: some_area_name\nDescription: On-the-fly ar...
    name:                 VIS006
    resolution:           3000.40316582
    calibration:          reflectance
    polarization:         None
    modifiers:            ()
    ancillary_variables:  []
<xarray.DataArray 'natural' (bands: 3, y: 3712, x: 3712)>
dask.array<concatenate, shape=(3, 3712, 3712), dtype=float64, chunksize=(1, 53, 3712)>
Coordinates:
  * x        (x) float64 -5.569e+06 -5.569e+06 -5.569e+06 -5.569e+06 ...
  * y        (y) float64 5.566e+06 5.566e+06 5.566e+06 5.566e+06 5.566e+06 ...
  * bands    (bands) |S1 'R' 'G' 'B'
Attributes:
    units:                   %
    platform_name:           Meteosat-10
    sensor:                  seviri
    satellite_longitude:     0.0
    satellite_latitude:      0.0
    satellite_altitude:      35785831.0
    start_time:              2017-01-19 09:30:00
    end_time:                2017-01-19 09:45:00
    area:                    Area ID: some_area_name\nDescription: On-the-fly...
    polarization:            None
    ancillary_variables:     []
    metadata_requirements:   []
    standard_name:           natural
    wavelength:              None
    optional_prerequisites:  []
    mode:                    RGB
    name:                    natural
    prerequisites:           [1.63, 0.85, 0.635]
    optional_datasets:       []
    resolution:              None
    calibration:             None
    modifiers:               None

Working with the data is quite straightforward as the datasets inside the scene are in effect DataArrays as per the xarray python package.

In [6]:
from satpy.dataset import combine_attrs
ndvi = (scn[0.8] - scn[0.6]) / (scn[0.8] + scn[0.6])
ndvi.attrs = combine_attrs(scn[0.8], scn[0.6])
scn['ndvi'] = ndvi
scn.show('ndvi')
/home/a001673/usr/src/dask/dask/local.py:271: RuntimeWarning: invalid value encountered in divide
  return func(*args2)

Resampling the data

Until now, we have used the channels directly as provided by the satellite, that is in satellite projection. Generating composites thus produces views in satellite projection, i.e. as viewed by the satellite.

Most often however, we will want to project the data onto a specific area so that only the area of interest is depicted in the RGB composites.

Here is how we do that:

In [7]:
local_scn = scn.resample("eurol")
print(local_scn)
<xarray.DataArray 'array-6cc44ebc64d08c71c8e71a0c2b12bf7a' (y: 2048, x: 2560)>
dask.array<transpose, shape=(2048, 2560), dtype=float64, chunksize=(1000, 1000)>
Coordinates:
  * y        (y) float64 -1.502e+06 -1.504e+06 -1.508e+06 -1.510e+06 ...
  * x        (x) float64 -3.778e+06 -3.776e+06 -3.772e+06 -3.77e+06 ...
Attributes:
    units:                %
    wavelength:           (0.74, 0.81, 0.88)
    standard_name:        toa_bidirectional_reflectance
    platform_name:        Meteosat-10
    sensor:               seviri
    satellite_longitude:  0.0
    satellite_latitude:   0.0
    satellite_altitude:   35785831.0
    start_time:           2017-01-19 09:30:10.553000
    end_time:             2017-01-19 09:42:41.403000
    area:                 Area ID: eurol\nDescription: Euro 3.0km area - Euro...
    name:                 VIS008
    resolution:           3000.40316582
    calibration:          reflectance
    polarization:         None
    modifiers:            ()
    ancillary_variables:  []
<xarray.DataArray 'array-c0f9f5bf6d34cf1e81f951008b71ce17' (y: 2048, x: 2560)>
dask.array<transpose, shape=(2048, 2560), dtype=float64, chunksize=(1000, 1000)>
Coordinates:
  * y        (y) float64 -1.502e+06 -1.504e+06 -1.508e+06 -1.510e+06 ...
  * x        (x) float64 -3.778e+06 -3.776e+06 -3.772e+06 -3.77e+06 ...
Attributes:
    units:                %
    wavelength:           (0.56, 0.635, 0.71)
    standard_name:        toa_bidirectional_reflectance
    platform_name:        Meteosat-10
    sensor:               seviri
    satellite_longitude:  0.0
    satellite_latitude:   0.0
    satellite_altitude:   35785831.0
    start_time:           2017-01-19 09:30:10.553000
    end_time:             2017-01-19 09:42:41.403000
    area:                 Area ID: eurol\nDescription: Euro 3.0km area - Euro...
    name:                 VIS006
    resolution:           3000.40316582
    calibration:          reflectance
    polarization:         None
    modifiers:            ()
    ancillary_variables:  []
<xarray.DataArray 'array-504d92bc846a1c334e0d3f610af353fc' (y: 2048, x: 2560)>
dask.array<transpose, shape=(2048, 2560), dtype=float64, chunksize=(1000, 1000)>
Coordinates:
  * y        (y) float64 -1.502e+06 -1.504e+06 -1.508e+06 -1.510e+06 ...
  * x        (x) float64 -3.778e+06 -3.776e+06 -3.772e+06 -3.77e+06 ...
Attributes:
    satellite_latitude:   0.0
    polarization:         None
    modifiers:            ()
    satellite_longitude:  0.0
    area:                 Area ID: eurol\nDescription: Euro 3.0km area - Euro...
    sensor:               seviri
    satellite_altitude:   35785831.0
    calibration:          reflectance
    ancillary_variables:  []
    platform_name:        Meteosat-10
    standard_name:        toa_bidirectional_reflectance
    end_time:             2017-01-19 09:42:41.403000
    units:                %
    resolution:           3000.40316582
    start_time:           2017-01-19 09:30:10.553000
    name:                 ndvi
<xarray.DataArray 'array-62f880414388c2e236e67307bea6deb6' (bands: 3, y: 2048, x: 2560)>
dask.array<transpose, shape=(3, 2048, 2560), dtype=float64, chunksize=(3, 1000, 1000)>
Coordinates:
  * bands    (bands) |S1 'R' 'G' 'B'
  * y        (y) float64 -1.502e+06 -1.504e+06 -1.508e+06 -1.510e+06 ...
  * x        (x) float64 -3.778e+06 -3.776e+06 -3.772e+06 -3.77e+06 ...
Attributes:
    units:                   %
    platform_name:           Meteosat-10
    sensor:                  seviri
    satellite_longitude:     0.0
    satellite_latitude:      0.0
    satellite_altitude:      35785831.0
    start_time:              2017-01-19 09:30:00
    end_time:                2017-01-19 09:45:00
    area:                    Area ID: eurol\nDescription: Euro 3.0km area - E...
    polarization:            None
    ancillary_variables:     []
    metadata_requirements:   []
    standard_name:           natural
    wavelength:              None
    optional_prerequisites:  []
    mode:                    RGB
    name:                    natural
    prerequisites:           [1.63, 0.85, 0.635]
    optional_datasets:       []
    resolution:              None
    calibration:             None
    modifiers:               None
In [8]:
local_scn.show('ndvi')

In [10]:
local_scn.show('natural')