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 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.
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='seviri_l1b_hrit')
scn = Scene(filenames=files)
We then define a composite to load and show it:
composite = 'natural'
scn.load([composite])
scn.show(composite)
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:
scn.available_composite_ids()
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:
scn.load(['VIS006', 0.8])
print(scn)
Working with the data is quite straightforward as the datasets inside the scene are in effect DataArrays as per the xarray python package.
from satpy.dataset import combine_metadata
ndvi = (scn[0.8] - scn[0.6]) / (scn[0.8] + scn[0.6])
ndvi.attrs = combine_metadata(scn[0.8], scn[0.6])
scn['ndvi'] = ndvi
scn.show('ndvi')
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:
local_scn = scn.resample("eurol")
print(local_scn)
local_scn.show('ndvi')
local_scn.show('natural')