« Plotting GINI Water... | Main | Simple Plotting in... »

The Python AWIPS Data Access Framework can be used to query available grid parameters and levels if given a known Grid name (as of AWIPS 15.1.3 we can not query derived parameters, only parameters which have been directly decoded).

This example requests the U and V wind components for the GFS 40km CONUS and plots the wind speed (total vector) as a gridded contour (color-filled isotachs, essentially):

from awips.dataaccess import DataAccessLayer    # Select grid  DataAccessLayer.changeEDEXHost('edex-cloud.unidata.ucar.edu')  request = DataAccessLayer.newDataRequest()  request.setDatatype('grid')  request.setLocationNames('RAP13')  

grid is given as the data type, and RAP13 is given as the grid name. The function getAvailableParameters() returns a list of parameters:

# Print parm list  available_parms = DataAccessLayer.getAvailableParameters(request)  available_parms.sort()  for parm in available_parms:      print parm  

printed as

    AV      BLI      CAPE      CFRZR6hr      CICEP6hr      CIn      CP6hr      CRAIN6hr      CSNOW6hr      GH      P      P6hr      PMSL      PVV      PW      RH      SLI      T      TP6hr      VSS      WEASD      WGH      uW      vW  

List Available Levels for Parameter

Selecting the u wind component (uW) first and printing the return from getAvailableLevels():

request.setParameters('uW')  available_levels = DataAccessLayer.getAvailableLevels(request)  available_levels.sort()  for level in available_levels:      print level  

prints as

    1000.0MB      950.0MB      925.0MB      900.0MB      875.0MB      850.0MB      825.0MB      800.0MB      775.0MB      725.0MB      600.0MB      575.0MB      0.0_30.0BL      60.0_90.0BL      90.0_120.0BL      0.5PV      2.0PV      30.0_60.0BL      1.0PV      750.0MB      120.0_150.0BL      975.0MB      700.0MB      675.0MB      650.0MB      625.0MB      550.0MB      525.0MB      500.0MB      450.0MB      400.0MB      300.0MB      250.0MB      200.0MB      150.0MB      100.0MB      0.0TROP      1.5PV      150.0_180.0BL      350.0MB      10.0FHAG      0.0MAXW  

Construct Wind Field from U and V Components

Calling setLevels("10.0FHAG") and selecting the last available time t[-1]:

import numpy  from metpy.units import units  request.setLevels('10.0FHAG')  t = DataAccessLayer.getAvailableTimes(request)  # Select last time for u-wind  response = DataAccessLayer.getGridData(request, [t[-1]])  data_uw = response[-1]  lons,lats = data_uw.getLatLonCoords()  

Requesting the grid for v wind componeents (vW):

request.setParameters('vW')  response = DataAccessLayer.getGridData(request, [t[-1]])  data_vw = response[-1]    print 'Time :', t[-1]  print 'Model:', data_vw.getLocationName()  print 'Unit :', data_vw.getUnit()  print 'Parms :', data_uw.getParameter(), data_vw.getParameter()  print data_vw.getRawData().shape  

prints as

    Time : 2016-04-20 18:00:00 (240)      Model: GFS40      Unit : m*sec^-1      Parms : vW vW      (185, 129)      windArray = [[ 1.47078204  1.69705617  0.69296461 ...,  6.98621511 ...,  0.91923875  1.24450791   1.28693426]]   

Calculate absolute wind speed from U and V

The data objects data_uw.getRawData() and data_vw.getRawData()

spd = numpy.sqrt( data_uw.getRawData()**2 + data_vw.getRawData()**2 )  spd = spd * units.knot  print 'windArray =', spd  

Plotting a Grid with Basemap

Using matplotlib, numpy, and basemap:

python  %matplotlib inline  import matplotlib.tri as mtri  import matplotlib.pyplot as plt  from matplotlib.transforms import offset_copy  from mpl_toolkits.basemap import Basemap, cm  import numpy as np  from numpy import linspace, transpose  from numpy import meshgrid    plt.figure(figsize=(12, 12), dpi=100)    map = Basemap(projection='cyl',        resolution = 'c',        llcrnrlon = lons.min(), llcrnrlat = lats.min(),        urcrnrlon =lons.max(), urcrnrlat = lats.max()  )  map.drawcoastlines()  map.drawstates()  map.drawcountries()    #   # We have to reproject our grid, see https://stackoverflow.com/questions/31822553/m  #  x = linspace(0, map.urcrnrx, data_uw.getRawData().shape[1])  y = linspace(0, map.urcrnry, data_uw.getRawData().shape[0])  xx, yy = meshgrid(x, y)  ngrid = len(x)  rlons = np.repeat(np.linspace(np.min(lons), np.max(lons), ngrid),        ngrid).reshape(ngrid, ngrid)  rlats = np.repeat(np.linspace(np.min(lats), np.max(lats), ngrid),        ngrid).reshape(ngrid, ngrid).T  tli = mtri.LinearTriInterpolator(mtri.Triangulation(lons.flatten(),        lats.flatten()), spd.flatten())  rdata = tli(rlons, rlats)  #cs = map.contourf(rlons, rlats, rdata, latlon=True)  cs = map.contourf(rlons, rlats, rdata, latlon=True, vmin=0, vmax=20, cmap='BuPu')    # Add colorbar  cbar = map.colorbar(cs,location='bottom',pad='5%')    cbar.set_label('Wind Speed (Knots)')    # Show plot  plt.show()  

or use pcolormesh rather than contourf

python  plt.figure(figsize=(12, 12), dpi=100)  map = Basemap(projection='cyl',        resolution = 'c',        llcrnrlon = lons.min(), llcrnrlat = lats.min(),        urcrnrlon =lons.max(), urcrnrlat = lats.max()  )  map.drawcoastlines()  map.drawstates()  map.drawcountries()  cs = map.pcolormesh(rlons, rlats, rdata, latlon=True, vmin=0, vmax=20, cmap='BuPu')  

Plotting a Grid with Cartopy

python  import os  import matplotlib.pyplot as plt  import numpy as np  import iris  import cartopy.crs as ccrs  from cartopy import config    lon,lat = data.getLatLonCoords()  plt.figure(figsize=(12, 12), dpi=100)  ax = plt.axes(projection=ccrs.PlateCarree())  cs = plt.contourf(rlons, rlats, rdata, 60, transform=ccrs.PlateCarree(), vmin=0, vmax=20, cmap='BuPu')  ax.coastlines()  ax.gridlines()    # add colorbar  cbar = plt.colorbar(orientation='horizontal')  cbar.set_label('Wind Speed (Knots)')  plt.show()