underscore needed in file names

This commit is contained in:
mjames-upc 2016-04-20 19:12:34 -05:00
parent dbcb2083f2
commit df8ee6e885
23 changed files with 1211 additions and 1 deletions

View file

@ -0,0 +1,3 @@
=========
Dev Guide
=========

View file

@ -11,7 +11,7 @@ List Available Parameters for a Grid Name
from awips.dataaccess import DataAccessLayer
# Select HRRR
DataAccessLayer.changeEDEXHost("edex-cloud.unidata.ucar.edu")
DataAccessLayer.changeEDEXHost("edex.unidata.ucar.edu")
request = DataAccessLayer.newDataRequest()
request.setDatatype("grid")
request.setLocationNames("GFS40")

View file

@ -0,0 +1,335 @@
============
Gridded Data
============
`Notebook <http://nbviewer.ipython.org/github/Unidata/python-awips/blob/master/examples/notebooks/Gridded_Data.ipynb>`_
EDEX Grid Inventory
-------------------
.. code:: python
from awips.dataaccess import DataAccessLayer
# Set host
DataAccessLayer.changeEDEXHost("edex-cloud.unidata.ucar.edu")
# Init data request
request = DataAccessLayer.newDataRequest()
# Set datatype
request.setDatatype("grid")
# Get a list of all available models
available_grids = DataAccessLayer.getAvailableLocationNames(request)
# Sort
available_grids.sort()
for grid in available_grids:
print grid
.. parsed-literal::
AUTOSPE
AVN211
AVN225
DGEX
ECMF-Global
ECMF1
ECMF10
ECMF11
ECMF12
ECMF2
ECMF3
ECMF4
ECMF5
ECMF6
ECMF7
ECMF8
ECMF9
ESTOFS
ETA
FFG-ALR
FFG-FWR
FFG-KRF
FFG-MSR
FFG-ORN
FFG-RHA
FFG-RSA
FFG-TAR
FFG-TIR
FFG-TUA
GFS
GFS40
GFSGuide
GFSLAMP5
GribModel:9:151:172
HFR-EAST_6KM
HFR-EAST_PR_6KM
HFR-US_EAST_DELAWARE_1KM
HFR-US_EAST_FLORIDA_2KM
HFR-US_EAST_NORTH_2KM
HFR-US_EAST_SOUTH_2KM
HFR-US_EAST_VIRGINIA_1KM
HFR-US_HAWAII_1KM
HFR-US_HAWAII_2KM
HFR-US_HAWAII_6KM
HFR-US_WEST_500M
HFR-US_WEST_CENCAL_2KM
HFR-US_WEST_LOSANGELES_1KM
HFR-US_WEST_LOSOSOS_1KM
HFR-US_WEST_NORTH_2KM
HFR-US_WEST_SANFRAN_1KM
HFR-US_WEST_SOCAL_2KM
HFR-US_WEST_WASHINGTON_1KM
HFR-WEST_6KM
HPCGuide
HPCqpf
HPCqpfNDFD
HRRR
LAMP2p5
MPE-Local-ORN
MPE-Local-RHA
MPE-Local-RSA
MPE-Local-TAR
MPE-Local-TIR
MPE-Mosaic-MSR
MPE-Mosaic-ORN
MPE-Mosaic-RHA
MPE-Mosaic-TAR
MPE-Mosaic-TIR
MRMS_1000
NAM12
NAM40
NCWF
NOHRSC-SNOW
NamDNG
NamDNG5
QPE-ALR
QPE-Auto-TUA
QPE-FWR
QPE-KRF
QPE-MSR
QPE-RFC-RSA
QPE-RFC-STR
QPE-TIR
QPE-TUA
QPE-XNAV-ALR
QPE-XNAV-FWR
QPE-XNAV-KRF
QPE-XNAV-MSR
QPE-XNAV-RHA
QPE-XNAV-SJU
QPE-XNAV-TAR
QPE-XNAV-TIR
QPE-XNAV-TUA
RAP13
RAP40
RCM
RFCqpf
RTMA
RTMA5
UKMET-Global
UKMET37
UKMET38
UKMET39
UKMET40
UKMET41
UKMET42
UKMET43
UKMET44
URMA25
estofsPR
fnmocWave
**LocationNames** is different for different plugins - radar is icao -
satellite is sector
Requesting a Grid
-----------------
.. code:: python
# Grid request
request.setLocationNames('RAP40')
request.setParameters("RH")
request.setLevels("850MB")
# Get available times
t = DataAccessLayer.getAvailableTimes(request)
# Select last available time [-1]
response = DataAccessLayer.getGridData(request, [t[0]])
data = response[0]
lon,lat = data.getLatLonCoords()
# Print info
print 'Time :', t[-1]
print 'Model:', data.getLocationName()
print 'Unit :', data.getUnit()
print 'Parm :', data.getParameter()
# Print data array
print data.getRawData().shape
print data.getRawData()
print "lat array =", lat
print "lon array =", lon
.. parsed-literal::
Time : 2016-02-23 15:00:00 (12)
Model: RAP40
Unit : %
Parm : RH
(151, 113)
[[ 93.05456543 93.05456543 87.05456543 ..., 73.05456543 72.05456543
71.05456543]
[ 70.05456543 70.05456543 67.05456543 ..., 69.05456543 46.05456924
37.05456924]
[ 40.05456924 56.05456924 68.05456543 ..., 51.05456924 73.05456543
74.05456543]
...,
[ 65.05456543 62.05456924 63.05456924 ..., 67.05456543 65.05456543
46.05456924]
[ 48.05456924 59.05456924 62.05456924 ..., 4.05456877 5.05456877
5.05456877]
[ 7.05456877 8.05456829 10.05456829 ..., 91.05456543 95.05456543
95.05456543]]
lat array = [[ 54.24940109 54.35071945 54.45080566 ..., 57.9545517 57.91926193
57.88272858]
[ 57.84495163 57.80593109 57.76566696 ..., 58.07667542 58.08861542
58.09931183]
[ 58.10876846 58.11697769 58.12394714 ..., 56.40270996 56.46187973
56.51980972]
...,
[ 19.93209648 19.89832115 19.86351395 ..., 20.054636 20.06362152
20.07156372]
[ 20.0784626 20.08431816 20.08912849 ..., 18.58354759 18.63155174
18.67854691]
[ 18.72453308 18.76950836 18.81346893 ..., 17.49624634 17.42861557
17.36001205]]
lon array = [[-139.83120728 -139.32348633 -138.81448364 ..., -79.26060486
-78.70166016 -78.14326477]
[ -77.58544922 -77.02822876 -76.47161865 ..., -100.70157623
-100.13801575 -99.57427216]
[ -99.01037598 -98.44634247 -97.88218689 ..., -121.69165039
-121.15060425 -120.60871887]
...,
[ -82.65139008 -82.26644897 -81.88170624 ..., -98.52494049
-98.13802338 -97.75105286]
[ -97.36403656 -96.97698212 -96.58989716 ..., -113.07767487
-112.69831085 -112.31866455]
[-111.93874359 -111.5585556 -111.17810822 ..., -69.85433197
-69.48160553 -69.10926819]]
Plotting a Grid with Basemap
----------------------------
Using **matplotlib**, **numpy**, and **basemap**:
.. code:: python
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)
lons,lats = data.getLatLonCoords()
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.getRawData().shape[1])
y = linspace(0, map.urcrnry, data.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()), data.getRawData().flatten())
rdata = tli(rlons, rlats)
cs = map.contourf(rlons, rlats, rdata, latlon=True, vmin=0, vmax=100, cmap='YlGn')
# add colorbar.
cbar = map.colorbar(cs,location='bottom',pad="5%")
cbar.set_label(data.getParameter() + data.getUnit() )
# Show plot
plt.show()
.. image:: Gridded_Data_files/Gridded_Data_5_0.png
or use **pcolormesh** rather than **contourf**
.. code:: 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=100, cmap='YlGn')
.. image:: Gridded_Data_files/Gridded_Data_7_0.png
Plotting a Grid with Cartopy
----------------------------
.. code:: 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=100, cmap='YlGn')
ax.coastlines()
ax.gridlines()
# add colorbar
cbar = plt.colorbar(orientation='horizontal')
cbar.set_label(data.getParameter() + data.getUnit() )
plt.show()
.. image:: Gridded_Data_files/Gridded_Data_9_0.png

Binary file not shown.

After

Width:  |  Height:  |  Size: 446 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 395 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 634 KiB

View file

@ -0,0 +1,112 @@
===================================
NEXRAD Level 3 Plot with Matplotlib
===================================
`Notebook <http://nbviewer.ipython.org/github/Unidata/python-awips/blob/master/examples/notebooks/NEXRAD_Level_3_Plot_with_Matplotlib.ipynb>`_
NEXRAD Level 3 Plot with Matplotlib
===================================
.. code:: python
%matplotlib inline
from awips.dataaccess import DataAccessLayer
from awips import ThriftClient, RadarCommon
from dynamicserialize.dstypes.com.raytheon.uf.common.time import TimeRange
from dynamicserialize.dstypes.com.raytheon.uf.common.dataplugin.radar.request import GetRadarDataRecordRequest
from datetime import datetime
from datetime import timedelta
import matplotlib.pyplot as plt
import numpy as np
from numpy import ma
# use metpy for color table
from metpy.plots import ctables
# Set EDEX server and radar site
DataAccessLayer.changeEDEXHost("edex.unidata.ucar.edu")
request = DataAccessLayer.newDataRequest()
request.setDatatype("radar")
request.setLocationNames("klzk")
datatimes = DataAccessLayer.getAvailableTimes(request)
# Get last available time
timerange = datatimes[-1].validPeriod
dateTimeStr = str(datatimes[-1])
# Buffer length in seconds
buffer = 60
dateTime = datetime.strptime(dateTimeStr, "%Y-%m-%d %H:%M:%S")
beginRange = dateTime - timedelta(0, buffer)
endRange = dateTime + timedelta(0, buffer)
timerange = TimeRange(beginRange, endRange)
print "using time",dateTimeStr
print "buffer by",buffer
print "using range",timerange
client = ThriftClient.ThriftClient(edex)
request = GetRadarDataRecordRequest()
request.setRadarId(site)
request.setPrimaryElevationAngle("0.5")
request.setTimeRange(timerange)
fig, axes = plt.subplots(1, 2, figsize=(15, 8))
for v, ctable, ax in zip((94, 99), ('NWSReflectivity', 'NWSVelocity'), axes):
request.setProductCode(v)
response = client.sendRequest(request)
if response.getData():
for record in response.getData():
idra = record.getHdf5Data()
rdat,azdat,depVals,threshVals = RadarCommon.get_hdf5_data(idra)
dim = rdat.getDimension()
yLen,xLen = rdat.getSizes()
array = rdat.getByteData()
# get data for azimuth angles if we have them.
if azdat :
azVals = azdat.getFloatData()
az = np.array(RadarCommon.encode_radial(azVals))
dattyp = RadarCommon.get_data_type(azdat)
az = np.append(az,az[-1])
print "found",v,record.getDataTime()
header = RadarCommon.get_header(record, format, xLen, yLen, azdat, "description")
rng = np.linspace(0, xLen, xLen + 1)
xlocs = rng * np.sin(np.deg2rad(az[:, np.newaxis]))
ylocs = rng * np.cos(np.deg2rad(az[:, np.newaxis]))
multiArray = np.reshape(array, (-1, xLen))
data = ma.array(multiArray)
data[data==0] = ma.masked
# Plot the data
norm, cmap = ctables.registry.get_with_steps(ctable, 16, 16)
ax.pcolormesh(xlocs, ylocs, data, norm=norm, cmap=cmap)
ax.set_aspect('equal', 'datalim')
multp = 100*(2*xLen/460)
ax.set_xlim(-multp,multp)
ax.set_ylim(-multp,multp)
# This is setting x/ylim on gate/pixel and not km
plt.show()
.. parsed-literal::
using time 2016-04-11 23:02:22
buffer by 60
using range (Apr 11 16 23:01:22 , Apr 11 16 23:03:22 )
found 94 2016-04-11 23:02:22
found 99 2016-04-11 23:02:22
.. image:: NEXRAD_Level_3_Plot_with_Matplotlib_files/NEXRAD_Level_3_Plot_with_Matplotlib_1_1.png

