Remote Sensing II

Advanced plotting methods

In the last session, we only plotted the satellite data of MSG itself. Unfortunately, we have to follow our instincts (and duty) as Geographers and add some spatial context to the plots so that the viewer knows what exactly she or he sees.

So lets start with the natural color composite of Greece.

First, we load MSG data again and resample the scene to our area definition of Greece and its surroundings:

from satpy import Scene, MultiScene
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from pyresample import geometry

dateien = ["../data/W_XX-EUMETSAT-Darmstadt,VIS+IR+HRV+IMAGERY,MSG1+SEVIRI_C_EUMG_20060117084508.nc"]
files = {"seviri_l1b_nc" : dateien}

scn = Scene(filenames=files)
scn.load(["natural_color"])
scn.load(scn.all_dataset_names()[1:])

area_id = "Griechenland"
description = "Griechenland und Umgebung in Mercator-Projektion"
proj_id = "Griechenland"
proj_dict = {"proj": "merc", "lat_ts": 38, 'lon_0': 25}

width = 800    # width of the result domain in pixels
height = 800   # height of the result domain in pixels

llx = -10E5   # projection x coordinate of lower left corner of lower left pixel
lly =  27E5   # projection y coordinate of lower left corner of lower left pixel
urx =  10E5   # projection x coordinate of upper right corner of upper right pixel
ury =  47E5   # projection y coordinate of upper right corner of upper right pixel

area_extent = (llx,lly,urx,ury)
area_def_greece = geometry.AreaDefinition(area_id, proj_id, description, proj_dict, width, height, area_extent)

local_scn = scn.resample(area_def_greece)

Now, we can plot the scene and add some coastlines and a lat/lon grid:

# First, we have to transpose the natural color composite values to a shape that can be interpreted 
# by the imshow method: (M,N,3)
image = np.asarray(local_scn["natural_color"]).transpose(1,2,0)

# Then we scale the values to the range between 0 and 1, clipping the lower and upper percentiles
# so that a potential contrast decrease caused by outliers is eliminated.
image = np.interp(image, (np.percentile(image,1), np.percentile(image,99)), (0, 1))

# Now we "copy" the coordinate reference system of our composite data set...
crs = local_scn["natural_color"].attrs["area"].to_cartopy_crs()

# ... and use it to generate an axes in our figure with the same CRS
fig = plt.subplots(figsize=(10,10))
ax = plt.axes(projection=crs)

# Now we can add some coastlines...
ax.coastlines(resolution="10m", color="white")

# ... and a lat/lon grid:
ax.gridlines(xlocs=range(10,45,5),ylocs=range(25,55,5))

# In the end, we can plot our image data...
ax.imshow(image, transform=crs, extent=crs.bounds, origin="upper")

# and add a title to our plot
plt.title("Natural color composite of Greece and surroundings, recorded by MSG at " + local_scn.attrs["start_time"].strftime("%Y-%m-%d %H:%M"))

# Finally, we can show the plot to the user:
plt.show()

png

If you want to add further shapes to the plot, you can do so by using the feature class of cartopy:

import cartopy.feature as cfeature

image = np.asarray(local_scn["natural_color"]).transpose(1,2,0)
image = np.interp(image, (np.percentile(image,1), np.percentile(image,99)), (0, 1))
crs = local_scn["natural_color"].attrs["area"].to_cartopy_crs()
fig = plt.subplots(figsize=(10,10))
ax = plt.axes(projection=crs)
ax.gridlines(xlocs=range(10,45,5),ylocs=range(25,55,5))

# ===================================================
states_provinces = cfeature.NaturalEarthFeature(
        category="cultural",
        name="admin_0_countries",
        scale="10m",
        facecolor="none")
ax.add_feature(states_provinces, edgecolor="white")
# ===================================================

ax.imshow(image, transform=crs, extent=crs.bounds, origin="upper")
plt.title("Natural color composite of Greece and surroundings, recorded by MSG at " + local_scn.attrs["start_time"].strftime("%Y-%m-%d %H:%M"))
plt.show()

png

Task
  1. Download this MSG data set using our usual password: https://box.uni-marburg.de/index.php/s/bCu96BUfVyXtFCg
  2. Resample the scene to Algeria and its surroundings using the Mercator projection. Use the following settings:
    • Latitude of true scale: 28
    • Longitude of projection center: 2
    • Domain size: 1000 x 1000
    • Area extent: [-13E5, 18E5, 12E5, 41E5]
  3. Plot the resampled scene as natural color composite, add a grid and country borders as well as a title.

Solution

  1. area_def_algeria = geometry.AreaDefinition("Algeria", "", "", 
                            {"proj": "merc", "lat_ts": 28, "lon_0": 2}, 
                            1000, 1000, 
                            (-13E5,18E5,12E5,41E5))
local_scn = scn.resample(area_def_algeria)

  2. image = np.asarray(local_scn["natural_color"]).transpose(1,2,0)
image = np.interp(image, (np.percentile(image,1), np.percentile(image,99)), (0, 1))
crs = local_scn["natural_color"].attrs["area"].to_cartopy_crs()
fig = plt.subplots(figsize=(10,10))
ax = plt.axes(projection=crs)
ax.gridlines(xlocs=range(-15,30,5),ylocs=range(15,47,5))
states_provinces = cfeature.NaturalEarthFeature(
  category="cultural",
  name="admin_0_countries",
  scale="50m",
  facecolor="none")
ax.add_feature(states_provinces, edgecolor="white")
ax.imshow(image, transform=crs, extent=crs.bounds, origin="upper")
plt.title("Natural color composite of Algeria and surroundings, recorded by MSG at " + local_scn.attrs["start_time"].strftime("%Y-%m-%d %H:%M"))
plt.show()

Working with multiple scenes (time series)

Often, we are interested not only in one specific scene but in the changes over consecutively recorded scenes.

For this purpose, we can make use of SatPy’s (experimental) MultiScene:

from satpy.multiscene import timeseries
from glob import glob

# Load all files in a given folder:
files = glob("../data/W_XX-EUMETSAT-Darmstadt,VIS+IR+HRV+IMAGERY,MSG1+SEVIRI_C_EUMG_*.nc")

# Create a MultiScene from the given list of scenes
scenes = [Scene(reader="seviri_l1b_nc", filenames=[f]) for f in files]
mscn = MultiScene(scenes)

# Load the 10.8 micrometer band of all scenes within the MultiScene
mscn.load(["IR_108"])

# Resample the MultiScene to our Greece domain
new_mscn = mscn.resample(area_def_greece)

# Blend the scenes to one single Scene with each dataset/channel extended by the time dimension
blended_scene = new_mscn.blend(blend_function=timeseries)

Now we can, for example, plot the loaded band of each scene:

blended_scene["IR_108"].plot(x="x", y="y", col="time", col_wrap=3)

<xarray.plot.facetgrid.FacetGrid at 0x7eff46a5cb00>

png

Loading MODIS data

MODIS data can be obtained at the LAADS DAAC (registration required): https://ladsweb.modaps.eosdis.nasa.gov/search/

MOD021KM

The MODIS MOD021KM data (Level 1B Calibrated Radiances at 1km resolution) can be read and displayed just like the MSG data using SatPy:

dateien = ["../data/MODIS/MOD021KM.A2006017.0855.006.2014218155904.hdf"]
files = {"modis_l1b" : dateien}

mod02_scn = Scene(filenames=files)
mod02_scn.load(["natural_color"],resolution=1000)
mod02_scn.show("natural_color")

This gives us an image of the swath data in original MODIS projection.

But where are we? To find that out, we can resample the data to a simple world map:

local_scn = mod02_scn.resample("worldeqc30km")

image = np.asarray(local_scn["natural_color"]).transpose(1,2,0)
image = np.interp(image, (np.nanpercentile(image,1), np.nanpercentile(image,99)), (0, 1))

crs = local_scn["natural_color"].attrs["area"].to_cartopy_crs()

