In [1]:
import datetime as dt

Extracting Gridded Forecasts Using Siphon


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]


/Users/charlie/dato-env/lib/python2.7/site-packages/siphon/catalog.py:62: UserWarning: URL http://nomads.ncdc.noaa.gov/thredds/catalog/gfs-004/201601/20160101/catalog.html?dataset=gfs-004/201601/20160101 returned HTML. Changing to: http://nomads.ncdc.noaa.gov/thredds/catalog/gfs-004/201601/20160101/catalog.xml?dataset=gfs-004/201601/20160101
  new_url))
Out[5]:
['GFS Grid 4  2016-01-01 18:00 UTC fct:120.grb2',
 'GFS Grid 4  2016-01-01 00:00 UTC fct:120.grb2',
 'GFS Grid 4  2016-01-01 18:00 UTC fct:183.grb2',
 'GFS Grid 4  2016-01-01 00:00 UTC fct:171.grb2',
 'GFS Grid 4  2016-01-01 12:00 UTC fct:201.grb2',
 'GFS Grid 4  2016-01-01 18:00 UTC fct:102.grb2',
 'GFS Grid 4  2016-01-01 18:00 UTC fct:219.grb2',
 'GFS Grid 4  2016-01-01 00:00 UTC fct:060.grb2',
 'GFS Grid 4  2016-01-01 00:00 UTC fct:012.grb2',
 'GFS Grid 4  2016-01-01 12:00 UTC fct:300.grb2']

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]:
[u'U-component_of_wind_height_above_ground',
 u'time',
 u'height_above_ground5',
 u'lat',
 u'lon',
 u'V-component_of_wind_height_above_ground']

In [9]:
data


Out[9]:
<type 'netCDF4._netCDF4.Dataset'>
root group (NETCDF3_CLASSIC data model, file format NETCDF3):
    Conventions: CF-1.0
    Originating_center: US National Weather Service - NCEP(WMC) (7)
    Generating_Model: Global Forecast System Model (formerly known as the Aviation)
    Product_Status: Operational products
    Product_Type: Forecast products
    title: US National Weather Service - NCEP(WMC) Global Forecast System Model (formerly known as the Aviation) Forecast products
    institution: Center US National Weather Service - NCEP(WMC) (7)
    source: Type: Forecast products Status: Operational products
    history: Direct read of GRIB-2 into NetCDF-Java 4 API
    CF:feature_type: GRID
    file_format: GRIB-2
    location: /nomads3_data/raid2/noaaport/merged/gfs4/201601/20160101/gfs_4_20160101_0000_024.grb2
    _CoordinateModelRunDate: 2016-01-01T00:00:00Z
    History: Translated to CF-1.0 Conventions by Netcdf-Java CDM (NetcdfCFWriter)
Original Dataset = /nomads3_data/raid2/noaaport/merged/gfs4/201601/20160101/gfs_4_20160101_0000_024.grb2; Translation Date = Thu Jan 14 23:28:50 EST 2016
    dimensions(sizes): time(1), height_above_ground5(3), lat(61), lon(121)
    variables(dimensions): float32 U-component_of_wind_height_above_ground(time,height_above_ground5,lat,lon), int32 time(time), float64 height_above_ground5(height_above_ground5), float64 lat(lat), float64 lon(lon), float32 V-component_of_wind_height_above_ground(time,height_above_ground5,lat,lon)
    groups: 

In [10]:
# dimension of dataset
data['U-component_of_wind_height_above_ground'][:].squeeze().shape


Out[10]:
(3, 61, 121)

Plotting Using Cartopy


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()


/Users/charlie/dato-env/lib/python2.7/site-packages/matplotlib/collections.py:590: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison
  if self._edgecolors == str('face'):

Data Manipulation Using Xray & Pandas


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]:
<xray.Dataset>
Dimensions:                                  (height_above_ground5: 3, lat: 61, lon: 121, time: 1)
Coordinates:
  * time                                     (time) datetime64[ns] 2016-01-02
  * height_above_ground5                     (height_above_ground5) float64 10.0 ...
  * lat                                      (lat) float64 50.0 49.5 49.0 ...
  * lon                                      (lon) float64 230.0 230.5 231.0 ...
Data variables:
    U-component_of_wind_height_above_ground  (time, height_above_ground5, lat, lon) float32 -9.47 ...
    V-component_of_wind_height_above_ground  (time, height_above_ground5, lat, lon) float32 12.09 ...
