Go to notebook file

[1]:

%reload_ext eradiate.notebook.tutorials


Last updated: 2022-10-05 13:34 (eradiate v0.22.4.post102+g84ef86b.d20221005)

Overview

In this tutorial, we introduce visualisation techniques for the data used and produced by Eradiate.

Prerequisites

What you will learn

How to visualise conveniently the hemispherical data used produced by Eradiate.

Let’s start by running a quick sample simulation. We use a hemispherical distant measure, which samples uniformly the hemisphere and stores radiance samples in a film data structure which maps the [0, 1]² space to the hemisphere. This coverage method is more efficient than uniform gridding of the (zenith, azimuth) space, but the output data can be tricky to manipulate and visualise.

[2]:

%load_ext eradiate
ert.set_mode("mono")

def run_exp(
film_resolution=(32, 32),
illumination_convention="east_right",
measure_convention="east_right",
):
exp = ert.experiments.OneDimExperiment(
surface={"type": "rpv"},
illumination={
"type": "directional",
"zenith": 30,
"azimuth": 90,
"azimuth_convention": illumination_convention,
},
measures={
"type": "hdistant",
"spp": 100000,
"film_resolution": film_resolution,
"azimuth_convention": measure_convention,
},
)
return ert.run(exp).squeeze(drop=True)  # Squeeze and drop scalar dimensions

result = run_exp()


## Inspecting and preparing the dataset#

The dataset itself has a comprehensive represention, which we can use to check what the data looks like. In the output of the following cell, pay attention to the Dimensions entry, then take a look at the dimensions associated with the vza, vaa and radiance variables.

[3]:

result

[3]:

<xarray.Dataset>
Dimensions:     (y_index: 32, x_index: 32)
Coordinates:
* y_index     (y_index) int64 0 1 2 3 4 5 6 7 8 ... 23 24 25 26 27 28 29 30 31
y           (y_index) float64 0.0 0.03226 0.06452 ... 0.9355 0.9677 1.0
* x_index     (x_index) int64 0 1 2 3 4 5 6 7 8 ... 23 24 25 26 27 28 29 30 31
x           (x_index) float64 0.0 0.03226 0.06452 ... 0.9355 0.9677 1.0
vza         (x_index, y_index) float64 86.47 86.47 86.47 ... 86.47 86.47
vaa         (x_index, y_index) float64 225.0 222.1 219.2 ... 39.19 42.1 45.0
Data variables:
radiance    (y_index, x_index) float64 0.2161 0.2131 ... 0.2257 0.2243
spp         int64 100000
srf         float64 1.0
brdf        (y_index, x_index) float64 0.1328 0.131 0.1315 ... 0.1387 0.1378
brf         (y_index, x_index) float64 0.417 0.4114 0.4132 ... 0.4357 0.433
Attributes:
convention:  CF-1.8
history:     2022-10-05 13:34:57 - data creation - AtmosphereExperiment.p...
references:
title:       Top-of-atmosphere simulation results

The ordering of the various dimensions is critical because it reflects how the data is laid out internally. Dimension ordering may be different from a variable to another; it is not automatically a problem, and xarray usually accomodates this well; it may however be a problem when flattening data buffers because variables with different dimension ordering will have different (and incompatible) data layouts. We can remedy that by forcing dimension ordering using Dataset.transpose() and force all variables to have the same dimension ordering. We arbitrarily decide to use a column-major dimension order (film width is the innermost dimension).

[4]:

result = result.transpose("y_index", "x_index")
result

[4]:

<xarray.Dataset>
Dimensions:     (y_index: 32, x_index: 32)
Coordinates:
* y_index     (y_index) int64 0 1 2 3 4 5 6 7 8 ... 23 24 25 26 27 28 29 30 31
y           (y_index) float64 0.0 0.03226 0.06452 ... 0.9355 0.9677 1.0
* x_index     (x_index) int64 0 1 2 3 4 5 6 7 8 ... 23 24 25 26 27 28 29 30 31
x           (x_index) float64 0.0 0.03226 0.06452 ... 0.9355 0.9677 1.0
vza         (y_index, x_index) float64 86.47 86.47 86.47 ... 86.47 86.47
vaa         (y_index, x_index) float64 225.0 227.9 230.8 ... 50.81 47.9 45.0
Data variables:
radiance    (y_index, x_index) float64 0.2161 0.2131 ... 0.2257 0.2243
spp         int64 100000
srf         float64 1.0
brdf        (y_index, x_index) float64 0.1328 0.131 0.1315 ... 0.1387 0.1378
brf         (y_index, x_index) float64 0.417 0.4114 0.4132 ... 0.4357 0.433
Attributes:
convention:  CF-1.8
history:     2022-10-05 13:34:57 - data creation - AtmosphereExperiment.p...
references:
title:       Top-of-atmosphere simulation results

## Visualising the film data#

Now that our data layout is nice and uniform, we can take a look at is the film data itself. The plotting facilities offered by xarray are very convenient to get a quick overview:

[5]:

result.radiance.squeeze().plot.imshow()


This is not a great plot, but it is already useful for a quick check. We can see that the hotspot we would expect from the RPV surface we used is correctly position (zenith 30°, azimuth 90° in the East right convention). The dimension ordering we have chosen is arbitrary, but not entirely random: it guarantees that xarray will map the film width to the x axis. A major issue with this plot is that it does not have an equal aspect ratio for the x and y axes: pixels are not squares.

Let’s build a cleaner plot. We will also add the viewing angles to see how they are mapped to the film.

[6]:

import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable

def plot_film(ds):
fig, axs = plt.subplots(1, 3, figsize=(4*3, 4))

for (var, title, cmap), ax in zip(
[
["vza", "Zenith", "cividis"],
["vaa", "Azimuth", "cividis"],
],
axs,
):
divider = make_axes_locatable(ax)
im = ax.imshow(
ds[var].values,
aspect="equal",  # Important to get square pixels
origin="lower",  # Position the origin as on the xarray plot
cmap=cmap,
)
fig.colorbar(im, cax=cax)
ax.set_title(title)

plt.tight_layout()
plt.show()
plt.close()

plot_film(result)


We enforced an equal aspect ratio for the x and y axes, and we can see that the viewing azimuth angles are mapped to film coordinates following with the East right convention.

## Gridded polar plots#

Film coordinates are not the canonical representation for hemispherical data in Earth observation. Instead, scientists usually prefer polar representations. We can make polar plots with our film data with a little bit of effort.

The first approach consists in gridding the data in the (zenith, azimuth) space. Eradiate includes the film_to_angular() function to automate this task. It is also accessible through the xarray accessor. In order to use this feature, we first have to make sure that the index coordinates for the film data are the film coordinates, and not pixel indices. Wethen use the accessor to grid the data on a 5°-spaced angular grid.

[7]:

import numpy as np

da = da.ert.to_angular(
)
da

[7]:

<xarray.DataArray (phi: 72, theta: 19)>
array([[0.17492927, 0.17486204, 0.17436774, ..., 0.18324441, 0.19365702,
0.2060571 ],
[0.17492927, 0.17530398, 0.17514415, ..., 0.18523042, 0.19542949,
0.20748091],
[0.17492927, 0.17574592, 0.17592056, ..., 0.18727253, 0.19709814,
0.20873257],
...,
[0.17492927, 0.17353621, 0.17201551, ..., 0.1781417 , 0.19038647,
0.20580826],
[0.17492927, 0.17397816, 0.17281492, ..., 0.17976667, 0.19145812,
0.20560341],
[0.17492927, 0.1744201 , 0.17359133, ..., 0.18137775, 0.19266521,
0.20614443]])
Coordinates:
* phi      (phi) float64 0.0 0.08727 0.1745 0.2618 ... 5.934 6.021 6.109 6.196
* theta    (theta) float64 0.0 0.08727 0.1745 0.2618 ... 1.396 1.484 1.571

We can now proceed with plotting using the pcolormesh() function.

[8]:

import numpy as np

def plot_polargrid(da):
da = da.transpose("theta", "phi")  # Ensure appropriate data layout
values = da.values

fig, ax = plt.subplots(1, 1, subplot_kw={"projection": "polar"})

ax.grid(False)  # Hide the grid
plt.colorbar()
ax.set_yticklabels([])  # No radial tick labels

plt.show()
plt.close()

plot_polargrid(da)


This method is costly because the gridding step interpolates the film data. However, it uses Matplotlib’s polar plotting features, which make it easy to customise the final plot. We can, for instance, change the azimuth convention very easily.

[9]:

def plot_polargrid_north_left(da):
da = da.transpose("theta", "phi")  # Ensure appropriate data layout
values = da.values

fig, ax = plt.subplots(1, 1, subplot_kw={"projection": "polar"})

ax.grid(False)  # Hide the grid
plt.colorbar()
ax.set_yticklabels([])  # No radial tick labels
ax.set_theta_direction(-1)
ax.set_theta_offset(np.pi / 2.0)

plt.show()
plt.close()

plot_polargrid_north_left(da)


## Polar film plots#

We do not necessarily had to grid the data to make a polar plot. Another approach is to use the tricontourf() function to make a contour plot with ungridded data. The following code optionally superimposes contour lines and azimuth values for checks.

[10]:

import matplotlib.pyplot as plt
import matplotlib.tri as tri

def plot_polarfilm(ds, levels=16, show_contour=True, show_azimuth=False):
zeniths = ds.vza.values.ravel()  # Degree

# Create triangulation
x = zeniths * np.cos(azimuths)
y = zeniths * np.sin(azimuths)
triangles = tri.Triangulation(x, y)

# Make plot
fig = plt.figure(0, figsize=(4, 4))
rect = [0, 0, 1, 1]

## Main plot in Cartesian coordinates
ax_cartesian.axis("off")  # Hide axis
ctr = ax_cartesian.tricontourf(triangles, radiances, levels=levels, cmap="viridis")

if show_contour:

if show_azimuth:
ax_cartesian.scatter(x, y, c=azimuths, cmap="plasma", s=3)

ax_cartesian.set_xlim([-90, 90])  ## Match limits with the full zenith range
ax_cartesian.set_ylim([-90, 90])

## Polar axes
ax_polar.set_rlim([0, 90]) # Cover the full zenith value range
ax_polar.grid(False)  # Hide the polar grid
ax_polar.set_yticklabels([])  # No radial tick labels

fig.colorbar(ctr, ax=[ax_cartesian, ax_polar])

plt.show()
plt.close()

plot_polarfilm(result)


We can refine the levels to get a smoother colour map:

[11]:

plot_polarfilm(result, levels=128, show_contour=False)


And we can also add the stored azimuth values:

[12]:

plot_polarfilm(result, levels=128, show_contour=False, show_azimuth=True)


Modifying the previous code to use the North left azimuth convention is similar to what we did for the gridded plot—with the important addition that angles must be transformed manually for the Cartesian part of the plot. The polar axes indeed display no data.

[13]:

def plot_polarfilm_north_left(ds, levels=16, show_contour=True, show_azimuth=False):
zeniths = ds.vza.values.ravel()  # Degree
azimuths = ert.frame.transform_azimuth(

# Create triangulation
x = zeniths * np.cos(azimuths)
y = zeniths * np.sin(azimuths)
triangles = tri.Triangulation(x, y)

# Make plot
fig = plt.figure(0, figsize=(4, 4))
rect = [0, 0, 1, 1]

## Main plot in Cartesian coordinates
ax_cartesian.axis("off")  # Hide axis
ctr = ax_cartesian.tricontourf(triangles, radiances, levels=levels, cmap="viridis")

if show_contour:

ax_cartesian.set_xlim([-90, 90])  ## Match limits with the full zenith range
ax_cartesian.set_ylim([-90, 90])

## Polar axes
ax_polar.set_rlim([0, 90]) # Cover the full zenith value range
ax_polar.grid(False)  # Hide the polar grid
ax_polar.set_yticklabels([])  # No radial tick labels
ax_polar.set_theta_direction(-1)
ax_polar.set_theta_offset(np.pi / 2.0)