MetPy Matching GOES16 Reflectance Brightness

Question

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...

Any help/guidance would be greatly appreciated!

Answer 1

After polling some meteorologists, I ended up making a color table that was in between the two images as the agreed general consensus was that they thought my version was too dark and the standard was too light.

I still used vmax and vmin for pcolormesh() and simplified my grayscale object to just two colors with a slightly darker gray than the standard.

Thanks to all who looked at this.

Answer 2

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.