In [1]:
#fetch GFS data and plot at ENA
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import netCDF4
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
import numpy as np
from siphon.catalog import TDSCatalog
from siphon.ncss import NCSS
from datetime import datetime, timedelta
from netCDF4 import num2date
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import os
import ftplib
import pygrib
import tempfile
import xarray
import boto3
%matplotlib inline
In [123]:
import imp
lib_loc = os.path.join(os.path.expanduser('~'), 'projects/ACE-ENA-EVA/code/ena_tools.py')
ena_tools = imp.load_source('ena_tools', lib_loc)
In [3]:
def get_ecmwf_137():
levels_loc = os.path.join(os.path.expanduser('~'), 'projects/ACE-ENA-EVA/ecmwf_137_levels.txt')
levels = np.genfromtxt(levels_loc, missing_values='-')
ht = levels[1::, 6]
pres = levels[1::, 3]
return ht, pres
def extract_3d_grib(grb_obj, search_term, skip_last=False):
grb_list = grb_obj.select(name=search_term)
level_nums = [this_grb['level'] for this_grb in grb_list]
order = np.array(level_nums).argsort()
lats, lons = grb_list[0].latlons()
shp = grb_list[order[0]].values.shape
transfer_array = np.empty([len(order), shp[0], shp[1]])
if skip_last:
llen = len(order) -1
else:
llen = len(order)
for i in range(llen):
transfer_array[i,:,:] = grb_list[order[i]].values
return transfer_array, lats, lons
def ecmwf_name_to_date(ename):
start_time = datetime.strptime('2017'+ename[3:11], '%Y%m%d%H%M')
valid_time = datetime.strptime('2017'+ename[11:19], '%Y%m%d%H%M')
return start_time, valid_time
def file_list_to_date(file_list):
start_times = [ecmwf_name_to_date(this_name)[0] for this_name in file_list]
end_times = [ecmwf_name_to_date(this_name)[1] for this_name in file_list]
return start_times, end_times
def get_run_hours(file_list):
gen_t, val_t = file_list_to_date(file_list)
gen_hour = np.array([dt.hour for dt in gen_t])
return np.unique(gen_hour)
def get_time_for_run(file_list, gen_hour):
gen_t, val_t = file_list_to_date(file_list)
gen_hours = np.array([dt.hour for dt in gen_t])
good_files = []
good_times_gen = []
good_times_val = []
for i in range(len(gen_hours)):
if gen_hours[i] == gen_hour:
good_files.append(file_list[i])
good_times_gen.append(gen_t[i])
good_times_val.append(val_t[i])
return good_files, good_times_gen, good_times_val
def create_bundle_latest(var_list, n=None, skippy = None):
username_file = os.path.join(os.path.expanduser('~'), 'ecmwf_username')
password_file = os.path.join(os.path.expanduser('~'), 'ecmwf_passwd')
uname_fh = open(username_file, 'r')
uname = uname_fh.readline()[0:-1]
uname_fh.close()
passwd_fh = open(password_file, 'r')
passwd = passwd_fh.readline()[0:-1]
passwd_fh.close()
host = 'dissemination.ecmwf.int'
#get ECMWF vert coord
ht, pres = get_ecmwf_137()
ftp = ftplib.FTP(host)
ftp.login(user=uname, passwd = passwd )
closest_now = datetime.utcnow().strftime('%Y%m%d')
ftp.cwd(closest_now)
lst = ftp.nlst()
lst.sort()
run_hours = get_run_hours(lst)
target_files, generated_times, valid_times = get_time_for_run(lst, run_hours[-1])
if skippy is None:
skippy = len(var_list)*[False]
if n is None:
these_target_files = target_files[1::]
these_run_times = generated_times[1::]
these_valid_times = valid_times[1::]
else:
these_target_files = target_files[1:n]
these_run_times = generated_times[1:n]
these_valid_times = valid_times[1:n]
bundle = {}
for var_name in var_list:
bundle.update({var_name: []})
for i in range(len(these_valid_times)):
print(these_target_files[i])
fh = tempfile.NamedTemporaryFile()
ftp.retrbinary('RETR ' + these_target_files[i], fh.write)
grbs = pygrib.open(fh.name)
grbs.seek(0)
for i in range(len(var_list)):
var_name = var_list[i]
this_step, lats, lons = extract_3d_grib(grbs, var_name, skip_last = skippy[i])
bundle[var_name].append(this_step)
grbs.close()
return bundle, these_valid_times, these_run_times, lats, lons
def concat_bundle(bundle):
varss = list(bundle.keys())
n_times = len(bundle[varss[0]])
hlatlon = bundle[varss[0]][0].shape
unwound = {}
for this_var in varss:
transfer = np.empty([n_times, 136, hlatlon[1], hlatlon[2]])
for time_step in range(n_times):
transfer[time_step, 0:136, :, :] = bundle[this_var][time_step][0:136, :, :]
unwound.update({this_var: transfer})
return unwound
def unwind_to_xarray(unwound, valid_times, lats, lons):
ds = xarray.Dataset()
for vvar in list(unwound.keys()):
my_data = xarray.DataArray(unwound[vvar],
dims = ('time', 'z', 'y', 'x'),
coords = {'time' : (['time'], valid_times),
'z' : (['z'], get_ecmwf_137()[0][0:136]),
'lat' :(['y','x'], lats),
'lon' : (['y','x'],my_lons)})
ds[vvar.replace(' ', '_')] = my_data
ds.lon.attrs = [('long_name', 'longitude of grid cell center'),
('units', 'degrees_east'),
('bounds', 'xv')]
ds.lat.attrs = [('long_name', 'latitude of grid cell center'),
('units', 'degrees_north'),
('bounds', 'yv')]
return ds
def save_one_ecmwf_clouds(dataset, gen_datetime):
start_str = gen_datetime.strftime('%Y%m%d_%H%M')
s3name = 'ecmwf_time_height_clwc_ciwc'
plt.figure(figsize=(17,5))
my_levels = [0.01, 0.05, 0.09, 0.13]
my_colors = ['white', 'yellow', 'cyan', 'pink']
(my_dataset.Specific_cloud_liquid_water_content*1000.0).mean(dim=('y','x')).plot.pcolormesh(x='time', y='z')
cs = (my_dataset.Specific_cloud_ice_water_content*1000.0).max(dim=('y','x')).plot.contour(x='time', y='z',
levels=my_levels,
colors=my_colors)
plt.clabel(cs, inline=1, fontsize=10, fmt='%0.3f')
plt.ylim([0,12000])
str1 = start_str + ' ECMWF Domain max cloud ice and domain mean cloud liquid water content (g/kg) \n'
str2 = 'ACE-ENA forecast guidence. ARM Climate Research Facility. scollis@anl.gov'
plt.title(str1+str2)
local_fig = tempfile.NamedTemporaryFile(suffix='.png')
fn = local_fig.name
plt.savefig(fn)
s3_key = s3name + '/' + ena_tools.gen_s3_key(gen_datetime, s3name)
s3 = boto3.resource('s3')
data = open(fn, 'rb')
s3.Bucket('aceena').put_object(Key=s3_key, Body=data, ACL='public-read')
data.close()
data2 = open(fn, 'rb')
s3.Bucket('aceena').put_object(Key='latest_' + s3name + '.png',
Body=data2, ACL='public-read')
data2.close()
def save_one_ecmwf_humidity(dataset, gen_datetime):
start_str = gen_datetime.strftime('%Y%m%d_%H%M')
s3name = 'ecmwf_time_height_clwc_ciwc'
plt.figure(figsize=(17,5))
my_levels = [0.01, 0.05, 0.09, 0.13]
my_colors = ['white', 'yellow', 'cyan', 'pink']
(my_dataset.Specific_humidity*1000.0).mean(dim=('y','x')).plot.pcolormesh(x='time', y='z')
cs = (my_dataset.Specific_cloud_ice_water_content*1000.0).max(dim=('y','x')).plot.contour(x='time', y='z',
levels=my_levels,
colors=my_colors)
plt.clabel(cs, inline=1, fontsize=10, fmt='%0.3f')
plt.ylim([0,12000])
str1 = start_str + ' ECMWF Domain max cloud ice and domain mean cloud liquid water content (g/kg) \n'
str2 = 'ACE-ENA forecast guidence. ARM Climate Research Facility. scollis@anl.gov'
plt.title(str1+str2)
local_fig = tempfile.NamedTemporaryFile(suffix='.png')
fn = local_fig.name
plt.savefig(fn)
s3_key = s3name + '/' + ena_tools.gen_s3_key(gen_datetime, s3name)
s3 = boto3.resource('s3')
data = open(fn, 'rb')
s3.Bucket('aceena').put_object(Key=s3_key, Body=data, ACL='public-read')
data.close()
data2 = open(fn, 'rb')
s3.Bucket('aceena').put_object(Key='latest_' + s3name + '.png',
Body=data2, ACL='public-read')
data2.close()
In [5]:
#https://www.ecmwf.int/en/forecasts/documentation-and-support/data-delivery/manage-transmission-ecpds/real-time-data-file
In [19]:
username_file = os.path.join(os.path.expanduser('~'), 'ecmwf_username')
password_file = os.path.join(os.path.expanduser('~'), 'ecmwf_passwd')
uname_fh = open(username_file, 'r')
uname = uname_fh.readline()[0:-1]
uname_fh.close()
passwd_fh = open(password_file, 'r')
passwd = passwd_fh.readline()[0:-1]
passwd_fh.close()
host = 'dissemination.ecmwf.int'
ftp = ftplib.FTP(host)
ftp.login(user=uname, passwd = passwd )
closest_now = datetime.utcnow().strftime('%Y%m%d')
ftp.cwd(closest_now)
lst = ftp.nlst()
lst.sort()
run_hours = get_run_hours(lst)
target_files, generated_times, valid_times = get_time_for_run(lst, run_hours[-1])
fh = tempfile.NamedTemporaryFile()
ftp.retrbinary('RETR ' + target_files[1], fh.write)
grbs = pygrib.open(fh.name)
In [21]:
grb_list = grbs.select(name='Temperature')
In [29]:
def extract_3d_grib(grb_obj, search_term, skip_last=False,
metadata_list = None):
if metadata_list is None:
metadata_list = ['units', 'cfName']
grb_list = grb_obj.select(name=search_term)
level_nums = [this_grb['level'] for this_grb in grb_list]
order = np.array(level_nums).argsort()
lats, lons = grb_list[0].latlons()
shp = grb_list[order[0]].values.shape
metadata = {}
for nm in metadata_list:
metadata.update({nm: grb_list[order[0]][nm]})
transfer_array = np.empty([len(order), shp[0], shp[1]])
if skip_last:
llen = len(order) -1
else:
llen = len(order)
for i in range(llen):
transfer_array[i,:,:] = grb_list[order[i]].values
return transfer_array, lats, lons, metadata
In [30]:
t, ll, ln, md = extract_3d_grib(grbs, 'Temperature')
In [31]:
md
Out[31]:
In [28]:
grb1 = grb_list[1]
print(grb1.keys())
print(grb1['cfName'])
In [ ]:
10:Specific rain water content:kg kg**-1 (instant):regular_ll:hybrid:level 1:fcst time 57 hrs:from 201706260000
11:Specific snow water content:kg kg**-1 (instant):regular_ll:hybrid:level 1:fcst time 57 hrs:from 201706260000
12:Temperature:K (instant):regular_ll:hybrid:level 1:fcst time 57 hrs:from 201706260000
13:V component of wind:m s**-1 (instant):regular_ll:hybrid:level 1:fcst time 57 hrs:from 201706260000
14:U component of wind:m s**-1 (instant):regular_ll:hybrid:level 1:fcst time 57 hrs:from 201706260000
15:Specific humidity:kg kg**-1 (instant):regular_ll:hybrid:level 1:fcst time 57 hrs:from 201706260000
16:Vertical velocity:Pa s**-1 (instant):regular_ll:hybrid:level 1:fcst time 57 hrs:from 201706260000
21:Specific cloud liquid water content:kg kg**-1 (instant):regular_ll:hybrid:level 1:fcst time 57 hrs:from 201706260000
22:Specific cloud ice water content:kg kg**-1 (instant):regular_ll:hybrid:level 1:fcst time 57 hrs:from 201706260000
In [16]:
var_list = ['Specific cloud liquid water content', 'Specific cloud ice water content',
'Specific rain water content', 'Specific snow water content',
'Temperature', 'V component of wind', 'U component of wind', 'Specific humidity', 'Vertical velocity']
skip_me = [False, False, False, False, False, True, True, False, False]
In [113]:
my_bundle, my_these_valid_times, my_these_run_times, my_lats, my_lons, metad = ena_tools.create_bundle_latest(var_list,
skippy=skip_me,
n=20)
In [129]:
import imp
lib_loc = os.path.join(os.path.expanduser('~'), 'projects/ACE-ENA-EVA/code/ena_tools.py')
ena_tools = imp.load_source('ena_tools', lib_loc)
In [130]:
my_unwind = ena_tools.concat_bundle(my_bundle)
trans = {'cfName': 'standard_name'}
my_dataset = ena_tools.unwind_to_xarray(my_unwind, my_these_valid_times, my_lats, my_lons, metad, trans=trans)
In [131]:
my_dataset.attrs['conventions'] = 'CF 1.6'
my_dataset.attrs['source'] = 'ECMWF 137 level 0.1 degree model'
my_dataset.attrs['conatact'] = 'Scott Collis, scollis@anl.gov'
st1 = "European Center for Medium range Weather Forecasting, "
st2 = "Atmospheric Climate Research Facility, A DoE User facility, "
st3 = "Scott Collis, Argonne National Laboratory"
my_dataset.attrs['attribution'] = st1 + st2 + st3
my_dataset.attrs['experiment'] ='ACE-ENA'
my_dataset.Specific_cloud_ice_water_content.attrs['standard_name'] = 'cloud_ice_mixing_ratio'
my_dataset.Specific_cloud_liquid_water_content.attrs['standard_name'] = 'cloud_liquid_water_mixing_ratio'
my_dataset.Specific_rain_water_content.attrs['standard_name'] = 'rain_water_content'
my_dataset.Specific_snow_water_content.attrs['standard_name'] = 'snow_water_content'
my_dataset.encoding['unlimited_dims'] = ['time']
my_dataset.z.encoding['_FillValue'] = -9999
my_dataset.lat.encoding['_FillValue'] = -9999
my_dataset.lon.encoding['_FillValue'] = -9999
In [ ]:
In [132]:
my_dataset
Out[132]:
In [ ]:
In [133]:
my_dataset.to_netcdf('/data/hard_stools.nc')
In [389]:
ds.Specific_cloud_liquid_water_content.mean(dim=('y','x')).plot?
In [392]:
my_dataset.Specific_cloud_liquid_water_content.mean(dim=('y','x')).plot()
Out[392]:
In [450]:
plt.figure(figsize=(17,5))
my_levels = [0.01, 0.05, 0.09, 0.13]
my_colors = ['white', 'yellow', 'cyan', 'pink']
(my_dataset.Specific_cloud_liquid_water_content*1000.0).mean(dim=('y','x')).plot.pcolormesh(x='time', y='z')
cs = (my_dataset.Specific_cloud_ice_water_content*1000.0).max(dim=('y','x')).plot.contour(x='time', y='z',
levels=my_levels,
colors=my_colors)
plt.clabel(cs, inline=1, fontsize=10, fmt='%0.3f')
plt.ylim([0,12000])
Out[450]:
In [55]:
pres2d.shape
Out[55]:
In [54]:
$RH=100wws≈0.263pq[exp(17.67(T−T0)T−29.65)]−1$
In [58]:
pres2d = np.tile(pres[0:136], [my_dataset.Specific_humidity.shape[0],1])
rh = my_dataset.Specific_humidity.mean(dim=('y','x')) * 0.263 * pres2d * \
np.exp((17.67*(my_dataset.Temperature.mean(dim=('y','x')) - 273.15))/(my_dataset.Temperature.mean(dim=('y','x'))-29.65))**(-1)
In [61]:
plt.figure(figsize=(17,5))
my_levels = [0.01, 0.05, 0.09, 0.13]
my_colors = ['white', 'yellow', 'cyan', 'pink']
rh.plot.pcolormesh(x='time', y='z')
cs = (my_dataset.Specific_cloud_ice_water_content*1000.0).max(dim=('y','x')).plot.contour(x='time', y='z',
levels=my_levels,
colors=my_colors)
plt.clabel(cs, inline=1, fontsize=10, fmt='%0.3f')
plt.ylim([0,12000])
Out[61]:
In [493]:
def save_one_ecmwf_clouds(dataset, gen_datetime):
start_str = gen_datetime.strftime('%Y%m%d_%H%M')
s3name = 'ecmwf_time_height_clwc_ciwc'
plt.figure(figsize=(17,5))
my_levels = [0.01, 0.05, 0.09, 0.13]
my_colors = ['white', 'yellow', 'cyan', 'pink']
(my_dataset.Specific_cloud_liquid_water_content*1000.0).mean(dim=('y','x')).plot.pcolormesh(x='time', y='z')
cs = (my_dataset.Specific_cloud_ice_water_content*1000.0).max(dim=('y','x')).plot.contour(x='time', y='z',
levels=my_levels,
colors=my_colors)
plt.clabel(cs, inline=1, fontsize=10, fmt='%0.3f')
plt.ylim([0,12000])
str1 = start_str + ' ECMWF Domain max cloud ice and domain mean cloud liquid water content (g/kg) \n'
str2 = 'ACE-ENA forecast guidence. ARM Climate Research Facility. scollis@anl.gov'
plt.title(str1+str2)
local_fig = tempfile.NamedTemporaryFile(suffix='.png')
fn = local_fig.name
plt.savefig(fn)
s3_key = s3name + '/' + ena_tools.gen_s3_key(gen_datetime, s3name)
s3 = boto3.resource('s3')
data = open(fn, 'rb')
s3.Bucket('aceena').put_object(Key=s3_key, Body=data, ACL='public-read')
data.close()
data2 = open(fn, 'rb')
s3.Bucket('aceena').put_object(Key='latest_' + s3name + '.png',
Body=data2, ACL='public-read')
data2.close()
In [494]:
save_one_ecmwf_clouds(my_dataset, my_these_run_times[0])
In [501]:
sstting = '/lcrc/group/earthscience/ecmwf/%Y%m%d/'
fst = 'ecmwf_%Y%m%d_%H%M.nc'
local_dir = my_these_run_times[0].strftime(sstting)
local_file = my_these_run_times[0].strftime(fst)
print(local_file)
if not os.path.exists(local_dir):
os.mkdirs(local_dir)
my_dataset.to_netcdf(local_dir+)
In [409]:
my_dataset.Specific_cloud_ice_water_content.max()
Out[409]:
In [318]:
def plot_xarray_slice(item, level, timestep, vmin, vmax, units):
plt.figure(figsize=(17,10))
ter_lat = 38.7216
ter_lon = -27.2206
gra_lat = 39.0525
gra_lon = -28.0069
ax = plt.axes(projection=ccrs.PlateCarree())
cf = item[timestep][level].plot.pcolormesh(ax=ax, transform=ccrs.PlateCarree(),
x='lon', y='lat', add_colorbar=False, vmin = vmin, vmax = vmax)
ax.set_xticks([-23, -24, -26, -28 ], crs=ccrs.PlateCarree())
ax.set_yticks([37,39,41], crs=ccrs.PlateCarree())
lon_formatter = LongitudeFormatter(zero_direction_label=True)
lat_formatter = LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
ax.plot([ter_lon, gra_lon], [ter_lat, gra_lat],
'ro', transform=ccrs.PlateCarree())
ax.text(ter_lon+.2, ter_lat+.2,
'Tericia', transform=ccrs.PlateCarree(), fontsize = 16)
ax.text(gra_lon+.2, gra_lat+.2,
'Graciosa', transform=ccrs.PlateCarree(), fontsize = 16)
coast = cfeature.NaturalEarthFeature(category='physical', scale='10m',
facecolor='none', name='coastline')
_ = ax.add_feature(coast, edgecolor='black')
plt.colorbar(cf, ax=ax, fraction=0.032, label = units)
In [438]:
plot_xarray_slice(my_dataset.Specific_cloud_ice_water_content,
100, 55, 0., 0.00002, 'CLWC (kg/kg)')
In [334]:
for i in range(30):
filename = '/data/aceena_x{0:02d}.png'.format(i)
print(filename)
plot_xarray_slice(my_dataset.Specific_cloud_liquid_water_content,
120, i, 0.0, 0.5/1000.0, 'CLWC (kg/kg)')
plt.savefig(filename)
In [81]:
In [ ]:
#my_bundle, my_these_valid_times, my_these_run_times, my_lats, my_lons
In [17]:
import smtplib
# Import the email modules we'll need
from email.mime.text import MIMEText
In [19]:
import smtplib
# Import the email modules we'll need
from email.mime.text import MIMEText
msg = MIMEText('New ECMWF data is in.')
me = 'scollis@blogin5.lcrc.anl.gov'
you = 'scollis@anl.gov'
msg['Subject'] = 'ECMWF file status'
msg['From'] = me
msg['To'] = you
# Send the message via our own SMTP server.
s = smtplib.SMTP('localhost')
s.send_message(msg)
s.quit()
In [ ]: