I am having an issue with matching up the color table/brightness on CMI01 through CMI06 when creating GOES16 imagery with MetPy. I've tried using stock color tables and using random vmin/vmax to try and get a match. I've also tried using custom made color tables and even tried integrating things like min_reflectance_factor && max_reflectance_factor as vmin/vmax values.
Maybe I'm making this way more difficult than it is? Is there something I'm missing? Below are excerpts of code helping to create the current image output that I have:
grayscale = {"colors": [(0,0,0),(0,0,0),(255,255,255),(255,255,255)], "position": [0, 0.0909, 0.74242, 1]}
CMI_C02 = {"name": "C02", "commonName": "Visible Red Band", "grayscale": True, "baseDir": "visRed", "colorMap": grayscale}
dat = data.metpy.parse_cf('CMI_'+singleChannel['name'])
proj = dat.metpy.cartopy_crs
maxConcat = "max_reflectance_factor_"+singleChannel['name']
vmax = data[maxConcat]
sat = ax.pcolormesh(x, y, dat, cmap=make_cmap(singleChannel['colorMap']['colors'], position=singleChannel['colorMap']['position'], bit=True), transform=proj, vmin=0, vmax=vmax)
make_cmap
is a handy dandy method I found that helps to create custom color tables. This code is part of a multiprocessing process, so singleChannel
is actually CMI_C02
.
For reference, the first image is from College of DuPage and the second is my output...
So your problem is, I believe, because there's a non-linear transformation being applied to the data on College of DuPage, in this case a square root (sqrt
). This has been applied to GOES imagery in the past, as mentioned in the GOES ABI documentation. I think that's what is being done by CoD.
Here's a script to compare with and without sqrt:
import cartopy.feature as cfeature
from datetime import datetime, timedelta
import matplotlib.pyplot as plt
import metpy
import numpy as np
from siphon.catalog import TDSCatalog
# Trying to find the most recent image from around ~18Z
cat = TDSCatalog('http://thredds.ucar.edu/thredds/catalog/satellite/goes16'
'/GOES16/CONUS/Channel02/current/catalog.xml')
best_time = datetime.utcnow().replace(hour=18, minute=0, second=0, microsecond=0)
if best_time > datetime.utcnow():
best_time -= timedelta(days=1)
ds = cat.datasets.filter_time_nearest(best_time)
# Open with xarray and pull apart with some help using MetPy
data = ds.remote_access(use_xarray=True)
img_data = data.metpy.parse_cf('Sectorized_CMI')
x = img_data.metpy.x
y = img_data.metpy.y
# Create a two panel figure: one with no enhancement, one using sqrt()
fig = plt.figure(figsize=(10, 15))
for panel, func in enumerate([None, np.sqrt]):
if func is not None:
plot_data = func(img_data)
title = 'Sqrt Enhancement'
else:
plot_data = img_data
title = 'No Enhancement'
ax = fig.add_subplot(2, 1, panel + 1, projection=img_data.metpy.cartopy_crs)
ax.imshow(plot_data, extent=(x[0], x[-1], y[-1], y[0]),
cmap='Greys_r', origin='upper')
ax.add_feature(cfeature.COASTLINE, edgecolor='cyan')
ax.add_feature(cfeature.BORDERS, edgecolor='cyan')
ax.add_feature(cfeature.STATES, edgecolor='cyan')
ax.set_title(title)
Which results in:
The lower image, with the sqrt
transformation applied seems to match the CoD image pretty well.