0

I am having an issue displaying satellite data on the correctly projection using cartopy and imshow. I am unsure of the initial data projection, but I am trying to put it into Plate Carree with a limited area extent using bbox. There are nans in the array and I tried to use np.ma.masked_invalid to mask.

https://scitools.org.uk/cartopy/docs/v0.15/examples/geostationary.html

I tried the following examples, but it did not resolve the issue. [https://stackoverflow.com/questions/52883594/projecting-goes-16-geostationary-data-into-plate-carree-cartopy]2

What happens is the data is correct, but it is not plotting it into the correct projection with axes. I overlayed states and coastlines and they are all out of whack. This is using SCMI data for GOES-16 from NOAAPort.

Here is my code and image that shows the correct display (using different approach; last image on bottom) and an image showing the incorrect display (includes little box in top left corner and lighter colors) from this code... for comparison. It should be noted that contourf works fine. I tried to use pcolormesh, but it takes forever and using runs into a mem error.

edexServer='awips'
DataAccessLayer.changeEDEXHost(edexServer)
request = DataAccessLayer.newDataRequest()
request.setDatatype("satellite")
sector="ECONUS"
request.setLocationNames(sector)
request.setParameters("CH-13-10.35um")
cycles = DataAccessLayer.getAvailableTimes(request, True)
times = DataAccessLayer.getAvailableTimes(request)
fcstRun = DataAccessLayer.getForecastRun(cycles[-1], times)
response = DataAccessLayer.getGridData(request, [fcstRun[0]])
grid = response[0]
data = grid.getRawData()
imdata = np.ma.masked_invalid(np.atleast_2d(data))
#where_are_NaNs = np.isnan(data)
#data[where_are_NaNs] = 0
#data = ndimage.gaussian_filter(data, sigma=0.50, order=0, mode='nearest')
lons, lats = grid.getLatLonCoords()

bbox = [-130.0, -65.0, 20.0, 55.0]
fig = plt.figure(figsize=(6,6),dpi=200)
ax = plt.axes(projection=ccrs.PlateCarree())
#ax.set_extent(bbox, ccrs.PlateCarree())
coastlines = '/home/awips/.local/share/cartopy/shapefiles/natural_earth/physical/ne_10m_coastline.shp'
states = '/home/awips/.local/share/cartopy/shapefiles/natural_earth/physical/ne_10m_admin_1_states_provinces.shp'
globe = ccrs.Globe(semimajor_axis=6378137.0, semiminor_axis=6356752.31414)
crs = ccrs.Geostationary(central_longitude=-75, 
                         satellite_height=35786000, globe=globe)
cpt=loadCPT('/home/awips/Downloads/IR4AVHRR6.cpt')
cpt_convert = LinearSegmentedColormap('cpt', cpt)
ir_norm, ir_cmap = colortables.get_with_range('ir_bd', 150, 357)
#cs = ax.pcolormesh(imdata,cmap='coolwarm', transform=ccrs.PlateCarree())
#cs = ax.contourf(lons,lats,data,cmap='Greys_r', transform=ccrs.PlateCarree())
#levels = np.arange(170, 357, 5)
#image = ax.imshow(data, origin='upper', extent=(-149,-53,14,51), cmap=cpt_convert, vmin=150, vmax=375)
image = ax.imshow(imdata, origin='upper',cmap=cpt_convert, vmin=150, vmax= 375, transform = ccrs.PlateCarree())

enter image description here

enter image description here

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.