Binary file not shown.

After

Width:  |  Height:  |  Size: 47 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 55 KiB

View file

@ -0,0 +1,134 @@
==============================
Plotting a Sounding with MetPy
==============================
`Notebook <http://nbviewer.ipython.org/github/Unidata/python-awips/blob/master/examples/notebooks/Plotting_a_Sounding_with_MetPy.ipynb>`_
Plotting a Sounding with MetPy
==============================
.. code:: python
%matplotlib inline
from awips.dataaccess import DataAccessLayer
import matplotlib.tri as mtri
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
import numpy as np
from metpy.calc import get_wind_components, lcl, dry_lapse, parcel_profile
from metpy.plots import SkewT, Hodograph
from metpy.units import units, concatenate
plt.rcParams['figure.figsize'] = (12, 14)
# Set EDEX host
DataAccessLayer.changeEDEXHost("edex.unidata.ucar.edu")
request = DataAccessLayer.newDataRequest()
# Data type bufrua
request.setDatatype("bufrua")
# Parameters
request.setParameters("tpMan","tdMan","prMan","htMan","wdMan","wsMan")
# Station ID (name doesn't work yet)
request.setLocationNames("72469")
datatimes = DataAccessLayer.getAvailableTimes(request)
# Get most recent record
response = DataAccessLayer.getGeometryData(request,times=datatimes[-1].validPeriod)
# Initialize data arrays
tpMan,tdMan,prMan,htMan,wdMan,wsMan = [],[],[],[],[],[]
# Build ordered arrays
for ob in response:
print float(ob.getString("prMan")), float(ob.getString("wsMan"))
tpMan.append(float(ob.getString("tpMan")))
tdMan.append(float(ob.getString("tdMan")))
prMan.append(float(ob.getString("prMan")))
htMan.append(float(ob.getString("htMan")))
wdMan.append(float(ob.getString("wdMan")))
wsMan.append(float(ob.getString("wsMan")))
# we can use units.* here...
T = np.array(tpMan)-273.15
Td = np.array(tdMan)-273.15
p = np.array(prMan)/100
height = np.array(htMan)
direc = np.array(wdMan)
spd = np.array(wsMan)
u, v = get_wind_components(spd, np.deg2rad(direc))
p = p * units.mbar
T = T * units.degC
Td = Td * units.degC
spd = spd * units.knot
direc = direc * units.deg
# Create a skewT plot
skew = SkewT()
# Plot the data using normal plotting functions, in this case using
# log scaling in Y, as dictated by the typical meteorological plot
skew.plot(p, T, 'r')
skew.plot(p, Td, 'g')
skew.plot_barbs(p, u, v)
skew.ax.set_ylim(1000, 100)
skew.ax.set_xlim(-40, 60)
# Calculate LCL height and plot as black dot
l = lcl(p[0], T[0], Td[0])
lcl_temp = dry_lapse(concatenate((p[0], l)), T[0])[-1].to('degC')
skew.plot(l, lcl_temp, 'ko', markerfacecolor='black')
# Calculate full parcel profile and add to plot as black line
prof = parcel_profile(p, T[0], Td[0]).to('degC')
skew.plot(p, prof, 'k', linewidth=2)
# Example of coloring area between profiles
skew.ax.fill_betweenx(p, T, prof, where=T>=prof, facecolor='blue', alpha=0.4)
skew.ax.fill_betweenx(p, T, prof, where=T<prof, facecolor='red', alpha=0.4)
# An example of a slanted line at constant T -- in this case the 0 isotherm
l = skew.ax.axvline(0, color='c', linestyle='--', linewidth=2)
# Draw hodograph
ax_hod = inset_axes(skew.ax, '40%', '40%', loc=3)
h = Hodograph(ax_hod, component_range=80.)
h.add_grid(increment=20)
h.plot_colormapped(u, v, spd)
# Show the plot
plt.show()
.. parsed-literal::
83900.0 1.5
100000.0 -9999998.0
92500.0 -9999998.0
85000.0 -9999998.0
70000.0 0.5
50000.0 6.09999990463
40000.0 3.0
30000.0 7.69999980927
25000.0 16.8999996185
20000.0 7.19999980927
15000.0 10.1999998093
10000.0 13.8000001907
7000.0 9.19999980927
5000.0 7.69999980927
3000.0 5.59999990463
2000.0 6.59999990463
1000.0 10.8000001907
700.0 5.09999990463
500.0 -9999.0
300.0 -9999.0
200.0 -9999.0
100.0 -9999.0
.. image:: Plotting_a_Sounding_with_MetPy_files/Plotting_a_Sounding_with_MetPy_1_1.png