fig = plt.subplots(figsize=(15,8))
ax = plt.axes(projection=crs)
ax.coastlines(resolution="10m",color="grey")
ax.imshow(image, transform=crs, extent=crs.bounds, origin="upper")
plt.show()

png

As MODIS has a very wide aperture angle, the resolution towards the scan margins drops considerably and we can see a lot of distortion at the left and right image borders of the original swath data.

Most often, however, you will want to reproject MODIS data to another projection that is suitable for the covered domain.

This can be done exactly like with the MSG data. So lets resample our MODIS scene to the Greece domain that we defined above:

local_scn = mod02_scn.resample(area_def_greece)
local_scn.show("natural_color")

png

Stack swath data

Well ok, that only covers parts of the domain. To better fill the domain, we have to load multiple swath sections.

For this, again, we can use the MultiScene object:

scenes = [
    Scene(reader='modis_l1b', filenames=["../data/MODIS/MOD021KM.A2006017.0855.006.2014218155904.hdf"]),
    Scene(reader='modis_l1b', filenames=["../data/MODIS/MOD021KM.A2006017.0850.006.2014218155843.hdf"])
    ]

mscn = MultiScene(scenes)
mscn.load(['natural_color'],resolution=1000)
mscn_res = mscn.resample(area_def_greece)
blended_scene = mscn_res.blend()
blended_scene.show("natural_color")

png

And again, with some tweaking, we can also add additional information to the image using CartoPy:

image = np.asarray(blended_scene["natural_color"]).transpose(1,2,0)
image = np.interp(image, (np.nanpercentile(image,1), np.nanpercentile(image,99)), (0, 1))
crs = blended_scene["natural_color"].attrs["area"].to_cartopy_crs()
fig = plt.subplots(figsize=(10,10))
ax = plt.axes(projection=crs)
ax.gridlines(xlocs=range(10,45,5),ylocs=range(25,55,5))
states_provinces = cfeature.NaturalEarthFeature(
    category="cultural",
    name="admin_0_countries",
    scale="50m",
    facecolor="none")
ax.add_feature(states_provinces, edgecolor="white")
ax.imshow(image, transform=crs, extent=crs.bounds, origin="upper")
plt.title("Natural color composite of Greece and surroundings, recorded by MODIS at " + local_scn.attrs["start_time"].strftime("%Y-%m-%d %H:%M"))
plt.show()

png

Task
  1. Download and read the MOD021KM data at: https://box.uni-marburg.de/index.php/s/9079aAAbGRYWurA
  2. Plot the natural color composite in the Algeria domain that you defined before.

Solution

  1. mod02_scn = Scene({"modis_l1b" : ["path/to/MODIS-data.hdf"]})
mod02_scn.load(["natural_color"],resolution=1000)
local_scn = mod02_scn.resample(area_def_algeria)
local_scn.show("natural_color")

MOD06_L2

Currently, MOD06_L2 data cannot be read by SatPy. Fortunately, we can circumvent this problem and use xarray directly.

So, let’s open an example data set:

import xarray as xr
modis_l2 = xr.open_dataset("../data/MODIS/MOD06_L2.A2006017.0855.061.2017271194536.hdf")

That’s it. There are many data arrays within this data set that you can list with the data_vars attribute. But let’s only plot the cloud fraction:

modis_l2["Cloud_Fraction"].plot.imshow(origin="upper")
Task
  1. Download and read the MOD06_L2 data at: https://box.uni-marburg.de/index.php/s/MXrVcvOj9AoFESz
  2. Find out the start and end times of the scan.
  3. Plot the cloud optical thickness data set
  1. modis_l2 = xr.open_dataset("../path/to/dataset.hdf")

  2. modis_l2.Scan_Start_Time[0,0].values
modis_l2.Scan_Start_Time[-1,-1].values

  3. modis_l2["Cloud_Optical_Thickness"].plot.imshow(origin="upper")

Ok, looks like it worked. To resample the data, we now have to use the lat/lon bands that come with the MOD06_L2 data set, create a Swath object and use that as area definition in a self-made SatPy Scene object:

