[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[AWIPS #LWN-904229]: Displaying Full Disk Through Python DAF

Hi Justin,

I just wanted to copy some of our correspondence into our support system so 
it's available to reference
later on.  Glad you got things working, and it's good to know that imshow can 
be a feasible alternative
for plotting larger imagery!


Justin:

I am running into the issue where EFD and ECONUS plot just fine, but plotting 
WFD and WCONUS plot blank images.
I assume it has to do with the lat/lon stuff, but am a bit at a loss. I can 
plot the exact same coordinates 
with ECONUS and EFD, but switch the sector and nada.

Here is my definition to plot cartography. I am plotting it using local 
shapefiles.

# This generates the map/cartography. It is using local files for now.
def make_map(bbox, projection=ccrs.PlateCarree()):    
    coastlines = 
'/home/awips/.local/share/cartopy/shapefiles/natural_earth/physical/ne_50m_coastline.shp'
    states = 
'/home/awips/.local/share/cartopy/shapefiles/natural_earth/physical/ne_110m_admin_1_states_provinces.shp'
    counties = 
'/home/awips/.local/share/cartopy/shapefiles/natural_earth/physical/cb_2018_us_county_5m.shp'
    lakes = 
'/home/awips/.local/share/cartopy/shapefiles/natural_earth/physical/ne_10m_lakes.shp'
    fig, ax = plt.subplots(figsize=(6.4,3.6),dpi=300,
            subplot_kw=dict(projection=projection))
    ax.set_extent(bbox, crs=ccrs.PlateCarree())
    ax.set_title(title,fontname='FreeSans',size=6, loc="left",y=0.980) # Title
    
ax.add_geometries(Reader(coastlines).geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='black',linewidth=0.25)
    # This adds the states from a shapefile.
    
#ax.add_geometries(Reader(states).geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='gray',linewidth=0.50)
    
ax.add_geometries(Reader(counties).geometries(),ccrs.PlateCarree(),facecolor='none',linewidth=0.25,edgecolor='black')
    
ax.add_geometries(Reader(lakes).geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='black',linewidth=0.25)

    gl = ax.gridlines(draw_labels=True, linewidth=0.15, 
linestyle='--',color='gray')
    gl.xlabels_top = gl.ylabels_right = False
    gl.xlabel_style = {'size': 4, 'color': 'gray'}
    gl.ylabel_style = {'size': 4, 'color': 'gray'}
    gl.xformatter = LONGITUDE_FORMATTER
    gl.yformatter = LATITUDE_FORMATTER
    return fig, ax

On the actual code to retrieve/plot data, here is what I have for imshow.
        edexServer='myserver'
        DataAccessLayer.changeEDEXHost(edexServer)
        request = DataAccessLayer.newDataRequest()
        request.setDatatype("satellite")
        request.setLocationNames(sec)
        request.setParameters(channel)
        #cycles = DataAccessLayer.getAvailableTimes(request, True)
        times = DataAccessLayer.getAvailableTimes(request)
        response = DataAccessLayer.getGridData(request, [times[Nframe]])
        grid = response[0]
        imdata = grid.getRawData()
        if (channel == 1 or channel == 2 or channel == 3):
            data = imdata
        else:
            data = imdata - 273.15

# Set up the projection using known parameters
        globe = ccrs.Globe(semimajor_axis=6378137.0, 
semiminor_axis=6356752.31414)
        crs = ccrs.Geostationary(central_longitude=-75,
                                 satellite_height=35786023, globe=globe, 
sweep_axis='x')

        # Convert grid lons/lats to native projected coordinates
        lons, lats = grid.getLatLonCoords()
        x, y, _ = crs.transform_points(ccrs.PlateCarree(), lons, lats).T

        bbox = [newBounds[0],newBounds[1],newBounds[2],newBounds[3]]

        fig, ax = make_map(bbox=bbox)
        # Converts the CPT file to be used in Python
        cpt = loadCPT(cptI)
        cpt_convert = LinearSegmentedColormap('cpt', cpt)
        check = str(times[Nframe])
        if check.find(".") == -1:
            validTime = datetime.datetime.strptime(str(times[Nframe]),'%Y-%m-%d 
%H:%M:%S')
        else:
            validTime = datetime.datetime.strptime(str(times[Nframe]),'%Y-%m-%d 
%H:%M:%S.%f')

        image = ax.imshow(data, origin='upper', extent=(np.nanmin(x), 
np.nanmax(x), np.nanmin(y),
                          
np.nanmax(y)),cmap=cpt_convert,vmin=vmn,vmax=vmx,transform=crs)
        cbar = fig.colorbar(image, shrink=0.7, orientation='vertical',pad=0.02)
        cbar.ax.tick_params(labelsize=4,right='off', 
labelleft='off',labelright='on',pad=0.03)


Shay:

One thing that might help: the West images cross the international dateline, so 
I think your assumption 
about the problem lying with lat/lons might be correct.  If you look at the 
notebook I sent you yesterday, 
in the last section dealing with WCONUS and WFD I had two lines of code that 
might be useful for your example:

    # Shift the lons, because we're going to shift the projection we're 
plotting into
    lons = lons+180
    ....
    fig3, ax3 = make_map(bbox=bbox, 
projection=ccrs.PlateCarree(central_longitude=180))

When I call the make_map method for WCONUS and WFD, I'm actually passing in a 
different projection, so the 
default one in the method/function is overridden.  This projection moves the 
center of the coordinate system
so that we're not crossing the dateline from 179.999 --> -179.999.
I also shifted all of the lons by 180 degrees to account for this shift in 
projection as well.
I **think** this was the right thing to do, and **may** be all you need to get 
your code to work?
Let me know, and if it doesn't work, please send me your local shape files and 
I'll try and get your code up
and running and working.
Thanks!


Justin:

Thanks Shay --  that actually helped me figure things out a bit further. I was 
using the center for GOES-East; 
where I needed to change it to GOES-West.

Before:
globe = ccrs.Globe(semimajor_axis=6378137.0, semiminor_axis=6356752.31414)
crs = ccrs.Geostationary(central_longitude=-75.0,
                         satellite_height=35786023, globe=globe, sweep_axis='x')

After:
globe = ccrs.Globe(semimajor_axis=6378137.0, semiminor_axis=6356752.31414)
crs = ccrs.Geostationary(central_longitude=-137.2,
                         satellite_height=35786023, globe=globe, sweep_axis='x')

Things are now plotting correctly.

--Shay Carter

AWIPS Software Engineer
UCAR - Unidata

If you're interested, please feel free to fill out a survey about the support 
you receive: 
https://docs.google.com/forms/d/e/1FAIpQLSeDIkdk8qUMgq8ZdM4jhP-ubJPUOr-mJMQgxInwoAWoV5QcOw/viewform

Ticket Details
===================
Ticket ID: LWN-904229
Department: Support AWIPS
Priority: Normal
Status: Open
===================
NOTE: All email exchanges with Unidata User Support are recorded in the Unidata 
inquiry tracking system and then made publicly available through the web.  If 
you do not want to have your interactions made available in this way, you must 
let us know in each email you send to us.