Binary file not shown.

After

Width:  |  Height:  |  Size: 104 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 110 KiB

View file

@ -0,0 +1,213 @@
===========
Surface Obs
===========
`Notebook <http://nbviewer.ipython.org/github/Unidata/python-awips/blob/master/examples/notebooks/Surface_Obs.ipynb>`_
.. code:: python
from awips.dataaccess import DataAccessLayer
# Set host
DataAccessLayer.changeEDEXHost("edex.unidata.ucar.edu")
# Init data request
request = DataAccessLayer.newDataRequest()
request.setDatatype("obs")
request.setLocationNames("KBJC")
datatimes = DataAccessLayer.getAvailableTimes(request)
time = datatimes[-1].validPeriod
# "presWeather","skyCover","skyLayerBase"
# are multi-dimensional... deal with these later
request.setParameters(
"stationName",
"timeObs",
"wmoId",
"autoStationType",
"elevation",
"reportType",
"temperature",
"tempFromTenths",
"dewpoint",
"dpFromTenths",
"windDir",
"windSpeed",
"windGust",
"visibility",
"altimeter",
"seaLevelPress",
"pressChange3Hour",
"pressChangeChar",
"maxTemp24Hour",
"minTemp24Hour",
"precip1Hour",
"precip3Hour",
"precip6Hour",
"precip24Hour"
)
response = DataAccessLayer.getGeometryData(request,times=time)
for ob in response:
print "getParameters is",ob.getParameters()
print len(ob.getParameters())
#getParameters
print ob.getString("stationName"), "from", ob.getDataTime().getRefTime()
print "stationName is",ob.getString("stationName")
print "timeObs is",ob.getString("timeObs")
print "wmoId is",ob.getString("wmoId")
print "autoStationType is",ob.getString("autoStationType")
print "elevation is",ob.getString("elevation")
print "reportType is",ob.getString("reportType")
print "temperature is",ob.getString("temperature")
print "tempFromTenths is",ob.getString("tempFromTenths")
print "dewpoint is",ob.getString("dewpoint")
print "dpFromTenths is",ob.getString("dpFromTenths")
print "windDir is",ob.getString("windDir")
print "windSpeed is",ob.getString("windSpeed")
print "windGust is",ob.getString("windGust")
print "visibility is",ob.getString("visibility")
print "altimeter is",ob.getString("altimeter")
print "seaLevelPress is",ob.getString("seaLevelPress")
print "pressChange3Hour is",ob.getString("pressChange3Hour")
print "pressChangeChar is",ob.getString("pressChangeChar")
print "maxTemp24Hour is",ob.getString("maxTemp24Hour")
print "minTemp24Hour is",ob.getString("minTemp24Hour")
print "precip1Hour is",ob.getString("precip1Hour")
print "precip3Hour is",ob.getString("precip3Hour")
print "precip6Hour is",ob.getString("precip6Hour")
print "precip24Hour is",ob.getString("precip24Hour")
.. code:: python
# multi-dimensional present WX
request = DataAccessLayer.newDataRequest()
request.setDatatype("obs")
request.setLocationNames("KBJC")
request.setParameters("presWeather")
response = DataAccessLayer.getGeometryData(request,times=time)
for ob in response:
print "getParameters is",ob.getParameters()
print ob.getString("presWeather")
# multi-dimensional Sky Condition
request.setParameters("skyCover", "skyLayerBase")
response = DataAccessLayer.getGeometryData(request,times=time)
for ob in response:
print ob.getString("skyCover")
print ob.getString("skyLayerBase")
Synop/Marine
------------
.. code:: python
from awips.dataaccess import DataAccessLayer
DataAccessLayer.changeEDEXHost("edex.unidata.ucar.edu")
request = DataAccessLayer.newDataRequest()
request.setDatatype("sfcobs")
request.setLocationNames("72421") # Covington, Kentucky (KCVG)
request.setParameters("stationId","timeObs","elevation","reportType",
"wx_present","visibility","seaLevelPress","stationPress",
"pressChange3Hour","pressChangeChar","temperature",
"dewpoint","seaSurfaceTemp","wetBulb","windDir",
"windSpeed","equivWindSpeed10m","windGust","precip1Hour",
"precip6Hour","precip24Hour" )
datatimes = DataAccessLayer.getAvailableTimes(request)
time = datatimes[-1].validPeriod
response = DataAccessLayer.getGeometryData(request,times=time)
print response
for ob in response:
print "getParameters is",ob.getParameters()
print len(ob.getParameters())
Profiler
--------
.. code:: python
MULTI_DIM_PARAMS = set(['vComponent', 'uComponent', 'peakPower',
'levelMode', 'uvQualityCode', 'consensusNum',
'HorizSpStdDev', 'wComponent', 'height',
'VertSpStdDev'])
request = DataAccessLayer.newDataRequest("profiler")
request.setParameters('numProfLvls', 'elevation', 'windDirSfc', 'validTime',
'windSpeedSfc', 'pressure', 'submode', 'relHumidity',
'profilerId', 'rainRate', 'temperature')
request.getParameters().extend(MULTI_DIM_PARAMS)
datatimes = DataAccessLayer.getAvailableTimes(request)
time = datatimes[-1].validPeriod
response = DataAccessLayer.getGeometryData(request,times=time)
print response
for ob in response:
print "getParameters is",ob.getParameters()
print len(ob.getParameters())
ACARS
-----
.. code:: python
request = DataAccessLayer.newDataRequest("acars")
request.setParameters("tailNumber", "receiver", "pressure", "flightPhase",
"rollAngleQuality", "temp", "windDirection", "windSpeed",
"humidity", "mixingRatio", "icing")
datatimes = DataAccessLayer.getAvailableTimes(request)
time = datatimes[-1].validPeriod
response = DataAccessLayer.getGeometryData(request,times=time)
print response
for ob in response:
print "getParameters is",ob.getParameters()
print len(ob.getParameters())
AIREP
-----
.. code:: python
request = DataAccessLayer.newDataRequest("airep")
request.setParameters("id", "flightLevel", "temp", "windDirection", "windSpeed",
"flightWeather", "flightHazard", "flightConditions")
datatimes = DataAccessLayer.getAvailableTimes(request)
time = datatimes[-1].validPeriod
response = DataAccessLayer.getGeometryData(request,times=time)
print response
for ob in response:
print "getParameters is",ob.getParameters()
print len(ob.getParameters())
PIREP
-----
.. code:: python
MULTI_DIM_PARAMS = set(["hazardType",
"turbType", "turbBaseHeight", "turbTopHeight",
"iceType", "iceBaseHeight", "iceTopHeight",
"skyCover1", "skyCover2", "skyBaseHeight", "skyTopHeight"
])
request = DataAccessLayer.newDataRequest("pirep")
request.setParameters('id', 'flightLevel', 'temp', 'windDirection', 'windSpeed',
'horzVisibility', 'aircraftType', 'weatherGroup')
request.getParameters().extend(MULTI_DIM_PARAMS)
datatimes = DataAccessLayer.getAvailableTimes(request)
time = datatimes[-1].validPeriod
response = DataAccessLayer.getGeometryData(request,times=time)
print response
for ob in response:
print "getParameters is",ob.getParameters()
print len(ob.getParameters())

