import os
from datetime import datetime, timedelta
import numpy as np
import numpy.ma as ma
import matplotlib.patheffects as mpatheffects
import matplotlib.pyplot as plt
import urllib
from urllib.request import urlopen
from urllib.request import urlretrieve
from metpy.plots import StationPlot, add_metpy_logo, add_unidata_logo, add_timestamp
from metpy.calc import wind_components
from metpy.plots import StationPlot, StationPlotLayout, wx_code_map, simple_layout
from metpy.plots.wx_symbols import sky_cover, current_weather
from metpy.units import units
import metpy.calc as mpcalc
from siphon.catalog import TDSCatalog
from siphon.simplewebservice.ndbc import NDBC
import pandas
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import cartopy.io.shapereader as cshaperead
import pyart
from awips.dataaccess import DataAccessLayer
from dynamicserialize.dstypes.com.raytheon.uf.common.time import TimeRange
import warnings
warnings.filterwarnings('ignore')
def rebin(a, shape):
sh = shape[0],a.shape[0]//shape[0],shape[1],a.shape[1]//shape[1]
return a.reshape(sh).mean(-1).mean(1)
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
now = datetime.now()
datetime_str = datetime.now().strftime('%Y %b %d - %H:%M UTC')
#USER INPUT
goes_channel = "02"
radar_str = "KIWX"
bound_box = [-90, -84, 40, 43]
central_latitude = 45.
central_longitude = -87.
print("Imports and User Input Complete.")
## You are using the Python ARM Radar Toolkit (Py-ART), an open source ## library for working with weather radar data. Py-ART is partly ## supported by the U.S. Department of Energy as part of the Atmospheric ## Radiation Measurement (ARM) Climate Research Facility, an Office of ## Science user facility. ## ## If you use this software to prepare a publication, please cite: ## ## JJ Helmus and SM Collis, JORS 2016, doi: 10.5334/jors.119 Imports and User Input Complete.
print("Retrieving GOES Channel " + goes_channel + " satellite data...",end="")
cat1 = TDSCatalog('http://thredds-test.unidata.ucar.edu/thredds/catalog/satellite/goes16/GOES16/CONUS/Channel'+goes_channel+'/current/catalog.xml')
sat1 = cat1.datasets[-1].remote_access(use_xarray=True)
var1 = sat1.metpy.parse_cf('Sectorized_CMI')
ref_gamma_1 = np.sqrt(var1.values)
print("Done.")
Retrieving GOES Channel 02 satellite data...Done.
print("Retrieving radar data...",end="")
YYYYMMDD = datetime.now().strftime("%Y%m%d")
directory = 'http://thredds.ucar.edu/thredds/catalog/nexrad/level2/'+radar_str+'/'+YYYYMMDD
cat = TDSCatalog(directory+'/catalog.html')
filename = str(cat.datasets[1])
#radar_datetime = datetime.strptime(filename.split('_')[2]+filename.split('_')[3][0:-5], '%Y%m%d%H%M')
http_directory = 'http://thredds.ucar.edu/thredds/fileServer/nexrad/level2/'+radar_str+'/'+YYYYMMDD
file2download = http_directory+'/'+filename
content = urlretrieve(file2download,filename) #download radar data
radar = pyart.io.read(filename)
display = pyart.graph.RadarMapDisplayCartopy(radar)
if os.path.exists(filename): #remove radar data from computer
os.remove(filename)
print("Done.")
Retrieving radar data...Done.
print("Retrieving METAR data...",end="")
single_value_params = ["timeObs", "stationName", "longitude", "latitude",
"temperature", "dewpoint", "windDir",
"windSpeed", "seaLevelPress"]
multi_value_params = ["presWeather", "skyCover", "skyLayerBase"]
all_params = single_value_params + multi_value_params
obs_dict = dict({all_params: [] for all_params in all_params})
pres_weather = []
sky_cov = []
sky_layer_base = []
LastHourDateTime = datetime.utcnow() - timedelta(hours = 3)
NowDateTime = datetime.utcnow()
timerange = TimeRange(LastHourDateTime, NowDateTime)
DataAccessLayer.changeEDEXHost("edex-cloud.unidata.ucar.edu")
request = DataAccessLayer.newDataRequest()
request.setDatatype("obs")
request.setParameters(*(all_params))
response = DataAccessLayer.getGeometryData(request,timerange)
print("Connected to Server...",end="")
for ob in response:
avail_params = ob.getParameters()
if "presWeather" in avail_params:
pres_weather.append(ob.getString("presWeather"))
elif "skyCover" in avail_params and "skyLayerBase" in avail_params:
sky_cov.append(ob.getString("skyCover"))
sky_layer_base.append(ob.getNumber("skyLayerBase"))
else:
for param in single_value_params:
if param in avail_params:
if param == 'timeObs':
obs_dict[param].append(datetime.fromtimestamp(ob.getNumber(param)/1000.0))
else:
try:
obs_dict[param].append(ob.getNumber(param))
except TypeError:
obs_dict[param].append(ob.getString(param))
else:
obs_dict[param].append(None)
obs_dict['presWeather'].append(pres_weather);
obs_dict['skyCover'].append(sky_cov);
obs_dict['skyLayerBase'].append(sky_layer_base);
pres_weather = []
sky_cov = []
sky_layer_base = []
print("Using Pandas...",end="")
df = pandas.DataFrame(data=obs_dict, columns=all_params)
#sort rows with the newest first
df = df.sort_values(by='timeObs', ascending=False)
#group rows by station
groups = df.groupby('stationName')
#create a new DataFrame for the most recent values
df_recent = pandas.DataFrame(columns=all_params)
#retrieve the first entry for each group, which will
#be the most recent observation
for rid, station in groups:
row = station.head(1)
df_recent = pandas.concat([df_recent, row])
print("Saving to Dictionary...",end="")
lats = np.asarray(df_recent['latitude'])
lons = np.asarray(df_recent['longitude'])
inds = np.where((lats>bound_box[2]) & (lats<bound_box[3]) &
(lons>bound_box[0]) & (lons<bound_box[1]))[0]
data = dict()
data['stid'] = np.array(df_recent["stationName"])[inds]
data['latitude'] = np.array(df_recent['latitude'])[inds]
data['longitude'] = np.array(df_recent['longitude'])[inds]
data['air_temperature'] = (np.array(df_recent['temperature'], dtype=float)[inds] * units.degC).to('degF')
data['air_temperature'] = ma.masked_less(data['air_temperature'], -9999.* units.degF)
data['dew_point'] = (np.array(df_recent['dewpoint'], dtype=float)[inds] * units.degC).to('degF')
data['dew_point'][np.where((np.array(data['dew_point'])<-50) | (np.array(data['dew_point'])>100))[0]] = np.nan
data['slp'] = np.array(df_recent['seaLevelPress'])[inds] * units('mbar')
#data['slp'] = ma.masked_values(data['slp'],-9999.*units('mbar'))
data['slp'][np.where(np.array(data['slp'])==-9999)[0]] = np.nan
u, v = wind_components(np.array(df_recent['windSpeed'])[inds] * units('knots'),
np.array(df_recent['windDir'])[inds] * units.degree)
data['eastward_wind'], data['northward_wind'] = u, v
data['cloud_frac'] = np.array([int(get_cloud_cover(x)*8) for x in df_recent['skyCover']])[inds]
data['presWeather'] = np.array([i[0] for i in np.array(df_recent['presWeather'])])[inds]
wx = []
for i in range(len(data['presWeather'])):
try:
wx.append(wx_code_map[data['presWeather'][i]])
except:
wx.append('')
stations = np.asarray(obs_dict['stationName'])[inds]
print("Done. "+ str(len(stations))+ " stations.")
Retrieving METAR data...Connected to Server...Using Pandas...Saving to Dictionary...Done. 57 stations.
print("Pulling in buoy data...",end="")
buoy = NDBC.latest_observations()
buoy.dropna(subset=['pressure', 'wind_speed', 'wind_direction'], inplace=True)
buoy_u, buoy_v = mpcalc.wind_components(buoy['wind_speed'].values, buoy['wind_direction'].values * units.degree)
xy = var1.metpy.cartopy_crs.transform_points(ccrs.PlateCarree(), buoy['longitude'].values, buoy['latitude'].values)
mask = mpcalc.reduce_point_density(xy[..., 0:2], 30000, priority=buoy['pressure'].values)
lats = np.asarray(buoy['latitude'].values)
lons = np.asarray(buoy['longitude'].values)
inds_buoy = np.where((lats>bound_box[2]) & (lats<bound_box[3]) &
(lons>bound_box[0]) & (lons<bound_box[1]))[0]
print("Done. "+ str(len(buoy_u))+ " stations.")
Pulling in buoy data...Done. 435 stations.
#%matplotlib inline
print("*Begin Plotting*")
print("Map Projection: Lambert Conformal.")
proj = ccrs.LambertConformal(central_longitude=central_longitude, central_latitude=central_latitude)
fig = plt.figure(figsize=(10, 8))#, dpi=100)
ax = fig.add_subplot(1, 1, 1, projection=proj)
print("Plotting satellite data...",end="")
ax.pcolormesh(var1.x.values, var1.y.values, ref_gamma_1, cmap='Greys_r', transform=var1.metpy.cartopy_crs)
print("Done.")
print("plotting radar data...",end="")
field='reflectivity'
title = 'TEST'
offset = 0.75
display.plot_ppi_map(field, 0, ax=ax,vmin=5, vmax=80,embelish=False, title_flag=False, colorbar_flag=False, cmap = pyart.graph.cm.NWSRef)
print("Done.")
print("plotting ASAS station data...",end="")
kwargs = dict(path_effects=[mpatheffects.withStroke(foreground='white', linewidth=1)], clip_on=True)
sp = StationPlot(ax, data['longitude'], data['latitude'], transform=ccrs.PlateCarree(), fontsize=10)
simple_layout.plot(sp, data)
sp.plot_barb(data['eastward_wind'], data['northward_wind'], edgecolor='k',zorder=2, clip_on=True)
#sp.plot_text('SE', data["stid"], fontsize=8,zorder=3, **kwargs)
sp.plot_parameter('NE', data['slp'], color='black', fontsize=8,zorder=4,
formatter=lambda v: format( v/100., '.0f')[-13:-9], **kwargs)
sp.plot_parameter('SW', data['dew_point'], color='darkgreen', fontsize=8,zorder=5, **kwargs)
sp.plot_parameter('NW', data['air_temperature'], color='red', fontsize=8,zorder=6, **kwargs)
sp.plot_symbol('W', wx, current_weather, color='m', clip_on=True, zorder=7)
#sp.plot_symbol('C', data['cloud_frac'], sky_cover,zorder=5,clip_on=True)
print("Done.")
print("Plotting buoy data...",end="")
kwargs = dict(path_effects=[mpatheffects.withStroke(foreground='white', linewidth=1)], clip_on=True)
sp = StationPlot(ax, buoy['longitude'].values[inds_buoy], buoy['latitude'].values[inds_buoy], spacing=8, transform=ccrs.PlateCarree())
sp.plot_barb(buoy_u[inds_buoy], buoy_v[inds_buoy], edgecolor='k',zorder=2, clip_on=True)
sp.plot_parameter('NE', buoy['pressure'].values[inds_buoy], color='black', fontsize=7,zorder=4, **kwargs)
sp.plot_parameter('SW', (9./5.)*(buoy['dewpoint'].values[inds_buoy])+32., color='darkgreen', fontsize=8,zorder=5, **kwargs)
sp.plot_parameter('NW', (9./5.)*(buoy['air_temperature'].values[inds_buoy])+32., color='red', fontsize=8,zorder=6, **kwargs)
print("Done.")
print("Adding map features...",end="")
ax.add_feature(cfeature.COASTLINE.with_scale('10m'), edgecolor='white')
ax.add_feature(cfeature.STATES.with_scale('10m'), edgecolor='white')
ax.set_extent([-90, -84, 40, 43], crs=ccrs.PlateCarree())
print("Done.")
print("Adding radar colorbar...",end="")
cint_refl = 5.0
cbar_refl_min = 5
cbar_refl_max = 80 + (cint_refl)
cflevs_refl = np.arange(cbar_refl_min, cbar_refl_max, cint_refl)
refl_ticks = np.arange(cbar_refl_min, cbar_refl_max, cint_refl*1.)
blah1 = np.random.rand(3,2)*10.
blah2 = np.random.rand(3,2)*10.
blah3 = np.random.rand(3,2)*10.
CS1 = plt.contourf(blah1,blah2,blah3,levels=cflevs_refl, cmap = pyart.graph.cm.NWSRef, extend='max')
cbar_ax = fig.add_axes([0.91, 0.20, 0.05, 0.45],frameon=False)
cbar = plt.colorbar(CS1, pad=0.15, fraction=0.9, orientation='vertical', extend='max')
cbar_ax.set_xticks([])
cbar_ax.set_yticks([])
cbar.set_ticks(refl_ticks)
cbar.set_label('Radar Reflectivity (dBZ)',size=10)
cbar.ax.tick_params(labelsize=8)
print("Done.")
print("Add logos, titles, etc...",end="")
title_text = 'GOES-16 Channel ' + goes_channel + ' - ' + radar_str + ' - NWS ASOS - NDBC Buoy'
ax.set_title(title_text)
ax.annotate(datetime_str, xy=(0.90, 1.01), xytext=(-15, -15), fontsize=9,
xycoords='axes fraction', textcoords='offset points',
bbox=dict(facecolor='#F1F0F2', edgecolor='black', alpha=1.0, boxstyle='round,pad=0.5'),
horizontalalignment='center', verticalalignment='top',zorder=8)
#add_unidata_logo(fig, y=2000, x=3250, size='large')
#add_metpy_logo(fig, y=1825, x=3250, size='large')
#ARMlogo = plt.imread('ARM_Logo.png')
#ax.figure.figimage(ARMlogo, 3220, 1720)
#ax.annotate('PyART', xy=(1.085, 0.868), xytext=(-15, -15), fontsize=7,
# xycoords='axes fraction', textcoords='offset points',
# horizontalalignment='center', verticalalignment='top',zorder=8)
print("Done.")
print("Saving file...",end="")
plt.savefig('/home/dgoines/Desktop/MetPy_blog.png', bbox_inches='tight', dpi=400)
print("Done.")
*Begin Plotting* Map Projection: Lambert Conformal. Plotting satellite data...Done. plotting radar data...Done. plotting ASAS station data...Done. Plotting buoy data...Done. Adding map features...Done. Adding radar colorbar...Done. Add logos, titles, etc...Done. Saving file...Done.