In [1]:
import datetime as dt
In [2]:
from siphon.catalog import TDSCatalog
from siphon.ncss import NCSS
In [3]:
# THREDDS base url for NCEP GFS 0.5deg operational model
baseurl = 'http://nomads.ncdc.noaa.gov/thredds/catalog/gfs-004/%s/%s/catalog.html?dataset=gfs-004/%s/%s'
In [4]:
# date of model initialization time
date = dt.datetime(2016,1,1)
# reformat date to conform with THREDDS format
month_format = dt.datetime.strftime(date, '%Y%m')
day_format = dt.datetime.strftime(date, '%Y%m%d')
thredds_url = baseurl % (month_format, day_format, month_format, day_format)
In [5]:
# get list of GFS 0.5deg grib files available for given date
gfs_gribs = TDSCatalog(thredds_url)
list(gfs_gribs.datasets)[:10]
Out[5]:
In [6]:
# specifiy model init time given date, and the fcst_hour (forecast horizon) of interest
init_time = '00'
fcst_hour = '024'
# generate actual grib filename
grid_name = 'GFS Grid 4 %s %s:00 UTC fct:%s.grb2' % (dt.datetime.strftime(date, '%Y-%m-%d'), init_time, fcst_hour)
# get grib file of interest access points
grid = gfs_gribs.datasets[grid_name]
In [7]:
# use siphon's ncss class to access grid via THREDD's netcdf server
ncss = NCSS(grid.access_urls['NetcdfServer'])
# create a query of interest with boundary box, include lat lon, and query variables of interest
query = ncss.query()
query.lonlat_box(-130, -70, 20, 50)
query.add_lonlat()
query.variables( 'U-component_of_wind_height_above_ground',
'V-component_of_wind_height_above_ground', )
# query the data
data = ncss.get_data(query)
In [8]:
data.variables.keys()
Out[8]:
In [9]:
data
Out[9]:
In [10]:
# dimension of dataset
data['U-component_of_wind_height_above_ground'][:].squeeze().shape
Out[10]:
In [11]:
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import numpy as np
In [12]:
# parsing out data of interest from netcdf format into numpy arrays
y = data.variables['lat'][:].squeeze()
x = data.variables['lon'][:].squeeze()
u = data['U-component_of_wind_height_above_ground'][:].squeeze()[0]
v = data['V-component_of_wind_height_above_ground'][:].squeeze()[0]
In [13]:
# gfs lon are degrees from 0, so need to normalize
x = x - 360
In [14]:
%matplotlib inline
plt.figure(figsize=(40,20))
ax = plt.subplot(1, 1, 1, projection=ccrs.PlateCarree())
ax.coastlines('10m')
ax.set_extent([-130, -120, 45, 50], ccrs.PlateCarree())
ax.quiver(x, y, u, v, regrid_shape=20, **{'scale':300})
plt.show()
In [15]:
dataraw = ncss.get_data_raw(query)
with open('tmp_grib', 'wb') as outf:
outf.write(dataraw)
In [16]:
import xray
ndf = xray.open_dataset('./tmp_grib')
ndf
Out[16]:
In [17]:
ndf['U-component_of_wind_height_above_ground'].mean(['lat', 'lon'])
Out[17]:
In [18]:
ndf['U-component_of_wind_height_above_ground'].max(['lat', 'lon'])
Out[18]:
In [19]:
monthly_ncdf = xray.open_mfdataset('./201512/*fct:024.grb2', concat_dim='time')
In [20]:
monthly_ncdf
Out[20]:
In [21]:
thtmeans = monthly_ncdf['U-component_of_wind_height_above_ground'].mean(['lat', 'lon', 'height_above_ground5'])
In [22]:
thtmeans.to_dataframe()
Out[22]:
In [67]:
monthly_ncdf.to_dataframe().head(10)
Out[67]:
In [ ]: