I'm trying to plot GOES-East full disk data using metpy, and Siphon to download the latest data from the THREDDS data server. However, after comparing my plots with the realtime imagery, ther seems to be a large difference.
Below is my code:
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import metpy.calc as mpcalc
from metpy.plots.ctables import registry
from metpy.plots import add_timestamp
from metpy.units import units
from siphon.catalog import TDSCatalog
import xarray as xr
import numpy as np
from xarray.backends import NetCDF4DataStore
from datetime import datetime, timedelta
dt = datetime.utcnow().date()
data = TDSCatalog(f'http://thredds.ucar.edu/thredds/catalog/satellite/goes/east/products/'
f'CloudAndMoistureImagery/FullDisk/Channel09/{dt:%Y%m%d}/catalog.xml')
sat_dataset = data.datasets[0].remote_access(use_xarray = True)
cmi = sat_dataset.metpy.parse_cf('Sectorized_CMI')
x = cmi.coords['x'][:]
y = cmi.coords['y'][:]
timestamp = datetime.strptime(str(cmi.time.values.astype('datetime64[s]')), '%Y-%m-%dT%H:%M:%S')
print(timestamp)
vtime = timestamp.strftime('%Y-%m-%d %H%M%S')
# Create the figure
fig = plt.figure(figsize = [16, 10])
ax = fig.add_subplot(1, 1, 1, projection = cmi.metpy.cartopy_crs)
ax.set_extent([-80, -45, -50, -15], crs = ccrs.PlateCarree())
ax.add_feature(cfeature.BORDERS.with_scale('50m'), edgecolor = 'black', linewidth = 1)
ax.add_feature(cfeature.COASTLINE.with_scale('50m'), edgecolor = 'black', linewidth = 1)
ax.add_feature(cfeature.STATES.with_scale('50m'), edgecolor = 'white', linewidth = 1)
# Add mapping information
ax.add_feature(cfeature.STATES)
ax.add_feature(cfeature.BORDERS, linewidth=2)
# Plot the image with our colormapping choices
wv_norm, wv_cmap = registry.get_with_range('WVCIMSS_r', 193, 283)
im = ax.imshow(cmi, extent=(x[0], x[-1], y[0], y[-1]), origin='upper',
cmap = wv_cmap, norm = wv_norm, transform = cmi.metpy.cartopy_crs)
plt.colorbar(im, ticks = np.arange(193, 293, 10), ax = ax)
plt.title(f'Vapor da Água em Níveis Médios [$K$] \nValid: {vtime} UTC', loc = 'left')
plt.savefig(f'/mnt/c/Users/vitor/Desktop/WV_{vtime}.jpg', bbox_inches = 'tight')
Also below, is a comparison between the output from my code and the actual water vapor imagery from the CODNEXLAB website. I also looked at the metadata of the downloaded files and everything seems to be fine. Not sure if I'm doing something wrong here.
What you're seeing is that your image is flipped (it's easier to identify if you look at the global plot of that data). What's happening is the origin you specified ('upper'/'lower') disagree with what you passed as extent. So either tweak your origin
parameter:
im = ax.imshow(cmi, extent=(x[0], x[-1], y[0], y[-1]),
origin='lower', cmap=wv_cmap, norm=wv_norm,
transform=cmi.metpy.cartopy_crs)
or flip the order of your y extents:
im = ax.imshow(cmi, extent=(x[0], x[-1], y[-1], y[0]),
origin='upper', cmap=wv_cmap, norm=wv_norm,
transform=cmi.metpy.cartopy_crs)