mirror of
https://github.com/Unidata/python-awips.git
synced 2025-02-24 06:57:56 -05:00
344 lines
13 KiB
ReStructuredText
344 lines
13 KiB
ReStructuredText
=========================
|
|
Regional Surface Obs Plot
|
|
=========================
|
|
`Notebook <http://nbviewer.ipython.org/github/Unidata/python-awips/blob/master/examples/notebooks/Regional_Surface_Obs_Plot.ipynb>`_
|
|
This exercise creates a surface observsation station plot for the state
|
|
of Florida, using both METAR (datatype *obs*) and Synoptic (datatype
|
|
*sfcobs*). Because we are using the AWIPS Map Database for state and
|
|
county boundaries, there is no use of Cartopy ``cfeature`` in this
|
|
exercise.
|
|
|
|
.. code:: ipython3
|
|
|
|
from awips.dataaccess import DataAccessLayer
|
|
from dynamicserialize.dstypes.com.raytheon.uf.common.time import TimeRange
|
|
from datetime import datetime, timedelta
|
|
import numpy as np
|
|
import cartopy.crs as ccrs
|
|
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
|
|
from cartopy.feature import ShapelyFeature
|
|
from shapely.geometry import Polygon
|
|
import matplotlib.pyplot as plt
|
|
from metpy.units import units
|
|
from metpy.calc import wind_components
|
|
from metpy.plots import simple_layout, StationPlot, StationPlotLayout
|
|
import warnings
|
|
%matplotlib inline
|
|
|
|
def get_cloud_cover(code):
|
|
if 'OVC' in code:
|
|
return 1.0
|
|
elif 'BKN' in code:
|
|
return 6.0/8.0
|
|
elif 'SCT' in code:
|
|
return 4.0/8.0
|
|
elif 'FEW' in code:
|
|
return 2.0/8.0
|
|
else:
|
|
return 0
|
|
|
|
.. code:: ipython3
|
|
|
|
# EDEX request for a single state
|
|
edexServer = "edex-cloud.unidata.ucar.edu"
|
|
DataAccessLayer.changeEDEXHost(edexServer)
|
|
request = DataAccessLayer.newDataRequest('maps')
|
|
request.addIdentifier('table', 'mapdata.states')
|
|
request.addIdentifier('state', 'FL')
|
|
request.addIdentifier('geomField', 'the_geom')
|
|
request.setParameters('state','name','lat','lon')
|
|
response = DataAccessLayer.getGeometryData(request)
|
|
record = response[0]
|
|
print("Found " + str(len(response)) + " MultiPolygon")
|
|
state={}
|
|
state['name'] = record.getString('name')
|
|
state['state'] = record.getString('state')
|
|
state['lat'] = record.getNumber('lat')
|
|
state['lon'] = record.getNumber('lon')
|
|
#state['geom'] = record.getGeometry()
|
|
state['bounds'] = record.getGeometry().bounds
|
|
print(state['name'], state['state'], state['lat'], state['lon'], state['bounds'])
|
|
print()
|
|
|
|
# EDEX request for multiple states
|
|
request = DataAccessLayer.newDataRequest('maps')
|
|
request.addIdentifier('table', 'mapdata.states')
|
|
request.addIdentifier('geomField', 'the_geom')
|
|
request.addIdentifier('inLocation', 'true')
|
|
request.addIdentifier('locationField', 'state')
|
|
request.setParameters('state','name','lat','lon')
|
|
request.setLocationNames('FL','GA','MS','AL','SC','LA')
|
|
response = DataAccessLayer.getGeometryData(request)
|
|
print("Found " + str(len(response)) + " MultiPolygons")
|
|
|
|
# Append each geometry to a numpy array
|
|
states = np.array([])
|
|
for ob in response:
|
|
print(ob.getString('name'), ob.getString('state'), ob.getNumber('lat'), ob.getNumber('lon'))
|
|
states = np.append(states,ob.getGeometry())
|
|
|
|
|
|
.. parsed-literal::
|
|
|
|
Found 1 MultiPolygon
|
|
Florida FL 28.67402 -82.50934 (-87.63429260299995, 24.521051616000022, -80.03199876199994, 31.001012802000048)
|
|
|
|
Found 6 MultiPolygons
|
|
Florida FL 28.67402 -82.50934
|
|
Georgia GA 32.65155 -83.44848
|
|
Louisiana LA 31.0891 -92.02905
|
|
Alabama AL 32.79354 -86.82676
|
|
Mississippi MS 32.75201 -89.66553
|
|
South Carolina SC 33.93574 -80.89899
|
|
|
|
|
|
Now make sure we can plot the states with a lat/lon grid.
|
|
|
|
.. code:: ipython3
|
|
|
|
def make_map(bbox, proj=ccrs.PlateCarree()):
|
|
fig, ax = plt.subplots(figsize=(16,12),subplot_kw=dict(projection=proj))
|
|
ax.set_extent(bbox)
|
|
gl = ax.gridlines(draw_labels=True, color='#e7e7e7')
|
|
gl.top_labels = gl.right_labels = False
|
|
gl.xformatter = LONGITUDE_FORMATTER
|
|
gl.yformatter = LATITUDE_FORMATTER
|
|
return fig, ax
|
|
|
|
# buffer our bounds by +/i degrees lat/lon
|
|
bounds = state['bounds']
|
|
bbox=[bounds[0]-3,bounds[2]+3,bounds[1]-1.5,bounds[3]+1.5]
|
|
|
|
fig, ax = make_map(bbox=bbox)
|
|
shape_feature = ShapelyFeature(states,ccrs.PlateCarree(),
|
|
facecolor='none', linestyle="-",edgecolor='#000000',linewidth=2)
|
|
ax.add_feature(shape_feature)
|
|
|
|
|
|
|
|
|
|
.. parsed-literal::
|
|
|
|
<cartopy.mpl.feature_artist.FeatureArtist at 0x11dcfedd8>
|
|
|
|
|
|
|
|
|
|
.. image:: Regional_Surface_Obs_Plot_files/Regional_Surface_Obs_Plot_4_1.png
|
|
|
|
|
|
--------------
|
|
|
|
Plot METAR (obs)
|
|
----------------
|
|
|
|
Here we use a spatial envelope to limit the request to the boundary or
|
|
our plot. Without such a filter you may be requesting many tens of
|
|
thousands of records.
|
|
|
|
.. code:: ipython3
|
|
|
|
# Create envelope geometry
|
|
envelope = Polygon([(bbox[0],bbox[2]),(bbox[0],bbox[3]),
|
|
(bbox[1], bbox[3]),(bbox[1],bbox[2]),
|
|
(bbox[0],bbox[2])])
|
|
|
|
# New obs request
|
|
DataAccessLayer.changeEDEXHost(edexServer)
|
|
request = DataAccessLayer.newDataRequest("obs", envelope=envelope)
|
|
availableProducts = DataAccessLayer.getAvailableParameters(request)
|
|
single_value_params = ["timeObs", "stationName", "longitude", "latitude",
|
|
"temperature", "dewpoint", "windDir",
|
|
"windSpeed", "seaLevelPress"]
|
|
multi_value_params = ["presWeather", "skyCover", "skyLayerBase"]
|
|
params = single_value_params + multi_value_params
|
|
request.setParameters(*(params))
|
|
|
|
# Time range
|
|
lastHourDateTime = datetime.utcnow() - timedelta(minutes = 60)
|
|
start = lastHourDateTime.strftime('%Y-%m-%d %H:%M:%S')
|
|
end = datetime.utcnow().strftime('%Y-%m-%d %H:%M:%S')
|
|
|
|
beginRange = datetime.strptime( start , "%Y-%m-%d %H:%M:%S")
|
|
endRange = datetime.strptime( end , "%Y-%m-%d %H:%M:%S")
|
|
timerange = TimeRange(beginRange, endRange)
|
|
# Get response
|
|
response = DataAccessLayer.getGeometryData(request,timerange)
|
|
# function getMetarObs was added in python-awips 18.1.4
|
|
obs = DataAccessLayer.getMetarObs(response)
|
|
print("Found " + str(len(response)) + " records")
|
|
print("Using " + str(len(obs['temperature'])) + " temperature records")
|
|
|
|
|
|
.. parsed-literal::
|
|
|
|
Found 3468 records
|
|
Using 152 temperature records
|
|
|
|
|
|
Next grab the simple variables out of the data we have (attaching
|
|
correct units), and put them into a dictionary that we will hand the
|
|
plotting function later:
|
|
|
|
- Get wind components from speed and direction
|
|
- Convert cloud fraction values to integer codes [0 - 8]
|
|
- Map METAR weather codes to WMO codes for weather symbols
|
|
|
|
.. code:: ipython3
|
|
|
|
data = dict()
|
|
data['stid'] = np.array(obs['stationName'])
|
|
data['latitude'] = np.array(obs['latitude'])
|
|
data['longitude'] = np.array(obs['longitude'])
|
|
tmp = np.array(obs['temperature'], dtype=float)
|
|
dpt = np.array(obs['dewpoint'], dtype=float)
|
|
|
|
# Suppress nan masking warnings
|
|
warnings.filterwarnings("ignore",category =RuntimeWarning)
|
|
|
|
tmp[tmp == -9999.0] = 'nan'
|
|
dpt[dpt == -9999.] = 'nan'
|
|
data['air_temperature'] = tmp * units.degC
|
|
data['dew_point_temperature'] = dpt * units.degC
|
|
data['air_pressure_at_sea_level'] = np.array(obs['seaLevelPress'])* units('mbar')
|
|
direction = np.array(obs['windDir'])
|
|
direction[direction == -9999.0] = 'nan'
|
|
u, v = wind_components(np.array(obs['windSpeed']) * units('knots'),
|
|
direction * units.degree)
|
|
data['eastward_wind'], data['northward_wind'] = u, v
|
|
data['cloud_coverage'] = [int(get_cloud_cover(x)*8) for x in obs['skyCover']]
|
|
data['present_weather'] = obs['presWeather']
|
|
proj = ccrs.LambertConformal(central_longitude=state['lon'], central_latitude=state['lat'],
|
|
standard_parallels=[35])
|
|
custom_layout = StationPlotLayout()
|
|
custom_layout.add_barb('eastward_wind', 'northward_wind', units='knots')
|
|
custom_layout.add_value('NW', 'air_temperature', fmt='.0f', units='degF', color='darkred')
|
|
custom_layout.add_value('SW', 'dew_point_temperature', fmt='.0f', units='degF', color='darkgreen')
|
|
custom_layout.add_value('E', 'precipitation', fmt='0.1f', units='inch', color='blue')
|
|
ax.set_title(str(response[-1].getDataTime()) + " | METAR Surface Obs | " + edexServer)
|
|
stationplot = StationPlot(ax, data['longitude'], data['latitude'], clip_on=True,
|
|
transform=ccrs.PlateCarree(), fontsize=10)
|
|
custom_layout.plot(stationplot, data)
|
|
fig
|
|
|
|
|
|
|
|
|
|
.. image:: Regional_Surface_Obs_Plot_files/Regional_Surface_Obs_Plot_8_0.png
|
|
|
|
|
|
|
|
--------------
|
|
|
|
Plot Synoptic (sfcobs)
|
|
----------------------
|
|
|
|
.. code:: ipython3
|
|
|
|
# New sfcobs/SYNOP request
|
|
DataAccessLayer.changeEDEXHost(edexServer)
|
|
request = DataAccessLayer.newDataRequest("sfcobs", envelope=envelope)
|
|
availableProducts = DataAccessLayer.getAvailableParameters(request)
|
|
# (sfcobs) uses stationId, while (obs) uses stationName,
|
|
# the rest of these parameters are the same.
|
|
single_value_params = ["timeObs", "stationId", "longitude", "latitude",
|
|
"temperature", "dewpoint", "windDir",
|
|
"windSpeed", "seaLevelPress"]
|
|
multi_value_params = ["presWeather", "skyCover", "skyLayerBase"]
|
|
pres_weather, sky_cov, sky_layer_base = [],[],[]
|
|
params = single_value_params + multi_value_params
|
|
request.setParameters(*(params))
|
|
|
|
# Time range
|
|
lastHourDateTime = datetime.utcnow() - timedelta(minutes = 60)
|
|
start = lastHourDateTime.strftime('%Y-%m-%d %H:%M:%S')
|
|
end = datetime.utcnow().strftime('%Y-%m-%d %H:%M:%S')
|
|
|
|
beginRange = datetime.strptime( start , "%Y-%m-%d %H:%M:%S")
|
|
endRange = datetime.strptime( end , "%Y-%m-%d %H:%M:%S")
|
|
timerange = TimeRange(beginRange, endRange)
|
|
|
|
# Get response
|
|
response = DataAccessLayer.getGeometryData(request,timerange)
|
|
# function getSynopticObs was added in python-awips 18.1.4
|
|
sfcobs = DataAccessLayer.getSynopticObs(response)
|
|
print("Found " + str(len(response)) + " records")
|
|
print("Using " + str(len(sfcobs['temperature'])) + " temperature records")
|
|
|
|
|
|
.. parsed-literal::
|
|
|
|
Found 260 records
|
|
Using 78 temperature records
|
|
|
|
|
|
.. code:: ipython3
|
|
|
|
data = dict()
|
|
data['stid'] = np.array(sfcobs['stationId'])
|
|
data['lat'] = np.array(sfcobs['latitude'])
|
|
data['lon'] = np.array(sfcobs['longitude'])
|
|
|
|
# Synop/sfcobs temps are stored in kelvin (degC for METAR/obs)
|
|
tmp = np.array(sfcobs['temperature'], dtype=float)
|
|
dpt = np.array(sfcobs['dewpoint'], dtype=float)
|
|
direction = np.array(sfcobs['windDir'])
|
|
# Account for missing values
|
|
tmp[tmp == -9999.0] = 'nan'
|
|
dpt[dpt == -9999.] = 'nan'
|
|
direction[direction == -9999.0] = 'nan'
|
|
|
|
data['air_temperature'] = tmp * units.kelvin
|
|
data['dew_point_temperature'] = dpt * units.kelvin
|
|
data['air_pressure_at_sea_level'] = np.array(sfcobs['seaLevelPress'])* units('mbar')
|
|
try:
|
|
data['eastward_wind'], data['northward_wind'] = wind_components(
|
|
np.array(sfcobs['windSpeed']) * units('knots'),direction * units.degree)
|
|
data['present_weather'] = sfcobs['presWeather']
|
|
except ValueError:
|
|
pass
|
|
|
|
fig_synop, ax_synop = make_map(bbox=bbox)
|
|
shape_feature = ShapelyFeature(states,ccrs.PlateCarree(),
|
|
facecolor='none', linestyle="-",edgecolor='#000000',linewidth=2)
|
|
ax_synop.add_feature(shape_feature)
|
|
|
|
custom_layout = StationPlotLayout()
|
|
custom_layout.add_barb('eastward_wind', 'northward_wind', units='knots')
|
|
custom_layout.add_value('NW', 'air_temperature', fmt='.0f', units='degF', color='darkred')
|
|
custom_layout.add_value('SW', 'dew_point_temperature', fmt='.0f', units='degF', color='darkgreen')
|
|
custom_layout.add_value('E', 'precipitation', fmt='0.1f', units='inch', color='blue')
|
|
ax_synop.set_title(str(response[-1].getDataTime()) + " | SYNOP Surface Obs | " + edexServer)
|
|
stationplot = StationPlot(ax_synop, data['lon'], data['lat'], clip_on=True,
|
|
transform=ccrs.PlateCarree(), fontsize=10)
|
|
custom_layout.plot(stationplot, data)
|
|
|
|
|
|
|
|
.. image:: Regional_Surface_Obs_Plot_files/Regional_Surface_Obs_Plot_11_0.png
|
|
|
|
|
|
--------------
|
|
|
|
Plot both METAR and SYNOP
|
|
-------------------------
|
|
|
|
.. code:: ipython3
|
|
|
|
custom_layout = StationPlotLayout()
|
|
custom_layout.add_barb('eastward_wind', 'northward_wind', units='knots')
|
|
custom_layout.add_value('NW', 'air_temperature', fmt='.0f', units='degF', color='darkred')
|
|
custom_layout.add_value('SW', 'dew_point_temperature', fmt='.0f', units='degF', color='darkgreen')
|
|
custom_layout.add_value('E', 'precipitation', fmt='0.1f', units='inch', color='blue')
|
|
ax.set_title(str(response[-1].getDataTime()) + " | METAR/SYNOP Surface Obs | " + edexServer)
|
|
stationplot = StationPlot(ax, data['lon'], data['lat'], clip_on=True,
|
|
transform=ccrs.PlateCarree(), fontsize=10)
|
|
custom_layout.plot(stationplot, data)
|
|
fig
|
|
|
|
|
|
|
|
|
|
.. image:: Regional_Surface_Obs_Plot_files/Regional_Surface_Obs_Plot_13_0.png
|
|
|
|
|