I am trying desperately to project some geostationary data from GOES-16 netCDF file to a different projection. I can get the background map to re-project but can't seem to get the data to follow.
I'm not super versed in this yet, but here is what I have thus far:
Reading the data through NetCDF4:
from netCDF4 import Dataset
nc = Dataset('OR_ABI-L1b-RadF-
M3C13_G16_s20182831030383_e20182831041161_c20182831041217.nc')
data = nc.variables['Rad'][:]
Here I'm trying to get the geostationary info:
sat_h = nc.variables['goes_imager_projection'].perspective_point_height
X = nc.variables['x'][:] * sat_h
Y = nc.variables['y'][:] * sat_h
# Satellite longitude
sat_lon =
nc.variables['goes_imager_projection'].longitude_of_projection_origin
# Satellite sweep
sat_sweep = nc.variables['goes_imager_projection'].sweep_angle_axis
Here I'm taking projection data from the .nc file:
proj_var = nc.variables['goes_imager_projection']
sat_height = proj_var.perspective_point_height
central_lon = proj_var.longitude_of_projection_origin
semi_major = proj_var.semi_major_axis
semi_minor = proj_var.semi_minor_axis
print proj_var
<type 'netCDF4._netCDF4.Variable'>
int32 goes_imager_projection()
long_name: GOES-R ABI fixed grid projection
grid_mapping_name: geostationary
perspective_point_height: 35786023.0
semi_major_axis: 6378137.0
semi_minor_axis: 6356752.31414
inverse_flattening: 298.2572221
latitude_of_projection_origin: 0.0
longitude_of_projection_origin: -75.0
sweep_angle_axis: x
unlimited dimensions:
current shape = ()
filling on, default _FillValue of -2147483647 used
And here is a small snippet of my code that's relevant:
fig = plt.figure(figsize=(30,20))
globe = ccrs.Globe(semimajor_axis=semi_major, semiminor_axis=semi_minor)
proj = ccrs.Geostationary(central_longitude=central_lon,
satellite_height=sat_height, globe=globe)
ax = fig.add_subplot(1, 1, 1, projection=proj)
IR_img = ax.imshow(data[:,:],origin='upper',extent=(X.min(), X.max(), Y.min(), Y.max()),
cmap=IR_cmap,interpolation='nearest',vmin=162.,vmax=330.)
And an image of everyone playing nicely: Data and map working
When I try and get say a Plate Carree projection I try:
proj = ccrs.PlateCarree(central_longitude=central_lon,globe=globe)
And an image of my failure: Data and map not working
I've tried messing with the extent in the imshow method, I've tried adding a
transform=proj
in the imshow and no luck, it just gets hung up and I have to restart the kernel.
Clearly it is a lack of understanding on my part. If anyone can quickly and easily help/explain the way I want to change my projection from geostationary, I would greatly appreciate that.
I'm running archaic python2.
Thanks for looking.
EDIT: Problem seems to be resolved thanks to insight from DopplerShift and ajdawson, I guess I was maybe a little impatient/ignorant of how long a full disk transformation would take.
It looks like you need to specify the transform keyword to imshow. This keyword tells cartopy what coordinates your data are in, which in this case should be geostationary.
I don't have your dataset so I cannot test this, but the snippet below illustrates the concept. The projection and the transform are independent so you should define both. The value of the transform argument (crs in the example below) is fixed for the data set, but the projection can be anything you like (including the same as crs).
See this example of reprojecting a geostationary image: https://scitools.org.uk/cartopy/docs/v0.16/gallery/geostationary.html#sphx-glr-gallery-geostationary-py. Also see the guide to projection and transform arguments here: https://scitools.org.uk/cartopy/docs/v0.16/tutorials/understanding_transform.html.
globe = ccrs.Globe(semimajor_axis=semi_major, semiminor_axis=semi_minor)
crs = ccrs.Geostationary(central_longitude=central_lon,
satellite_height=sat_height, globe=globe)
proj = ccrs.PlateCarree(central_longitude=central_lon, globe=globe)
ax = fig.add_subplot(1, 1, 1, projection=proj)
IR_img = ax.imshow(data[:,:], origin='upper',
extent=(X.min(), X.max(), Y.min(), Y.max()),
transform=crs,
cmap=IR_cmap,
interpolation='nearest', vmin=162., vmax=330.)