Attributes:
    Conventions: CF-1.0
    Originating_center: US National Weather Service - NCEP(WMC) (7)
    Generating_Model: Global Forecast System Model (formerly known as the Aviation)
    Product_Status: Operational products
    Product_Type: Forecast products
    title: US National Weather Service - NCEP(WMC) Global Forecast System Model (formerly known as the Aviation) Forecast products
    institution: Center US National Weather Service - NCEP(WMC) (7)
    source: Type: Forecast products Status: Operational products
    history: Direct read of GRIB-2 into NetCDF-Java 4 API
    CF:feature_type: GRID
    file_format: GRIB-2
    location: /nomads3_data/raid2/noaaport/merged/gfs4/201601/20160101/gfs_4_20160101_0000_024.grb2
    _CoordinateModelRunDate: 2016-01-01T00:00:00Z
    History: Translated to CF-1.0 Conventions by Netcdf-Java CDM (NetcdfCFWriter)
Original Dataset = /nomads3_data/raid2/noaaport/merged/gfs4/201601/20160101/gfs_4_20160101_0000_024.grb2; Translation Date = Thu Jan 14 23:28:55 EST 2016

In [17]:
ndf['U-component_of_wind_height_above_ground'].mean(['lat', 'lon'])


Out[17]:
<xray.DataArray 'U-component_of_wind_height_above_ground' (time: 1, height_above_ground5: 3)>
array([[ 0.5291627 ,  1.28073144,  1.39978862]], dtype=float32)
Coordinates:
  * height_above_ground5  (height_above_ground5) float64 10.0 80.0 100.0
  * time                  (time) datetime64[ns] 2016-01-02

In [18]:
ndf['U-component_of_wind_height_above_ground'].max(['lat', 'lon'])


Out[18]:
<xray.DataArray 'U-component_of_wind_height_above_ground' (time: 1, height_above_ground5: 3)>
array([[ 14.23999977,  15.88000011,  16.03000069]], dtype=float32)
Coordinates:
  * height_above_ground5  (height_above_ground5) float64 10.0 80.0 100.0
  * time                  (time) datetime64[ns] 2016-01-02

In [19]:
monthly_ncdf = xray.open_mfdataset('./201512/*fct:024.grb2', concat_dim='time')

In [20]:
monthly_ncdf


Out[20]:
<xray.Dataset>
Dimensions:                                  (height_above_ground: 3, height_above_ground2: 1, height_above_ground3: 1, height_above_ground5: 3, lat: 2, lon: 2, pressure: 26, time: 74)
Coordinates:
  * height_above_ground2                     (height_above_ground2) float64 2.0
  * lat                                      (lat) float64 48.0 47.5
  * lon                                      (lon) float64 237.5 238.0
  * height_above_ground3                     (height_above_ground3) float64 80.0
  * height_above_ground5                     (height_above_ground5) float64 10.0 ...
  * pressure                                 (pressure) float64 1e+05 ...
  * height_above_ground                      (height_above_ground) float64 2.0 ...
  * time                                     (time) datetime64[ns] 2015-12-02 ...
Data variables:
    Dew_point_temperature                    (time, height_above_ground2, lat, lon) float32 278.3 ...
    Pressure_height_above_ground             (time, height_above_ground3, lat, lon) float32 100322.0 ...
    Temperature_height_above_ground          (time, height_above_ground, lat, lon) float32 282.33 ...
    Relative_humidity                        (time, pressure, lat, lon) float32 67.3 ...
    U-component_of_wind_height_above_ground  (time, height_above_ground5, lat, lon) float32 -3.37 ...
    V-component_of_wind_height_above_ground  (time, height_above_ground5, lat, lon) float32 5.55 ...
Attributes:
    Conventions: CF-1.0
    Originating_center: US National Weather Service - NCEP(WMC) (7)
    Generating_Model: Global Forecast System Model (formerly known as the Aviation)
    Product_Status: Operational products
    Product_Type: Forecast products
    title: US National Weather Service - NCEP(WMC) Global Forecast System Model (formerly known as the Aviation) Forecast products
    institution: Center US National Weather Service - NCEP(WMC) (7)
    source: Type: Forecast products Status: Operational products
    history: Direct read of GRIB-2 into NetCDF-Java 4 API
    CF:feature_type: GRID
    file_format: GRIB-2
    location: /nomads3_data/raid2/noaaport/merged/gfs4/201512/20151201/gfs_4_20151201_0000_024.grb2
    _CoordinateModelRunDate: 2015-12-01T00:00:00Z
    History: Translated to CF-1.0 Conventions by Netcdf-Java CDM (NetcdfCFWriter)
