Satellite Data Visualization

Satellite Data Visualization#

Tutorial Leads: Carina Poulin (NASA, SSAI), Ian Carroll (NASA, UMBC) and Sean Foley (NASA, MSU)

The following notebooks are prerequisites for this tutorial:

Earthdata Cloud Access

Summary#

This is an introductory tutorial to the visualization possibilities arising from PACE data, meant to give you ideas and tools to develop your own scientific data visualizations.

Learning Objectives#

At the end of this notebook you will know:

How to create an easy global map from OCI data from the cloud
How to create a true-color image from OCI data processed with OCSSW
How to make a false color image to look at clouds
How to make an interactive tool to explore OCI data
What goes into an animation of multi-angular HARP2 data

1. Setup#

Begin by importing all of the packages used in this notebook. If your kernel uses an environment defined following the guidance on the tutorials page, then the imports will be successful.

import os

from holoviews.streams import Tap
from matplotlib import animation
from matplotlib.colors import ListedColormap
from PIL import Image, ImageEnhance
from scipy.ndimage import gaussian_filter1d
from xarray.backends.api import open_datatree
import cartopy.crs as ccrs
import cmocean
import earthaccess
import holoviews as hv
import matplotlib.pyplot as plt
import matplotlib.pylab as pl
import numpy as np
import panel.widgets as pnw
import xarray as xr

options = xr.set_options(display_expand_attrs=False)

In this tutorial, we suppress runtime warnings that show up when calculating log for negative values, which is common with our datasets.

Define a function to apply enhancements, our own plus generic image enhancements from the Pillow package.

def enhance(rgb, scale = 0.01, vmin = 0.01, vmax = 1.04, gamma=0.95, contrast=1.2, brightness=1.02, sharpness=2, saturation=1.1):
    """The SeaDAS recipe for RGB images from Ocean Color missions.

    Args:
        rgb: a data array with three dimensions, having 3 or 4 bands in the third dimension
        scale: scale value for the log transform
        vmin: minimum pixel value for the image
        vmax: maximum pixel value for the image
        gamma: exponential factor for gamma correction
        contrast: amount of pixel value differentiation 
        brightness: pixel values (intensity)
        sharpness: amount of detail
        saturation: color intensity

    Returns:
       a transformed data array better for RGB display
    """
    rgb = rgb.where(rgb > 0)
    rgb = np.log(rgb / scale) / np.log(1 / scale)
    rgb = rgb.where(rgb >= vmin, vmin)
    rgb = rgb.where(rgb <= vmax, vmax)    
    rgb_min = rgb.min(("number_of_lines", "pixels_per_line"))
    rgb_max = rgb.max(("number_of_lines", "pixels_per_line"))
    rgb = (rgb - rgb_min) / (rgb_max - rgb_min)
    rgb = rgb * gamma
    img = rgb * 255
    img = img.where(img.notnull(), 0).astype("uint8")
    img = Image.fromarray(img.data)
    enhancer = ImageEnhance.Contrast(img)
    img = enhancer.enhance(contrast)
    enhancer = ImageEnhance.Brightness(img)
    img = enhancer.enhance(brightness)
    enhancer = ImageEnhance.Sharpness(img)
    img = enhancer.enhance(sharpness)
    enhancer = ImageEnhance.Color(img)
    img = enhancer.enhance(saturation)
    rgb[:] = np.array(img) / 255
    return rgb

def enhancel3(rgb, scale = .01, vmin = 0.01, vmax = 1.02, gamma=.95, contrast=1.5, brightness=1.02, sharpness=2, saturation=1.1):
   
    rgb = rgb.where(rgb > 0)
    rgb = np.log(rgb / scale) / np.log(1 / scale)
    rgb = (rgb -  rgb.min()) / (rgb.max() - rgb.min())
    rgb = rgb * gamma
    img = rgb * 255
    img = img.where(img.notnull(), 0).astype("uint8")
    img = Image.fromarray(img.data)
    enhancer = ImageEnhance.Contrast(img)
    img = enhancer.enhance(contrast)
    enhancer = ImageEnhance.Brightness(img)
    img = enhancer.enhance(brightness)
    enhancer = ImageEnhance.Sharpness(img)
    img = enhancer.enhance(sharpness)
    enhancer = ImageEnhance.Color(img)
    img = enhancer.enhance(saturation)
    rgb[:] = np.array(img) / 255    
    return rgb

def pcolormesh(rgb):
    fig = plt.figure()
    axes = plt.subplot()
    artist = axes.pcolormesh(
        rgb["longitude"],
        rgb["latitude"],
        rgb,
        shading="nearest",
        rasterized=True,
    )
    axes.set_aspect("equal")

back to top

2. Easy Global Chlorophyll-a Map#

Let’s start with the most basic visualization. First, get level-3 chlorophyll map product, which is already gridded on latitudes and longitudes.

results = earthaccess.search_data(
    short_name="PACE_OCI_L3M_CHL_NRT",
    granule_name="*.MO.*.0p1deg.*",
)
paths = earthaccess.open(results)

dataset = xr.open_dataset(paths[-1])
chla = dataset["chlor_a"]

Now we can already create a global map. Let’s use a special colormap from cmocean to be fancy, but this is not necessary to make a basic map. Use the robust="true" option to remove outliers.

artist = chla.plot(cmap=cmocean.cm.algae, robust="true")
plt.gca().set_aspect("equal")

../../_images/b301c115a0b8bb2c9567a280da8332d4b26a207a0ba13dfc1141221f5adcfa48.png

5. False Color for Ice Clouds#

We can use the same RGB image method, this time with different bands, to create false-color images that highlight specific elements in an image.

For example, using a combination of infrared bands can highlight ice clouds in the atmosphere, versus regular water vapor clouds. Let’s try it:

rhos_ice = dataset["rhos"].sel({"wavelength_3d": [1618, 2131, 2258]})

rhos_ice_enhanced = enhance(rhos_ice, vmin=0, vmax=0.9, scale=.1, gamma =1, contrast=1, brightness=1, saturation=1)
pcolormesh(rhos_ice_enhanced)

/tmp/ipykernel_246/1825938188.py:63: UserWarning: The input coordinates to pcolormesh are interpreted as cell centers, but are not monotonically increasing or decreasing. This may lead to incorrectly calculated cell edges, in which case, please supply explicit cell edges to pcolormesh.
  artist = axes.pcolormesh(

../../_images/891b4d2df24fc28526b3afc16fc1e5fb0af8364295192fad548521174b530f17.png

Here, the ice clouds are purple and water vapor clouds are white, like we can see in the northwestern region of the scene.

back to top

7. Full Spectra from Global Oceans#

This section contains interactive plots that require a live Python kernel. Not all the widgets will work when viewing the notebook; the cells have to be executed and the kernel still running. Download the notebook and give it a try!

The holoview library and its bokeh extension allow us to explore datasets interactively.

hv.extension("bokeh")

Let’s open a level 3 Rrs map.

results = earthaccess.search_data(
    short_name="PACE_OCI_L3M_RRS_NRT",
    granule_name="*.MO.*.1deg.*",
)
paths = earthaccess.open(results)

We can create a map from a single band in the dataset and see the Rrs value by hovering over the map.

dataset = xr.open_dataset(paths[-1])
def single_band(w):
    array = dataset.sel({"wavelength": w})
    return hv.Image(array, kdims=["lon", "lat"], vdims=["Rrs"]).opts(aspect="equal", frame_width=450, frame_height=250, tools=['hover'])

single_band(368)

In order to explore the hyperspectral datasets from PACE, we can have a look at the Rrs Spectrum at a certain location.

def spectrum(x, y):
    array = dataset.sel({"lon": x, "lat": y}, method="nearest")
    return hv.Curve(array, kdims=["wavelength"]).redim.range(Rrs=(-0.01, 0.04))

spectrum(0, 0)

Now we can build a truly interactive map with a slider to change the mapped band and a capability to show the spectrum where we click on the map. First, we build the slider and its interactivity with the map:

slider = pnw.DiscreteSlider(name="wavelength", options=list(dataset["wavelength"].data))
band_dmap = hv.DynamicMap(single_band, streams={"w": slider.param.value})

Then we build the capability to show the Rrs spectrum at the point where we click on the map.

points = hv.Points([])
tap = Tap(source=points, x=0, y=0)
spectrum_dmap = hv.DynamicMap(spectrum, streams=[tap])

Now try it!

slider

(band_dmap * points + spectrum_dmap).opts(shared_axes=False)

back to top

Satellite Data Visualization

Contents

Satellite Data Visualization#

Summary#

Learning Objectives#

Contents#

1. Setup#

2. Easy Global Chlorophyll-a Map#

3. Global Oceans in Quasi True Color#

4. Complete Scene in True Color#

5. False Color for Ice Clouds#

6. Phytoplankton in False Color#

7. Full Spectra from Global Oceans#

8. Animation from Multiple Angles#

Understanding Multi-Angle Data#

Radiance to Reflectance#

Animating an Overpass#