View file

@ -0,0 +1,144 @@
===========================
Surface Obs Plot with MetPy
===========================
`Notebook <http://nbviewer.ipython.org/github/Unidata/python-awips/blob/master/examples/notebooks/Surface_Obs_Plot_with_MetPy.ipynb>`_
Based on the MetPy example `"Station Plot with
Layout" <http://metpy.readthedocs.org/en/latest/examples/generated/Station_Plot_with_Layout.html>`_
.. code:: python
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from awips.dataaccess import DataAccessLayer
from metpy.calc import get_wind_components
from metpy.cbook import get_test_data
from metpy.plots import StationPlot, StationPlotLayout, simple_layout
from metpy.units import units
# Initialize
DataAccessLayer.changeEDEXHost("edex.unidata.ucar.edu")
data,latitude,longitude,stationName,temperature,dewpoint,seaLevelPress,windDir,windSpeed = [],[],[],[],[],[],[],[],[]
request = DataAccessLayer.newDataRequest()
request.setDatatype("obs")
#
# we need to set one station to query latest time. this is hack-y and should be fixed
# because when you DON'T set a location name, you tend to get a single observation
# that came in a second ago, so your "latest data for the last time for all stations"
# data array consists of one village in Peru and time-matching is suspect right now.
#
# So here take a known US station (OKC) and hope/assume that a lot of other stations
# are also reporting (and that this is a 00/20/40 ob).
#
request.setLocationNames("KOKC")
datatimes = DataAccessLayer.getAvailableTimes(request)
# Get most recent time for location
time = datatimes[0].validPeriod
# "presWeather","skyCover","skyLayerBase"
# are multi-dimensional(??) and returned seperately (not sure why yet)... deal with those later
request.setParameters("presWeather","skyCover", "skyLayerBase","stationName","temperature","dewpoint","windDir","windSpeed",
"seaLevelPress","longitude","latitude")
request.setLocationNames()
response = DataAccessLayer.getGeometryData(request,times=time)
print time
PRES_PARAMS = set(["presWeather"])
SKY_PARAMS = set(["skyCover", "skyLayerBase"])
# Build ordered arrays
wx,cvr,bas=[],[],[]
for ob in response:
#print ob.getParameters()
if set(ob.getParameters()) & PRES_PARAMS :
wx.append(ob.getString("presWeather"))
continue
if set(ob.getParameters()) & SKY_PARAMS :
cvr.append(ob.getString("skyCover"))
bas.append(ob.getNumber("skyLayerBase"))
continue
latitude.append(float(ob.getString("latitude")))
longitude.append(float(ob.getString("longitude")))
#stationName.append(ob.getString("stationName"))
temperature.append(float(ob.getString("temperature")))
dewpoint.append(float(ob.getString("dewpoint")))
seaLevelPress.append(float(ob.getString("seaLevelPress")))
windDir.append(float(ob.getString("windDir")))
windSpeed.append(float(ob.getString("windSpeed")))
print len(wx)
print len(temperature)
# Convert
data = dict()
data['latitude'] = np.array(latitude)
data['longitude'] = np.array(longitude)
data['air_temperature'] = np.array(temperature)* units.degC
data['dew_point_temperature'] = np.array(dewpoint)* units.degC
#data['air_pressure_at_sea_level'] = np.array(seaLevelPress)* units('mbar')
u, v = get_wind_components(np.array(windSpeed) * units('knots'),
np.array(windDir) * units.degree)
data['eastward_wind'], data['northward_wind'] = u, v
# Convert the fraction value into a code of 0-8, which can be used to pull out
# the appropriate symbol
#data['cloud_coverage'] = (8 * data_arr['cloud_fraction']).astype(int)
# Map weather strings to WMO codes, which we can use to convert to symbols
# Only use the first symbol if there are multiple
#wx_text = make_string_list(data_arr['weather'])
#wx_codes = {'':0, 'HZ':5, 'BR':10, '-DZ':51, 'DZ':53, '+DZ':55,
# '-RA':61, 'RA':63, '+RA':65, '-SN':71, 'SN':73, '+SN':75}
#data['present_weather'] = [wx_codes[s.split()[0] if ' ' in s else s] for s in wx]
# Set up the map projection
import cartopy.crs as ccrs
import cartopy.feature as feat
from matplotlib import rcParams
rcParams['savefig.dpi'] = 255
proj = ccrs.LambertConformal(central_longitude=-95, central_latitude=35,
standard_parallels=[35])
state_boundaries = feat.NaturalEarthFeature(category='cultural',
name='admin_1_states_provinces_lines',
scale='110m', facecolor='none')
# Create the figure
fig = plt.figure(figsize=(20, 10))
ax = fig.add_subplot(1, 1, 1, projection=proj)
# Add map elements
ax.add_feature(feat.LAND, zorder=-1)
ax.add_feature(feat.OCEAN, zorder=-1)
ax.add_feature(feat.LAKES, zorder=-1)
ax.coastlines(resolution='110m', zorder=2, color='black')
ax.add_feature(state_boundaries)
ax.add_feature(feat.BORDERS, linewidth='2', edgecolor='black')
ax.set_extent((-118, -73, 23, 50))
# Start the station plot by specifying the axes to draw on, as well as the
# lon/lat of the stations (with transform). We also the fontsize to 12 pt.
stationplot = StationPlot(ax, data['longitude'], data['latitude'],
transform=ccrs.PlateCarree(), fontsize=12)
# The layout knows where everything should go, and things are standardized using
# the names of variables. So the layout pulls arrays out of `data` and plots them
# using `stationplot`.
simple_layout.plot(stationplot, data)
.. parsed-literal::
(Apr 10 16 12:52:00 , Apr 10 16 12:52:00 )
425
85
.. image:: Surface_Obs_Plot_with_MetPy_files/Surface_Obs_Plot_with_MetPy_1_1.png

Binary file not shown.

After

Width:  |  Height:  |  Size: 753 KiB

View file

@ -0,0 +1,269 @@
==========================
Grid Levels and Parameters
==========================
`Notebook <http://nbviewer.ipython.org/github/Unidata/python-awips/blob/master/examples/notebooks/Grid_Levels_and_Parameters.ipynb>`_
List Available Parameters for a Grid Name
-----------------------------------------
.. code:: python
from awips.dataaccess import DataAccessLayer
# Select HRRR
DataAccessLayer.changeEDEXHost("edex-cloud.unidata.ucar.edu")
request = DataAccessLayer.newDataRequest()
request.setDatatype("grid")
request.setLocationNames("GFS40")
# Print parm list
available_parms = DataAccessLayer.getAvailableParameters(request)
available_parms.sort()
for parm in available_parms:
print parm
.. parsed-literal::
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
-----------------------------------
.. code:: python
# Set parm to u-wind
request.setParameters("uW")
# Print level list
available_levels = DataAccessLayer.getAvailableLevels(request)
available_levels.sort()
for level in available_levels:
print level
.. parsed-literal::
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
--------------------------------------------
.. code:: python
import numpy
from metpy.units import units
# Set level for u-wind
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()
# Select v-wind
request.setParameters("vW")
# Select last time for v-wind
response = DataAccessLayer.getGridData(request, [t[-1]])
data_uv = response[-1]
# Print
print 'Time :', t[-1]
print 'Model:', data_uv.getLocationName()
print 'Unit :', data_uv.getUnit()
print 'Parms :', data_uw.getParameter(), data_uv.getParameter()
print data_uv.getRawData().shape
# Calculate total wind speed
spd = numpy.sqrt( data_uw.getRawData()**2 + data_uv.getRawData()**2 )
spd = spd * units.knot
print "windArray =", spd
data = data_uw
.. parsed-literal::
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 ..., 9.390378 9.14996147 8.55599213] [ 8.23072243 8.20243835 8.31557465 ..., 1.48492408 0.56568539 0.39597979] [ 0.49497473 0.52325904 0.1979899 ..., 2.67286372 2.63043714 2.65872145] ..., [ 2.17788887 2.20617294 2.13546252 ..., 1.01823378 0.62225395 0.39597979] [ 0.02828427 0.8768124 1.51320839 ..., 6.47709799 6.68922997 6.84479332] [ 6.92964649 7.02864122 6.98621511 ..., 0.91923875 1.24450791 1.28693426]] knot
Plotting a Grid with Basemap
----------------------------
Using **matplotlib**, **numpy**, and **basemap**:
.. code:: 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.getRawData().shape[1])
y = linspace(0, map.urcrnry, data.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()
.. image:: Grid_Levels_and_Parameters_files/Grid_Levels_and_Parameters_7_0.png
or use **pcolormesh** rather than **contourf**
.. code:: 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')
.. image:: Grid_Levels_and_Parameters_files/Grid_Levels_and_Parameters_9_0.png
Plotting a Grid with Cartopy
----------------------------
.. code:: 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()
.. image:: Grid_Levels_and_Parameters_files/Grid_Levels_and_Parameters_11_0.png

Binary file not shown.

After

Width:  |  Height:  |  Size: 112 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 84 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 105 KiB