from datetime import datetime

# we need to rename dimensions to x,y because satpy expects names to be that way
new_dim_names = {"Cell_Along_Swath_5km:mod06": "y", "Cell_Across_Swath_5km:mod06": "x"}
modis_l2 = modis_l2.rename(new_dim_names)

# For now, we only want to take a look at the cloud fraction
data = modis_l2["Cloud_Fraction"]

# optionally similar to information found in the mod021km we can assign some more informational attributes
data.attrs["coordinates"] = ("longitude", "latitude")
data.attrs["resolution"] = 5000

# set the start time of the scan (we need it for plotting later)
# since xarray create a np.datetime datatype we convert it back to a datetime.datetime object which
# Satpy understands (important for saving).
st = modis_l2.Scan_Start_Time[0,0].values
data.attrs["start_time"] = datetime.utcfromtimestamp((st - np.datetime64("1970-01-01T00:00:00")) / np.timedelta64(1, "s"))

# create a swath area from lon/lats
lons = modis_l2["Longitude"]
lats = modis_l2["Latitude"]
swath = geometry.SwathDefinition(lons,lats)

# and assign the SwathDefinition to the area attribute of the data array
data.attrs["area"] = swath

# init an empty satpy Scene
mod06_scn = Scene()

# assign a new Scene dataset
mod06_scn["cloud_fraction"] = data

Now that we have generated a SatPy Scene object, we can use all its functionality, e.g. resampling:

mod06_scn_res = mod06_scn.resample(area_def_greece, resampler="nearest", radius_of_influence=10000)

mod06_scn_res.show("cloud_fraction")

And then we can plot it as usual:

crs = mod06_scn_res["cloud_fraction"].attrs["area"].to_cartopy_crs()

fig = plt.subplots(figsize=(10,10))
ax = plt.axes(projection=crs)
#ax.gridlines(xlocs=range(-15,30,5),ylocs=range(15,47,5))
states_provinces = cfeature.NaturalEarthFeature(
    category="cultural",
    name="admin_0_countries",
    scale="50m",
    facecolor="none")
ax.add_feature(states_provinces, edgecolor="black")
ax.imshow(mod06_scn_res["cloud_fraction"], transform=crs, extent=crs.bounds, origin="upper")
plt.title("Cloud fraction over Greece and surroundings, recorded by MODIS at " + np.datetime_as_string(mod06_scn_res.attrs["start_time"], timezone='UTC'))
plt.show()

png

Band calculations

As SatPy’s Scene class internally uses the xarray data structure to store its data sets, you can do various calculations on these bands with ease.

Often, you will need the difference valus of two bands, which help to highlight certain atmospheric features. For example, the difference values between the bands at 10.8μm and 3.9μm help to distinguish low stratus layers from higher clouds at night.

So, let’s first load the MSG scene at 00:00 UTC on 17. January 2006 and resample it to the area of Greece again:

dateien = ["../data/W_XX-EUMETSAT-Darmstadt,VIS+IR+HRV+IMAGERY,MSG1+SEVIRI_C_EUMG_20060117000010.nc"]
files = {"seviri_l1b_nc" : dateien}

scn = Scene(filenames=files)
scn.load(scn.all_dataset_names()[1:])

local_scn = scn.resample(area_def_greece)

We can now calculate the difference values by simply subtracting the bands from each other with the - operator and we can store them in a new data set within the Scene object:

local_scn["band_diff"] = local_scn["IR_108"]-local_scn["IR_039"]

Let’s see what we got:

local_scn["band_diff"].plot.imshow()

<matplotlib.image.AxesImage at 0x7fe77afe56a0>

png

High difference values indicate low stratus clouds. This means, the dark-red regions in the plot above are covered by low stratus.

By applying a threshold value (~5K), we can generate a low stratus cloud mask for our current scene:

local_scn["stratus_mask"] = local_scn["band_diff"] > 5
local_scn["stratus_mask"].plot.imshow(cmap="Greys_r")

<matplotlib.image.AxesImage at 0x7fe77af30630>

png

Exercise 7