Displaying GOES in mercator projection

[blaylockbk]

For posterity, this question was asked on the Unidata python-users email list ([email protected])

I am trying to display GOES True Color imagery in a Mercator projection. I thought I would do this with a line such as:

mc=cartopy.ccrs.Mercator(central_longitude=0.0, min_latitude=-80.0, max_latitude=84.0, globe=None)

ax = fig.add_subplot(1, 1, 1, projection=mc)

ax.set_extent([-135, -60, 10, 65], crs=ccrs.Mercator())

but this does not work.

Can someone point out what I am doing wrong or point me to the correct resource so I might understand my problem better?

---

Code snippet from the GOES-2-go True Color recipe:

fig = plt.figure(figsize=(15, 12))

# Generate an Cartopy projection
lc = ccrs.LambertConformal(central_longitude=-97.5, standard_parallels=(38.5, 38.5)) 

ax = fig.add_subplot(1, 1, 1, projection=lc)
ax.set_extent([-135, -60, 10, 65], crs=ccrs.PlateCarree())

ax.imshow(
    RGB,
    origin='upper',
    extent=(x.min(), x.max(), y.min(), y.max()),
    transform=geos,
    interpolation='none'
)

ax.coastlines(resolution='50m', color='black', linewidth=0.5)
ax.add_feature(ccrs.cartopy.feature.STATES, linewidth=0.5)
 
plt.title('GOES-16 True Color', loc='left', fontweight='bold', fontsize=15)
plt.title('{}'.format(scan_start.strftime('%d %B %Y %H:%M UTC ')), loc='right')

[blaylockbk]

Answered by Ryan May

So when you type:

ax.set_extent([-135, -60, 10, 65], crs=ccrs.Mercator())

You are telling cartopy "set the extents of the axes to a box with lower left corner (-135, 10) and upper right corner (-60, 65); those points are native coordinates in the Mercator() projection". You essentially specified a 75m x 55m box. To plot with Mercator you want:

# Tell cartopy the axes we're plotting on should be in Mercator
ax = fig.add_subplot(1, 1, 1, projection=Mercator())
# Specify the extent as lat/lon
ax.set_extent([-135, -60, 10, 65], crs=ccrs.PlateCarree())

[csteele2]

Hmmm. When I use an axes with ANY map projection other than the native GOES projection, I only get all the map layers I add after, and either just a blank satellite image over the region, or a severely displaced image when using goes2go.

[csteele2]

So, this happens in multiple environments for me, at least one of which was setup from scratch with goes2go. The dataset that is returned has x and y values that look like the scan angle, but the rgb looks like it has those multiplied by the satellite height. Still not exactly sure how it can display just fine using the returned geostationary crs, or even making a geostationary crs from scratch in cartopy, but changing the axes projection to just about anything else doesn't work.

[msw17002]

It doesn't look like this was completed.

For the record, this works for me:

#===mercator projection with lon/lat bounds of N,E,S,W: (43.16672846653851, -117.33568318037872, 37.3774715334615, -124.91871681962125)
lc = crs.Mercator(central_longitude=(W+E)/2, min_latitude=S, max_latitude=N)
#===convert bounds from rads to distance in meters...
#        ds was imported as, **_ds = xr.open_dataset(path_to_goes_nc_file,engine='netcdf4')_**
ps = ds["goes_imager_projection"]
xs = ds["x"][:] * ps.perspective_point_height
ys = ds["y"][:] * ps.perspective_point_height
#===generate figure
fig = plt.figure(figsize=(8, 9))
#---generate axis
ax = fig.add_subplot(1, 1, 1, projection=lc)
ax.set_extent([W, E, S, N], crs=crs.PlateCarree())
#---plot goes
RGB = ds.rgb.TrueColor(gamma=2.2, pseudoGreen=True, night_IR=True)
ax.imshow(RGB,origin='upper',extent=(xs.min(), xs.max(), ys.min(), ys.max()),transform=ds.rgb.crs,interpolation='none')
##---add mapping resources of the USA and its counties...
#ax.add_geometries(fname_cou['geometry'], crs.PlateCarree(), facecolor="none", edgecolor='silver', lw=0.2)
#ax.add_geometries(fname_usa['geometry'], crs.PlateCarree(), facecolor="none", edgecolor='silver', lw=1.0)
#---add titles
ax.set_title(f"{ds.orbital_slot} {product}", loc="left", fontweight="bold")
ax.set_title(f"{ds.t.dt.strftime('%H:%M UTC %d-%b-%Y').item()}", loc="right")
##---save and close figure
#plt.savefig(outimg, bbox_inches="tight")
#plt.close(fig=None)
plt.show()

Hope this helps!

[csteele2]

I swore I tried both multiplying and dividing that before, both to no avail, but come to think of it, I think I was creating a new ds. This works, and since it does, doesn't that mean that **rgbkwargs aren't really functioning as intended?

[uni-byte]

Hi, I am trying the exact same method above but I am getting the following:

ax.imshow(RGB, origin='upper', extent=(xs.min(), xs.max(), ys.min(), ys.max()), transform=G.rgb.crs, interpolation='none')
ValueError: dimensions () must have the same length as the number of data dimensions, ndim=1

It appears to have something to with the extent, but when I leave it as default, the plot is shown blank. I dont really know how to fix this, can anyone help?