Original Dataset = /nomads3_data/raid2/noaaport/merged/gfs4/201512/20151201/gfs_4_20151201_0000_024.grb2; Translation Date = Sun Jan 03 13:46:15 EST 2016

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]:
U-component_of_wind_height_above_ground
time
2015-12-02 00:00:00 -3.175833
2015-12-02 06:00:00 -2.867500
2015-12-02 12:00:00 -2.893333
2015-12-02 18:00:00 -2.309167
2015-12-03 00:00:00 -3.525000
2015-12-03 06:00:00 -5.040833
2015-12-03 12:00:00 -3.519167
2015-12-03 18:00:00 -5.461667
2015-12-05 06:00:00 -3.254167
2015-12-05 12:00:00 -3.729167
2015-12-05 18:00:00 -4.625000
2015-12-06 00:00:00 -4.208333
2015-12-06 06:00:00 -4.567500
2015-12-06 12:00:00 -3.484167
2015-12-06 18:00:00 -0.488333
2015-12-07 00:00:00 -2.208333
2015-12-07 06:00:00 -1.858333
2015-12-07 12:00:00 -0.632500
2015-12-07 18:00:00 -1.279167
2015-12-08 00:00:00 -2.269167
2015-12-08 06:00:00 0.775000
2015-12-08 12:00:00 -4.192500
2015-12-08 18:00:00 5.130833
2015-12-09 00:00:00 3.809167
2015-12-09 06:00:00 -3.902500
2015-12-09 12:00:00 7.096667
2015-12-09 18:00:00 6.219167
2015-12-10 00:00:00 0.794167
2015-12-10 06:00:00 -5.673333
2015-12-10 12:00:00 -3.980833
... ...
2015-12-14 06:00:00 -0.556667
2015-12-14 12:00:00 -0.640000
2015-12-14 18:00:00 -0.025000
2015-12-15 00:00:00 0.075000
2015-12-15 06:00:00 -0.599167
2015-12-15 12:00:00 -0.319167
2015-12-15 18:00:00 -0.738333
2015-12-16 00:00:00 -1.080833
2015-12-16 06:00:00 -1.035000
2015-12-16 12:00:00 0.504167
2015-12-16 18:00:00 -0.188333
2015-12-17 00:00:00 -2.236667
2015-12-17 06:00:00 -3.033333
2015-12-17 12:00:00 -3.712500
2015-12-17 18:00:00 -3.839167
2015-12-18 00:00:00 -3.900000
2015-12-18 06:00:00 -4.537500
2015-12-18 18:00:00 1.812500
2015-12-22 00:00:00 2.363333
2015-12-22 18:00:00 -1.485000
2015-12-23 00:00:00 -2.120000
2015-12-29 00:00:00 0.152500
2015-12-30 12:00:00 -1.935000
2015-12-30 18:00:00 -1.610833
2015-12-31 00:00:00 -0.357500
2015-12-31 06:00:00 -2.195833
2015-12-31 12:00:00 -2.481667
2015-12-31 18:00:00 -2.305000
2016-01-01 00:00:00 -1.681667
2016-01-01 06:00:00 -2.665000

74 rows × 1 columns


In [67]:
monthly_ncdf.to_dataframe().head(10)


Out[67]:
Dew_point_temperature Pressure_height_above_ground Temperature_height_above_ground Relative_humidity U-component_of_wind_height_above_ground V-component_of_wind_height_above_ground
height_above_ground height_above_ground2 height_above_ground3 height_above_ground5 lat lon pressure time
2 2 80 10 48 237.5 100000 2015-12-02 00:00:00 278.299988 100321.531250 282.329987 67.300003 -3.37 5.55
2015-12-02 06:00:00 278.200012 100310.359375 279.510010 88.400002 -3.33 6.95
2015-12-02 12:00:00 277.899994 100422.828125 279.269989 87.599998 -3.19 5.16
2015-12-02 18:00:00 279.299988 100627.937500 281.140015 85.699997 -3.10 4.57
2015-12-03 00:00:00 280.399994 100252.875000 282.500000 79.800003 -2.75 4.24
2015-12-03 06:00:00 277.799988 99659.406250 281.369995 72.400002 -3.03 4.10
2015-12-03 12:00:00 280.000000 99015.921875 281.230011 89.000000 -2.82 5.58
2015-12-03 18:00:00 279.799988 98530.625000 281.619995 84.099998 -4.52 5.22
2015-12-05 06:00:00 278.399994 100684.500000 280.589996 80.900002 -2.96 6.75
2015-12-05 12:00:00 276.899994 100549.726562 280.429993 69.300003 -3.34 5.03

In [ ]: