Source code for typhon.plots.ppath

# -*- coding: utf-8 -*-

"""Functions to create plots using matplotlib.
"""

import numpy as np
import matplotlib.pyplot as plt


__all__ = [
    'plot_ppath',
    'plot_ppath_field',
    'ppath_field_minmax_posbox',
    'adjust_ppath_field_box_by_minmax',
    'plot_ppath_field_zenith_coverage_per_gp_p',
]


[docs] def plot_ppath(ppath, planetary_radius=1, lat_is_x=True, scale_alt=1000, ax=None): """ Return the ARTS Ppath plotted on the surface Parameters: ppath: (Ppath or Workspace) Propagation path or Wosrkspace with the path planetary_radius: (Numeric) Spherical radius of the planetary body (if ppath_field is not Workspace) lat_is_x: (Boolean) Have lat or lon act as the x-coordinate? scale_alt: (Numeric) Divide altitude and planetary_radius by this value ax: (AxesSubplot) axes to plot into. Assumes polar projection. If None, creates a plt.figure(projection='polar') and draws on its axis Example: >>> arts = typhon.arts.workspace.Workspace() >>> ################################################### >>> # Run ARTS simulation that sets an external ppath # >>> ################################################### >>> plot_ppath(arts) """ try: planetary_radius = ppath.refellipsoid.value[0] except: pass try: ppath = ppath.ppath.value except: pass alt, lat, lon, za, aa = ppath.alt_lat_lon_za_aa() if lat_is_x: x = lat else: x = lon if ax is None: ax = plt.subplot(111, projection='polar') ax.plot(np.deg2rad(x), alt/scale_alt) ax.set_rorigin(-planetary_radius/scale_alt) return ax
[docs] def plot_ppath_field(ppath_field, planetary_radius=1, lat_is_x=True, scale_alt=1000, subplots=1, auto_adjust=True, axes=None): """ Return the ARTS Ppath plotted on the surface Parameters: ppath_field: (Ppath or Workspace) Propagation path or Wosrkspace with the path planetary_radius: (Numeric) Spherical radius of the planetary body (if ppath_field is not Workspace) lat_is_x: (Boolean) Have lat or lon act as the x-coordinate? scale_alt: (Numeric) Divide altitude and planetary_radius by this value subplots: (Index) Divides the ppath_field into this many equally long segments, so that the first len(ppath_field) // subplots plots are drawn on the same surface and so forth auto_adjust: (Boolean) If true calls adjust_ppath_field_box_by_minmax axes: Either list of subplots, or None. Returns: list of axis of the length of subplots """ try: planetary_radius = ppath_field.refellipsoid.value[0] except: pass try: ppath_field = ppath_field.ppath_field.value except: pass N = len(ppath_field) if N % subplots != 0: raise RuntimeError("Bad inputs. Only for evenly divisible numbers") # Number of drawings per subplot n = N // subplots n1 = int(np.ceil(np.sqrt(subplots))) n2 = n1 if np.isclose(n1*n1, subplots) or n1*n1 > subplots else n1 + 1 if axes is None: fig = plt.gcf() axes = [] for i in range(subplots): axes.append(fig.add_subplot(n1, n2, i+1, projection='polar')) assert subplots == len(axes), "Must have same number of axes as subplots" for i in range(subplots): for j in range(n): plot_ppath(ppath_field[j+n*i], planetary_radius, lat_is_x, scale_alt, ax=axes[i]) if auto_adjust: adjust_ppath_field_box_by_minmax(axes, ppath_field, lat_is_x, scale_alt, 0.05) return axes
[docs] def ppath_field_minmax_posbox(ppath_field): """ Return the minimum and maximum of all pos variables of ppath_field Parameters: ppath_field: (Ppath or Workspace) Propagation path or Wosrkspace with the path Returns: min(alt), max(alt), min(lat), max(lat), min(lon), max(lon) """ try: ppath_field = ppath_field.ppath_field.value except: pass minalt = minlat = minlon = np.infty maxalt = maxlat = maxlon = -np.infty for ppath in ppath_field: alt, lat, lon, za, aa = ppath.alt_lat_lon_za_aa() minalt = min(alt.min(), minalt) maxalt = max(alt.max(), maxalt) minlat = min(lat.min(), minlat) maxlat = max(lat.max(), maxlat) minlon = min(lon.min(), minlon) maxlon = max(lon.max(), maxlon) return minalt, maxalt, minlat, maxlat, minlon, maxlon
[docs] def adjust_ppath_field_box_by_minmax(axes, ppath_field, lat_is_x=True, scale_alt=1000, extrap=0.05): """ Adjusts the axis of plotted ppath_field to Parameters: axes: Plotted axes of ppath_field ppath_field: (Ppath or Workspace) Propagation path or Wosrkspace with the path planetary_radius: (Numeric) Spherical radius of the planetary body (if ppath_field is not Workspace) lat_is_x: (Boolean) Have lat or lon act as the x-coordinate? scale_alt: (Numeric) Divide altitude and planetary_radius by this value extrap: (Numeric) How large an extrapolation of the box in terms of full posbox minmax """ minalt, maxalt, minlat, maxlat, minlon, maxlon = ppath_field_minmax_posbox(ppath_field) if lat_is_x: xmin = minlat xmax = maxlat max_diff = 180 else: xmin = minlon xmax = maxlon max_diff = 360 dx = (xmax-xmin) * extrap xmin -= dx xmax += dx if (xmax-xmin) > max_diff: return # no adjustment possible for ax in axes: ax.set_thetamin(xmin) ax.set_thetamax(xmax)
def wzeniths(zeniths): n = len(zeniths) if not n: return np.array([0, 180]), np.array([1, 1]) inds = np.argsort(zeniths) zaz = np.deg2rad(zeniths[inds]) cz = np.cos(zaz) wz = np.zeros((2*n)) za = np.zeros((2*n)) for i in range(n-1): N = i*2 za[N:N+2] = zaz[i:i+2] wz[0+N] = wz[1+N] = 0.5 * (cz[i] - cz[i+1]) return za, wz
[docs] def plot_ppath_field_zenith_coverage_per_gp_p(ppath_field, scale_alt=1000, axes=None): """Plots the zenith angle coverage of a ppath_field for all the altitudes in the field. Parameters: ppath_field: (Ppath or Workspace) Propagation path or Wosrkspace with the path axes: Either list of subplots, or None. Returns: list of axis, list of altitudes, list of maximum weights """ try: ppath_field = ppath_field.ppath_field.value except: pass alt = [] gps = [] for path in ppath_field: for i in range(path.np): if path.gp_p[i] not in gps: gps.append(path.gp_p[i]) alt.append(path.pos[i, 0]) alt = np.array(alt) # Size of subplots N = len(alt) zas = np.full((N), None) for path in ppath_field: for i in range(path.np): ip = np.where(path.pos[i, 0] == alt)[0][0] if zas[ip]: zas[ip].append(path.los[i, 0]) else: zas[ip] = [path.los[i, 0]] # Sub-plot grid n1 = int(np.ceil(np.sqrt(N))) n2 = n1 if np.isclose(n1*n1, N) or n1*n1 > N else n1 + 1 if axes is None: fig = plt.gcf() axes = [] for i in range(N): axes.append(fig.add_subplot(n1, n2, i+1, projection='polar')) assert N == len(axes), "Must have same number of axes as subplots" inds = np.argsort(alt) wei = [] alts = [] for i in range(N): pos = inds[i] za, wz = wzeniths(np.array(zas[pos])) wei.append(wz.max()) alts.append(alt[pos]) axes[i].plot(za, wz,'k') axes[i].set_rorigin(-0.25/np.pi) axes[i].set_thetamin(0) axes[i].set_thetamax(180) axes[i].set_rgrids([0, wz.max()*1.05], ['', '']) axes[i].set_thetagrids([0, 90, 180], ['Zenith', 'Limb', 'Nadir']) axes[i].set_title("Alt {} km".format(alt[pos]/scale_alt)) return axes, alts, wei