Reprojecting Satellite Image

[sgdecker]

Changing map projections always trips me up, but I thought this would work:

import xarray as xr
import metpy
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt

sat_file = '/nfs/wxdata1/ldmdata/GOES/GOES18/Products/CloudAndMoistureImagery/FullDisk/Channel09/20220815/OR_ABI-L2-CMIPF-M6C09_G18_s20222272030200_e20222272030200_c20222272030200.nc'

ds = xr.open_dataset(sat_file)

ds = ds.coarsen({'x': 2, 'y': 2}, boundary='trim').mean()

dat = ds.metpy.parse_cf('Sectorized_CMI')
proj = dat.metpy.cartopy_crs
x = dat['x']
y = dat['y']

fig = plt.figure()
ax = fig.add_subplot(1, 1, 1, projection=proj)
ax.add_feature(cfeature.COASTLINE.with_scale('50m'), edgecolor='cyan')
ax.add_feature(cfeature.STATES.with_scale('50m'), linestyle=':', edgecolor='cyan')
ax.add_feature(cfeature.BORDERS.with_scale('50m'), edgecolor='cyan')
ax.pcolormesh(x, y, dat, transform=proj)

plt.show()


map_proj = ccrs.Mercator()
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1, projection=map_proj)
ax.add_feature(cfeature.COASTLINE.with_scale('50m'), edgecolor='cyan')
ax.add_feature(cfeature.STATES.with_scale('50m'), linestyle=':', edgecolor='cyan')
ax.add_feature(cfeature.BORDERS.with_scale('50m'), edgecolor='cyan')
ax.set_extent((-159,-112.2,12.5,60.5))
ax.pcolormesh(x, y, dat, transform=proj)

plt.show()

The first plot does give me a beautiful full disk image, but my attempt to zoom in and display the plot on a Mercator projection does not. I've attached a screenshot showing the results. Is this a MetPy deficiency, or did I miss something?

[dopplershift]

No MetPy issue here, since it's only involved in the parsing of the netCDF metadata into the cartopy CRS.

Your problem is this call:

map_proj = ccrs.Mercator()

That's not specifying a central_longitude, so it uses the default of 0. My guess is using a projection so far from the region of interest/valid range of data causes some weird "seam-crossing" issue. If I instead use:

map_proj = ccrs.Mercator(central_longitude=-130)

I get the expected image: