# HAMP comparison#

This notebook shows a comparison of cloud features as detected by different parts of the Halo Microwave Package (HAMP) and includes data from WALES and specMACS for additional context. Another application of the HAMP data is included in the chapter HALO UNIFIED dataset.

import numpy as np
import eurec4a
cat = eurec4a.get_intake_catalog(use_ipfs="QmahMN2wgPauHYkkiTGoG2TpPBmj3p5FoYJAq9uE9iXT9N")


We’ll define a helper-function to approximately convert cell centers to edges:

def center_to_edge(center):
center = np.asarray(center)
edge =  np.empty(center.shape[0]+1, dtype=center.dtype)
# the strange notation in the following lines is used to make the calculations datetime compatible
edge[1:-1] = center[1:] + (center[:-1] - center[1:])/2.
edge[0] = edge[1] + (edge[1] - edge[2])
edge[-1] = edge[-2] + (edge[-2] - edge[-3])
return edge


And define a multi-linear scale that we will use for the LWP plot later on. This scale allows us to use zoom in a certain y-range. In our case, we will later increase the scale from 20 to 20 by a factor of 10.

import matplotlib as mpl
from matplotlib.scale import ScaleBase, register_scale
from matplotlib.transforms import Transform
from matplotlib.ticker import (
NullFormatter, ScalarFormatter,
NullLocator, AutoLocator, AutoMinorLocator,
)

class MultiLinearTransform(Transform):
input_dims = output_dims = 1
is_separable = True

def __init__(self, segments):
super().__init__(self)
self.segments = segments
self.inverse_segments = [
(
self.transform_non_affine(np.array([p0]))[0],
self.transform_non_affine(np.array([p1]))[0],
1.0/scale
)
for (p0, p1, scale) in self.segments
]

def transform_non_affine(self, a):
result = a.copy().astype(float)
with np.errstate(invalid="ignore"):
for (p0, p1, scale) in self.segments:
result[p1<a] += (p1 - p0) * (scale - 1)
return result

def inverted(self):
return MultiLinearTransform(self.inverse_segments)

class MultiLinearScale(ScaleBase):
name = "multilinear"

def __init__(self, axis, segments):
"""
A MPL ax scale that magnifies given segments

Parameters
----------
axis : {matplotlib axis}
The axis handle
segments : {[tuples of 3-tuples]}
Format: [(p0, p1, scale), ...]
Magnify each segment from p0 to p1 by given scale factor.
p0 <= p1

Example
-------
ax.set_yscale('magnifylinear', segments=[[-20., 20., 10.]])

"""
for (p0, p1, scale) in segments:
if p0 > p1:
raise ValueError("p0 must be less or equal p1")
if scale <= 0:
raise ValueError("scale must be positive")

super().__init__(axis)
self.segments = segments

def get_transform(self):
return MultiLinearTransform(self.segments)

def set_default_locators_and_formatters(self, axis):
# docstring inherited
axis.set_major_locator(AutoLocator())
axis.set_major_formatter(ScalarFormatter())
axis.set_minor_formatter(NullFormatter())
# update the minor locator for x and y axis based on rcParams
if (axis.axis_name == 'x' and mpl.rcParams['xtick.minor.visible'] or
axis.axis_name == 'y' and mpl.rcParams['ytick.minor.visible']):
axis.set_minor_locator(AutoMinorLocator())
else:
axis.set_minor_locator(NullLocator())

# register scale in matplotlib
register_scale(MultiLinearScale)


## Collect datasets#

We want to look at data from several sensors within a spacific timeframe on flight HALO-0205:

flight_id = 'HALO-0205'
timeslice = slice('2020-02-05T13:07', '2020-02-05T13:11')



Ideally, we’d be handling the cloud flag meanings properly, but for now, let’s just check if they are defined similarly:

assert np.all(radar_cm.cloud_mask.flag_values == radiometer_cm.cloud_mask.flag_values)
assert wales.cloud_mask.flag_meanings == 'cloud_free probably_cloudy most_likely_cloudy'


Let’s prepare the radar reflectivity (dBZ) time-height plot. Get the dBZ in right shape and adjust the x and y axis according to the requirements of pcolormesh, this means x and y define the pixel edges in the respective direction. Also we convert the height coordinate from the WGS84 ellipsoid height to height above the Geoid. For this, we need y (and x) as 2D fields, as we have to correct the height in every time step as function of lat and lon. The Radar and Lidar cloud top heigh products are already defined above the Geoid or sea surface.

def wgs84_height(lon, lat):
#TODO: find wgs84_height(lon, lat) function that works with the licence
# this is a good average number for the EUREC4A circle area.
return np.zeros_like(lon+lat) - 47.5



## Plot timeseries#

%matplotlib inline
import matplotlib.pyplot as plt
plt.style.use("./mplstyle/book")

fig, (ax3, ax2, ax1, ax,) = plt.subplots(
nrows=4, sharex=True,
gridspec_kw=dict(height_ratios=[1, 1, .75, 0.5]),
)

ax.plot(wales.time, wales.cloud_mask + 0.2, '.', label='WALES', markersize=3)

# Make nice labels
ax.legend(ncol=3, loc='lower right')
ax.grid()

ax2.plot(
wales.time,  wales.cloud_top,
'.', markersize=3,
label='WALES cloud top'
)
ax2.plot(
'.', markersize=3,
)
ax2.set_ylabel('Height above Geoid (m)')
ax2.set_ylim(0, 2500)
ax2.grid()
ax2.legend(ncol=3)
ax2.minorticks_on()

# specMACS
ax3.set_title('specMACS')

# Plot MWR Retrievals