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
import glob
import warnings
# Date format in YYYYMMDDhhmm
fnames = glob.glob('/path/to/data/H*202009060000*__')
scn = Scene(reader='seviri_l1b_hrit', filenames=fnames)
We then define a composite to load and show it.
You may see warnings due to deprecated pyproj features and when numpy/dask process NaNs values. NaNs are used to mark invalid pixels and space pixels (like you'll see in a full disk image).
You can silence these with warnings.filterwarnings('ignore')
.
composite = 'natural_color'
scn.load([composite])
scn.show(composite)
<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:
scn.available_composite_ids()
[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)]
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)
<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.
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')
/home/a001673/usr/src/dask/dask/local.py:271: RuntimeWarning: invalid value encountered in divide return func(*args2)
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)
<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
local_scn.show('ndvi')
local_scn.show('natural')