Advanced visualisation

Advanced visualisation#

Overview

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

Prerequisites

Basic understanding of Eradiate’s angular conventions (see Conventions in Eradiate for details).
Basic knowledge of xarray’s data structures (see the xarray documentation for an introduction).

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
import eradiate as ert
ert.set_mode("mono")

def run_exp(
    film_resolution=(32, 32),
    illumination_convention="east_right",
    measure_convention="east_right",
):
    exp = ert.experiments.AtmosphereExperiment(
        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()

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()

[5]:

<matplotlib.image.AxesImage object at 0x7cfacc4f6a00>

../../_images/tutorials_howto_advanced_visualisation_9_3.png

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(
        [
            ["radiance", "Radiance", "viridis"],
            ["vza", "Zenith", "cividis"],
            ["vaa", "Azimuth", "cividis"],
        ],
        axs,
    ):
        divider = make_axes_locatable(ax)
        cax = divider.append_axes("right", size="5%", pad=0.1)
        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)

../../_images/tutorials_howto_advanced_visualisation_11_1.png

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.

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):
    radiances = ds.radiance.values.ravel()
    zeniths = ds.vza.values.ravel()  # Degree
    azimuths = np.deg2rad(ds.vaa.values).ravel()  # Radian

    # 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 = fig.add_axes(rect, aspect="equal")
    ax_cartesian.axis("off")  # Hide axis
    ctr = ax_cartesian.tricontourf(triangles, radiances, levels=levels, cmap="viridis")

    if show_contour:
        ax_cartesian.tricontour(triangles, radiances, levels=16, linewidths=0.5, colors="k")

    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 = fig.add_axes(rect, polar=True, facecolor="none")
    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

    # Add the color bar (important: both axes must be adjusted)
    fig.colorbar(ctr, ax=[ax_cartesian, ax_polar])

    plt.show()
    plt.close()

plot_polarfilm(result)

../../_images/tutorials_howto_advanced_visualisation_21_1.png

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

[11]:

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

../../_images/tutorials_howto_advanced_visualisation_23_1.png

And we can also add the stored azimuth values:

[12]:

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

../../_images/tutorials_howto_advanced_visualisation_25_1.png

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):
    radiances = ds.radiance.values.ravel()
    zeniths = ds.vza.values.ravel()  # Degree
    azimuths = ert.frame.transform_azimuth(
        np.deg2rad(ds.vaa.values).ravel(), to_convention="north_left"
    )  # Radian

    # 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 = fig.add_axes(rect, aspect="equal")
    ax_cartesian.axis("off")  # Hide axis
    ctr = ax_cartesian.tricontourf(triangles, radiances, levels=levels, cmap="viridis")

    if show_contour:
        ax_cartesian.tricontour(triangles, radiances, levels=16, linewidths=0.5, colors="k")

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

    ## Polar axes
    ax_polar = fig.add_axes(rect, polar=True, facecolor="none")
    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)

    # Add the color bar (important: both axes must be adjusted)
    fig.colorbar(ctr, ax=[ax_cartesian, ax_polar])

    plt.show()
    plt.close()

plot_polarfilm_north_left(result)

../../_images/tutorials_howto_advanced_visualisation_27_1.png

Final words#

Plotting data is an important part of the scientific workflow and doing so in polar coordinates is a challenge. This tutorial provided you with some guidance, but you will very likely encounter situations where you will have to go beyond. Good luck!

Advanced visualisation

Contents

Advanced visualisation#

Inspecting and preparing the dataset#

Visualising the film data#

Gridded polar plots#

Polar film plots#

Final words#