Remote Sensing I

SatPy

A Python Library for Weather Satellite Processing

From the documentation:

SatPy is a python library for reading and manipulating meteorological remote sensing data and writing it to various image and data file formats.

SatPy comes with the ability to make various RGB composites directly from satellite instrument channel data or higher level processing output.

The pyresample package is used to resample data to different uniform areas or grids.

Various atmospheric corrections and visual enhancements are also provided, either directly in Satpy or from those in the PySpectral and TrollImage packages.

Load & view MSG data

SatPy is well suited for working with data of Meteosat Second Generation Satellites:

First, we specify the file(s) of the scene we want to open together with the reader of the respective file type (in our case we use the netcdf reader for the SEVIRI instrument):

Then we instantiate a Scene object with the given files:

The scene attributes tell us what instrument we are dealing with and the start- and end-times of the scene scan:

To get an overview of all available datasets (=channels/bands) of the scene, we can call the all_dataset_names() function:

Before we are able to work with the datasets, we first have to load the bands we are actually interested in. Let’s say, we are only interested in the visible wavelength at 0.6µm:

We can have a look at all loaded bands by printing out the keys:

Let’s have a first look at what we are dealing with:

Ok, thats nice :)

However, SatPy’s show() function doesn’t offer the option to depict a colorbar (so that we get an impression of the value range of the depicted colors).

Fortunately, we can easily achieve this by using the imshow() or pcolormesh() functions from xarray (which are just thin wrappers around the same-named matplotlib functions). For most of our use-cases it’s better to use imshow() as it is much faster when we are dealing with large data sets.

Well, this is rahter ugly and the colorbar title is misleading (although that’s not xarray’s fault).

Like always, if we want to control specific plot elements, we can just use matplotlib itself. We need more lines of code with this approach but the result looks better:

It seems as if we are dealing with reflectance values here (range goes from 0% to 100%). Let’s go sure…

To get an overview of the metadata of a loaded band, you can just use the print() function on the band:

This shows us some valuable information about the data:

• The data is stored within an xarray.DataArray (which itself uses a dask.array to store the sensor data)
• The band has 2 dimensions x and y (it is an image)
• Over both dimensions, the image spreads over 3712 pixels
• Pixel values are stored as float32
• Values are already converted to “reflectance” between 0% and 100%. So SatPy already converted the original raw counts to reflectance values.
• And a lot more (like satellite and sensor name, resolution, …)

Calibration

By default Satpy will provide the version of the dataset with the highest resolution and the highest level of calibration (brightness temperature or reflectance over radiance). It is also possible to request one or more of these exact versions of a dataset by specifying the calibration argument of the load() function:

Coordinates

Let’s go back to the visible band at 0.6μm: From the metadata above, we can also derive that there is coordinate information provided with the data:

lat and lon are coordinate variables, each consisting of a full 3712x3712 layer, spanning along the dimensions x and y.

We can have a closer look at the coordinate layers by printing them to the console:

So here we see, that lat is given in units of “degrees_north” and that valid values are between -90° and +90°.

For the sake of completeness we can also plot the coordinate layers:

Composites

Built-in composites

SatPy comes with many built-in utility functions to create composite images of the data.

To get an overview of all available composites, you can do:

As with the single bands, you have to load a composite into the Scene object first, if you want to work with it:

Let’s have a look at the metadata of the composite:

• We see, that the composite consists of 3 bands, stored in a dask-array with the shape (3, 3712, 3712)
• Each band has 2 dimensions x and y
• Over both dimensions, the images spread over 3712 pixels
• Pixel values are stored as float32
• We also see, that the 3 bands correspond to the red, green and blue layers of the image that we plottetd above.
• These bands are: IR_016, VIS008 and VIS006 (prerequisites)

Guides on how to interpret these composites can be found on the Eumetsat website:

https://www.eumetsat.int/website/home/Data/Training/TrainingLibrary/DAT_2044069.html

Let’s load all available bands (except HRV) and the natural color composite so we can work with these data sets later on:

User-defined composites

Additionally to the many built-in composites, you can define your own composites using the Compositer classes of SatPy.

Let’s create an infrared overview using the bands at 3.9μm, 10.8μm and 12.0μm:

That’s it.

However, a composite itself is just defined as a special dataset which may have several bands (like R, G and B bands). The data isn’t stretched, or clipped or gamma filtered until an image is generated. To get an image out of the above composite, you have to convert it to an image object and tweak it’s display properties:

Cropping

Often, you are only interested in a specific region of the fulldisk.

Let’s say we want to take a closer look at Greece. We can crop the fulldisk to a specific region using the crop() function of SatPy’s Scene class:

The cropped scene keeps all original metadata. Only the bands, composites and coordinates/dimensions are cropped to the specified extents:

Resampling

Cropping is nice and all, but due to the satellite viewing angle, the cropped scene can look very distorted, especially when we are working with the margin regions of the scene. To get a better looking view of the domain of interest, you can use the resampling utilities of SatPy (or pyresample).

To resample the data to a specific area and projection, you can specify one of the many predefined areas of SatPy. All predefined areas can be found in SatPy’s satpy/etc/areas.yaml file.

For example, let’s take libya:

This resamples the scene and all loaded bands/composites to the extents of the “libya” area. If you look into the areas.yaml file, you see how it is defined:

libya:
  description: libya area
  projection:
    proj: merc
    lat_ts: 31.625
    lon_0: 17.875
    ellps: WGS84
  shape:
    height: 1024
    width: 1280
  area_extent:
    lower_left_xy: [-1921632.0902750609, 1725320.2028891125]
    upper_right_xy: [1918367.9097249391, 4797320.202889113]
    units: m

Let’s plot the natural color composite again, to get an overview of where we are:

Nice, but let’s go back to Greece!

Unfortunately, there is no predefined area definition in the yaml-File for Greece.

Luckily, we can also define our own area of interest! So, let’s resample the data using the Mercator projection:

For more information on resampling, visit the pyresample and PROJ.4 documentations.

Writing to disk

SatPy makes it possible to save datasets in multiple formats. Most of the time, the library is used to save the bands or composites of a satellite scene as images (e.g. PNG or TIFF) which are then provided to end users. You can do this using one of the built-in SatPy writers.

For example, let’s write the resampled band at 0.6μm to a simple greyscale PNG image:

You can also save multiple or all datasets of a scene simultaneously using the save_datasets() method:

You can also easily save your data in formats like NetCDF or GeoTIFF for further manipulation in other software tools or for distribution to third parties that are not used to work with satellite data types. So let’s save the “VIS006”, “IR_108” and the “natural_color” bands of the original scene in a netcdf-file:

And in the same way you can also save the data in geotiff files:

Exercise 6

• Complete the sixth assignment and push your results until tuesday 14:00 next week
• Assignment link: https://classroom.github.com/a/1-NlM77w