Source code for eradiate.tutorials

"""
Extensions and helpers for tutorials.
"""

from __future__ import annotations

import numpy as np
import xarray as xr
from matplotlib import pyplot as plt
from numpy.typing import ArrayLike

from .config import settings
from .frame import AzimuthConvention
from .scenes.atmosphere import AbstractHeterogeneousAtmosphere
from .spectral import SpectralIndex


[docs] def plot_polarfilm( da: xr.DataArray, levels: int | ArrayLike = 16, show_contour: bool = True, show_azimuth: bool = False, figsize: tuple = (3, 3), azimuth_convention: AzimuthConvention | str | None = None, vmin: float | None = None, vmax: float | None = None, theta_max: float = 90.0, ): """ Plot film data generated by a hemispherical distant sensor as a polar contour plot. Parameters ---------- da : xarray.DataArray Data array to be plotted. levels : int or array-like, default: 16 Number of levels on the contour plot. show_contour : bool, default: True If ``True``, show contour outlines. show_azimuth : bool, default: False [Debug option] If ``True``, show azimuth values as points. figsize : tuple, default: (3, 3) Figure size. azimuth_convention : .AzimuthConvention or str, optional Azimuth convention of plotted data (will adjust the azimuth of the plot). vmin, vmax : float, optional Range of values to plot. If unset, automatically from the data. theta_max : float, default: 90.0 Maximum value of the azimuth angle on plots. Returns ------- matplotlib.Figure """ import matplotlib.pyplot as plt import matplotlib.tri as tri # Filter grazing angles da = ( da.where(da.vza <= theta_max) .dropna("x_index", how="all") .dropna("y_index", how="all") ) theta_max = float(da.vza.max()) rscale = 90.0 / theta_max # Prepare data for plotting azimuth_convention = ( AzimuthConvention.convert(azimuth_convention) if azimuth_convention is not None else settings.AZIMUTH_CONVENTION ) values = da.transpose("x_index", "y_index").values.ravel() zeniths = da["vza"].values.ravel() # Degree azimuths = np.deg2rad(da["vaa"].values).ravel() # Radian # Create triangulation x = zeniths * np.cos(azimuths) * rscale y = zeniths * np.sin(azimuths) * rscale triangles = tri.Triangulation(x, y) # Make plot fig = plt.figure(0, figsize=figsize) 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, values, levels=levels, cmap="viridis", vmin=vmin, vmax=vmax, ) if show_contour: ax_cartesian.tricontour( triangles, values, levels=levels, linewidths=0.5, colors="w", alpha=0.5, vmin=vmin, vmax=vmax, ) 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, theta_max]) # Cover the full zenith value range ax_polar.grid(False, axis="x") ax_polar.grid(axis="y", alpha=0.25) ax_polar.tick_params(axis="y", colors="white") ax_polar.set_rlabel_position(120) ax_polar.set_theta_offset(azimuth_convention.value[0]) ax_polar.set_theta_direction(azimuth_convention.value[1]) # Add the color bar (important: both axes must be adjusted) fig.colorbar(ctr, ax=[ax_cartesian, ax_polar], shrink=0.75, pad=0.15) return fig
[docs] def plot_sigma_t( *atmospheres: AbstractHeterogeneousAtmosphere, labels: list[str] | None = None, altitude_extent: tuple[float, float] | None = None, si: SpectralIndex | None = None, show: bool = True, ) -> tuple[plt.Figure, plt.Axes]: """ Display the extinction coefficient of one or several atmosphere objects for a single spectral context. Parameters ---------- *atmospheres : .AbstractHeterogeneousAtmosphere One or several atmosphere objects for which to plot the extinction coefficient. labels : list of strings, optional Labels associated with the passed atmosphere objects. If unset, no legend will be added to the plot. altitude_extent : tuple of float A (min, max) altitude pair (in km) to which the plot is restricted. si : .SpectralIndex, optional The spectral index at which the extinction coefficient is evaluated. If unset, a default spectral index is created using :meth:`.SpectralIndex.new`. show : bool, optional If ``True``, return ``None`` and display the plot. Otherwise, return a (``Figure``, ``Axes``) pair. """ from matplotlib.ticker import ScalarFormatter from eradiate.units import to_quantity if si is None: # If none is found, fall back to the default (550 nm) si = SpectralIndex.new() else: # Otherwise, apply conversion protocol si = SpectralIndex.convert(si) if labels is None: label_iter = iter([None for _ in atmospheres]) else: label_iter = iter(labels) fig, ax = plt.subplots(1, 1) with plt.rc_context({"lines.linestyle": ":", "lines.marker": "."}): for atmosphere in atmospheres: altitude = to_quantity(atmosphere.eval_radprops(si=si).z_layer).m_as("km") sigma_t = to_quantity(atmosphere.eval_radprops(si=si).sigma_t).m_as("1/m") ax.plot(altitude, sigma_t, label=next(label_iter)) formatter = ScalarFormatter(useMathText=True) formatter.set_powerlimits((-3, 2)) ax.yaxis.set_major_formatter(formatter) ax.set_ylabel("Extinction coefficient [1/m]") if labels is not None: fig.legend( bbox_to_anchor=(1.0, 0.5), loc="center left", borderaxespad=0.0, ) if altitude_extent is not None: ax.set_xlim(altitude_extent) ax.set_xlabel("Altitude [km]") plt.tight_layout() if show: plt.show() else: return fig, ax
[docs] def load_ipython_extension(ipython): """ IPython extension for Eradiate tutorials. Should be loaded at the top of the tutorial notebook. It sets the Matplotlib style and prints the current date and Eradiate version as markdown. See Also -------- :func:`eradiate.plot.set_style` Notes ----- This extension should be loaded using the IPython magic: .. code:: %load_ext eradiate.tutorials """ import datetime from IPython.display import Markdown, display from . import __version__ from .plot import set_style set_style() xr.set_options(display_expand_data=False) display( Markdown( f"*Last updated: {datetime.datetime.now():%Y-%m-%d %H:%M} " f"(eradiate v{__version__})*" ) )