I'm trying to reproject GOES16 Full Disk imagery using cartopy or pyproj. I'd like to get it into a different projection. For this example, I try to reproject it to Mercator. However, when I run the code I get a full globe image of the data not in Mercator projection and without any cartopy features. I feel like I'm missing a simple step but can't figure out just what it is. Below is a reproducible example - I'm using Python 3.5.
import matplotlib
import matplotlib.pyplot as plt
from siphon.catalog import TDSCatalog, get_latest_access_url
import numpy as np
from datetime import datetime, timedelta
import cartopy.crs as ccrs
import cartopy.feature as cfeature
# query data
nowdate = datetime.utcnow()
cat = TDSCatalog('http://thredds-jumbo.unidata.ucar.edu/thredds/catalog/satellite/goes16/GOES16/Products/SeaSurfaceTemperature/FullDisk/' + \
str(nowdate.year) + str("%02d"%nowdate.month) + str("%02d"%nowdate.day) + '/catalog.xml')
dataset_name = sorted(cat.datasets.keys())[-1]
dataset = cat.datasets[dataset_name]
# load netcdf and read variables
nc = dataset.remote_access()
sst = np.array(nc.variables['SST'][:,:])
sst[np.isnan(sst)] = -1
sst = np.ma.array(sst)
sst[sst < 0] = np.ma.masked
X = nc.variables['x'][:]
Y = nc.variables['y'][:]
# define projections
proj_var = nc.variables['goes_imager_projection']
globe = ccrs.Globe(ellipse='sphere', semimajor_axis=proj_var.semi_major_axis,
semiminor_axis=proj_var.semi_minor_axis)
# define reprojection target
proj = ccrs.Mercator(central_longitude=proj_var.longitude_of_projection_origin,
latitude_true_scale=proj_var.latitude_of_projection_origin,
globe=globe)
# Plot
fig = plt.figure(figsize=(10, 10))
ax = fig.add_subplot(1, 1, 1, projection=proj)
ax.coastlines(resolution='50m', color='black')
ax.add_feature(cfeature.STATES, linestyle=':', edgecolor='black')
ax.add_feature(cfeature.BORDERS, linewidth=2, edgecolor='black')
im = ax.imshow(sst, extent=(X.min(), X.max(), Y.min(), Y.max()), origin='upper')
# try again, this time with pyproj
from pyproj import Proj
p = Proj(proj='geos', h=proj_var.perspective_point_height, lon_0=proj_var.longitude_of_projection_origin, sweep=proj_var.sweep_angle_axis)
X = nc.variables['x'][:] * proj_var.perspective_point_height
Y = nc.variables['y'][:] * proj_var.perspective_point_height
XX, YY = np.meshgrid(X,Y)
lons, lats = p(XX, YY, inverse=True)
fig = plt.figure(figsize=(10, 10))
ax = fig.add_subplot(1, 1, 1, projection=proj)
ax.coastlines(resolution='50m', color='black')
ax.add_feature(cfeature.STATES, linestyle=':', edgecolor='black')
ax.add_feature(cfeature.BORDERS, linewidth=2, edgecolor='black')
im = ax.imshow(sst, extent=(lons.min(), lons.max(), lats.min(), lats.max()), origin='upper')
Your methodology was correct, but you have to use pcolormesh instead of imshow.
This should work:
from datetime import datetime
import cartopy.feature as cfeature
from siphon.catalog import TDSCatalog
import matplotlib.pyplot as plt
from matplotlib import patheffects
import metpy
from metpy.plots import colortables
import xarray as xr
from xarray.backends import NetCDF4DataStore
%matplotlib inline
nowdate = datetime.utcnow()
cat = TDSCatalog('http://thredds-jumbo.unidata.ucar.edu/thredds/catalog/satellite/goes16/GOES16/Products/SeaSurfaceTemperature/FullDisk/' + \
str(nowdate.year) + str("%02d"%nowdate.month) + str("%02d"%nowdate.day) + '/catalog.xml')
dataset_name = sorted(cat.datasets.keys())[-1]
dataset = cat.datasets[dataset_name]
ds = dataset.remote_access(service='OPENDAP')
ds = NetCDF4DataStore(ds)
ds = xr.open_dataset(ds)
dqf = ds.metpy.parse_cf('DQF')
dat = ds.metpy.parse_cf('SST')
proj = dat.metpy.cartopy_crs
dat = dat.where(dqf == 0)
dat = dat.where(dat.variable > 274)
dat = dat.where(dat.variable < 310)
dat = dat - 273.15
# Plot in Mercator
import cartopy.crs as ccrs
newproj = ccrs.Mercator()
fig = plt.figure(figsize=[12, 12], dpi=100)
ax = fig.add_subplot(1,1,1, projection=newproj)
im = ax.pcolormesh(dat['x'], dat['y'], dat, cmap='jet', transform=proj, vmin=-2, vmax=38)
ax.set_extent((dat['x'].min() + 4000000, dat['x'].max()- 3200000, dat['y'].min()+ 5500000, dat['y'].max()- 650000), crs=proj)