import matplotlib.pyplot as plt
import matplotlib.colors as pltcolors
from matplotlib.backends.backend_agg import FigureCanvasAgg as FigureCanvas
from matplotlib.figure import Figure
%config InlineBackend.figure_format='retina'
from colorspacious import cspace_convert
import numpy as np
import pyart
from skimage.io import imread
import skimage.color as skcolor
import xarray as xr
## You are using the Python ARM Radar Toolkit (Py-ART), an open source ## library for working with weather radar data. Py-ART is partly ## supported by the U.S. Department of Energy as part of the Atmospheric ## Radiation Measurement (ARM) Climate Research Facility, an Office of ## Science user facility. ## ## If you use this software to prepare a publication, please cite: ## ## JJ Helmus and SM Collis, JORS 2016, doi: 10.5334/jors.119
<frozen importlib._bootstrap>:283: DeprecationWarning: the load_module() method is deprecated and slated for removal in Python 3.12; use exec_module() instead
This colormap is built with four configurable segments:
Control points include:
# Control
min_value = -20
start_spectral = 10
start_pink = 60
start_purple = 70
linearize_greys_with_spectral = True
grey_point_lightness = 0
resolution_factor = 4
# Define initial colormaps
spectral_rgb_initial = plt.cm.Spectral_r(
np.linspace(0, 1, (start_pink - start_spectral) * resolution_factor)
)[:, :3]
pink_rgb_initial = plt.cm.PiYG(
np.linspace(0, .25, (start_purple - start_pink) * resolution_factor)
)[:, :3]
purple_rgb_initial = plt.cm.PRGn_r(
np.linspace(0.75, 1, (start_purple - start_pink) * resolution_factor)
)[:, :3]
spectral_lab = skcolor.rgb2lab(spectral_rgb_initial)
pink_lab = skcolor.rgb2lab(pink_rgb_initial)
purple_lab = skcolor.rgb2lab(purple_rgb_initial)
greyscale_lab = np.stack([
np.linspace(
grey_point_lightness,
spectral_lab[0, 0],
(start_spectral - min_value) * resolution_factor
),
np.linspace(0, spectral_lab[0, 1], (start_spectral - min_value) * resolution_factor),
np.linspace(0, spectral_lab[0, 2], (start_spectral - min_value) * resolution_factor)
], axis=-1)
# Piecewise linearization of lightness
peak_idx = np.argmax(spectral_lab[:, 0])
peak_value = np.linspace(
start_spectral, start_pink, (start_pink - start_spectral) * resolution_factor
)[peak_idx]
print(f"Spectral inflection point at value: {peak_value}")
if linearize_greys_with_spectral:
greys_and_spectral_increasing_lab = np.concatenate([greyscale_lab, spectral_lab[:peak_idx]])
greys_and_spectral_increasing_lab[:, 0] = np.linspace(
greys_and_spectral_increasing_lab[0, 0],
greys_and_spectral_increasing_lab[-1, 0],
greys_and_spectral_increasing_lab.shape[0]
)
else:
spectral_increasing_lab = spectral_lab[:peak_idx].copy()
spectral_increasing_lab[:, 0] = np.linspace(
spectral_increasing_lab[0, 0],
spectral_increasing_lab[-1, 0],
spectral_increasing_lab.shape[0]
)
greys_and_spectral_increasing_lab = np.concatenate([greyscale_lab, spectral_increasing_lab])
spectral_decreasing_lab = spectral_lab[peak_idx:].copy()
spectral_decreasing_lab[:, 0] = np.linspace(
spectral_decreasing_lab[0, 0],
spectral_decreasing_lab[-1, 0],
spectral_decreasing_lab.shape[0]
)
pink_lab[:, 0] = np.linspace(
spectral_decreasing_lab[-1, 0],
pink_lab[-1, 0],
pink_lab.shape[0]
)
purple_lab[:, 0] = np.linspace(
pink_lab[-1, 0],
purple_lab[-1, 0],
purple_lab.shape[0]
)
combined_lab = np.concatenate([
greys_and_spectral_increasing_lab,
spectral_decreasing_lab,
pink_lab,
purple_lab
])
combined_rgb = skcolor.lab2rgb(combined_lab)
# Save txt
np.savetxt("chase-spectral-202211091600.txt", combined_rgb)
# Colormap!
spectral_refl = pltcolors.LinearSegmentedColormap.from_list("ChaseSpectral", combined_rgb)
spectral_refl
Spectral inflection point at value: 34.87437185929648
Now, let's explore some visualizations!
def colormap_evaluation(sample, cmap, vmin=0, vmax=80):
# Sample image in original RGB
fig = Figure()
canvas = FigureCanvas(fig)
ax = fig.gca()
sample.plot.imshow(
vmin=vmin,
vmax=vmax,
cmap=cmap,
add_labels=False,
ax=ax
)
fig.tight_layout(pad=0)
ax.margins(0)
canvas.draw() # draw the canvas, cache the renderer
image_from_plot = np.frombuffer(canvas.tostring_rgb(), dtype=np.uint8)
image_from_plot = image_from_plot.reshape(canvas.get_width_height()[::-1] + (3,))
image_from_plot = image_from_plot / 255.
# Convert to greyscale
img_greyscale_JCh = cspace_convert(image_from_plot, "sRGB1", "JCh")
img_greyscale_JCh[..., 1] = 0
img_greyscale_sRGB = cspace_convert(img_greyscale_JCh, "JCh", "sRGB1")
# Deuteranomaly
cvd_space = {
"name": "sRGB1+CVD",
"cvd_type": "deuteranomaly",
"severity": 50
}
img_deuteranomaly_sRGB = cspace_convert(image_from_plot, cvd_space, "sRGB1")
# Deuteranopia
cvd_space = {
"name": "sRGB1+CVD",
"cvd_type": "deuteranomaly",
"severity": 100
}
img_deuteranopia_sRGB = cspace_convert(image_from_plot, cvd_space, "sRGB1")
# Protanomaly
cvd_space = {
"name": "sRGB1+CVD",
"cvd_type": "protanomaly",
"severity": 50
}
img_protanomaly_sRGB = cspace_convert(image_from_plot, cvd_space, "sRGB1")
# Deuteranopia
cvd_space = {
"name": "sRGB1+CVD",
"cvd_type": "protanomaly",
"severity": 100
}
img_protanopia_sRGB = cspace_convert(image_from_plot, cvd_space, "sRGB1")
# Overall Figure
fig, axes = plt.subplot_mosaic(
[
['l', 'l', 'l'],
['orig', 'd_anom', 'd_anop'],
['bw', 'p_anom', 'p_anop']
],
figsize=(12, 9)
)
# Color lightness scale
x = np.linspace(vmin, vmax, int((vmax - vmin) * 4))
color_range_sRGB = cmap(np.linspace(0, 1, len(x)))[:, :3]
color_range_Lab = cspace_convert(color_range_sRGB, "sRGB1", "CIELab")
axes['l'].scatter(
x,
color_range_Lab[:, 0],
c = [tuple(row) for row in color_range_sRGB],
s = 100
)
axes['l'].set_title(f"Radar colormap test for '{cmap.name}'")
axes['l'].set_ylabel("Lightness")
# Sample images
axes['orig'].imshow(image_from_plot)
axes['orig'].axis('off')
axes['orig'].set_title("Full Color")
axes['bw'].imshow(img_greyscale_sRGB)
axes['bw'].axis('off')
axes['bw'].set_title("Greyscale")
axes['d_anom'].imshow(img_deuteranomaly_sRGB)
axes['d_anom'].axis('off')
axes['d_anom'].set_title("Deuteranomaly")
axes['d_anop'].imshow(img_deuteranopia_sRGB)
axes['d_anop'].axis('off')
axes['d_anop'].set_title("Deuteranopia")
axes['p_anom'].imshow(img_protanomaly_sRGB)
axes['p_anom'].axis('off')
axes['p_anom'].set_title("Protanomaly")
axes['p_anop'].imshow(img_protanopia_sRGB)
axes['p_anop'].axis('off')
axes['p_anop'].set_title("Protanopia")
return fig
ds = refl = xr.open_zarr(
"s3://svrimgdm/v0/svrimg_dm_v0_testing_full.zarr/",
storage_options={
'anon': True,
'client_kwargs': {
'endpoint_url': 'https://sfo3.digitaloceanspaces.com'
}
},
consolidated=True
)
ds
/home/jthielen/mambaforge/envs/main/lib/python3.10/site-packages/cfgrib/xarray_plugin.py:11: DeprecationWarning: distutils Version classes are deprecated. Use packaging.version instead. if LooseVersion(xr.__version__) <= "0.17.0": /home/jthielen/mambaforge/envs/main/lib/python3.10/site-packages/setuptools/_distutils/version.py:346: DeprecationWarning: distutils Version classes are deprecated. Use packaging.version instead. other = LooseVersion(other)
<xarray.Dataset> Dimensions: (sample: 1100, time_offset: 5, y_index: 256, x_index: 256) Coordinates: (12/13) confidence (sample) int64 dask.array<chunksize=(1,), meta=np.ndarray> elapsed_time (sample) float64 dask.array<chunksize=(1,), meta=np.ndarray> mode (sample) <U2 dask.array<chunksize=(1,), meta=np.ndarray> * sample (sample) object 'spc_hail_2011_341389' ... 'spc_hail_20... source (sample) <U7 dask.array<chunksize=(1,), meta=np.ndarray> time_of_analysis (sample, time_offset) datetime64[ns] dask.array<chunksize=(1100, 5), meta=np.ndarray> ... ... * time_offset (time_offset) int64 -2 -1 0 1 2 user_id (sample) int64 dask.array<chunksize=(1100,), meta=np.ndarray> x_coord (sample, x_index) float32 dask.array<chunksize=(1, 256), meta=np.ndarray> * x_index (x_index) int64 0 1 2 3 4 5 6 ... 250 251 252 253 254 255 y_coord (sample, y_index) float32 dask.array<chunksize=(1, 256), meta=np.ndarray> * y_index (y_index) int64 0 1 2 3 4 5 6 ... 250 251 252 253 254 255 Data variables: echo_top_height (sample, time_offset, y_index, x_index) float32 dask.array<chunksize=(1, 1, 256, 256), meta=np.ndarray> projection int64 ... reflectivity (sample, time_offset, y_index, x_index) float32 dask.array<chunksize=(1, 1, 256, 256), meta=np.ndarray> Attributes: author: Jonathan Thielen (https://thielen.science) description: SVRIMG for Detailed Morphology, Version 0, Testing Set: Lab... reference: Thielen (2021): Two Tasks to Improve Mesoscale Weather Fore...
# Choose a mode, get a random sample with that mode label
mode = "BE"
idx = np.random.choice(np.flatnonzero(ds['mode'].values == mode))
sample = ds['reflectivity'].isel(sample=idx, time_offset=2).load()
sample
<xarray.DataArray 'reflectivity' (y_index: 256, x_index: 256)> array([[nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan], ..., [nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan]], dtype=float32) Coordinates: (12/13) confidence int64 5 elapsed_time float64 5.0 mode <U2 'BE' sample <U20 'spc_wind_2011_339632' source <U7 'myrorss' time_of_analysis datetime64[ns] 2011-08-12T01:00:00 ... ... time_offset int64 0 user_id int64 2 x_coord (x_index) float32 -2.62e+05 -2.6e+05 ... 2.46e+05 2.48e+05 * x_index (x_index) int64 0 1 2 3 4 5 6 ... 250 251 252 253 254 255 y_coord (y_index) float32 2.72e+05 2.74e+05 ... 7.8e+05 7.82e+05 * y_index (y_index) int64 0 1 2 3 4 5 6 ... 250 251 252 253 254 255 Attributes: grid_mapping: projection long_name: column_maximum_of_equivalent_reflectivity_factor units: dBZ
colormap_evaluation(sample, spectral_refl, vmin=min_value, vmax=2 * start_purple - start_pink)
plt.show()
/home/jthielen/mambaforge/envs/main/lib/python3.10/site-packages/matplotlib_inline/config.py:68: DeprecationWarning: InlineBackend._figure_format_changed is deprecated in traitlets 4.1: use @observe and @unobserve instead. def _figure_format_changed(self, name, old, new): Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers). Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers). Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers). Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers). Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).