In [1]:

    

%matplotlib inline
from matplotlib import pyplot as plt


    

In [2]:

    

import os, sys
import numpy as np
from numpy import ma
import xray


    

In [3]:

    

dpath = os.path.join(os.environ['HOME'], 'data/NCEP1')


    

In [ ]:

    

dset_u = xray.open_dataset(os.path.join(dpath, 'wind/uwnd.mon.mean.nc'))
dset_v = xray.open_dataset(os.path.join(dpath, 'wind/vwnd.mon.mean.nc'))


    

In [ ]:

    

dset_u = dset_u.sel(time=slice('1948','2014'))
dset_v = dset_v.sel(time=slice('1948','2014'))


    

In [ ]:

    

lat = dset_u['lat'].data
lon = dset_u['lon'].data


    

In [ ]:

    

dates = dset_u['time'].data


    

In [ ]:

    

uwnd_1000 = dset_u.sel(level=1000)['uwnd'].data


    

In [ ]:

    

vwnd_1000 = dset_v.sel(level=1000)['vwnd'].data


    

In [ ]:

    

uwnd_850 = dset_u.sel(level=850)['uwnd'].data


    

In [ ]:

    

vwnd_850 = dset_v.sel(level=850)['vwnd'].data


    

In [ ]:

    

uwnd_200 = dset_u.sel(level=200)['uwnd'].data


    

In [ ]:

    

vwnd_200 = dset_v.sel(level=200)['vwnd'].data


    

In [ ]:

    

d = {}


    

In [ ]:

    

d['time'] = ('time',dates)
d['latitudes'] = ('latitudes',lat)
d['longitudes'] = ('longitudes', lon)
d['uwnd_1000'] = (['time','latitudes','longitudes'], uwnd_1000)
d['vwnd_1000'] = (['time','latitudes','longitudes'], vwnd_1000)
d['uwnd_850'] = (['time','latitudes','longitudes'], uwnd_850)
d['vwnd_850'] = (['time','latitudes','longitudes'], vwnd_850)
d['uwnd_200'] = (['time','latitudes','longitudes'], uwnd_200)
d['vwnd_200'] = (['time','latitudes','longitudes'], vwnd_200)


    

In [ ]:

    

wind = xray.Dataset(d)


    

In [ ]:

    

wind


    

In [ ]:

    

wind.to_netcdf('./outputs/NCEP1_monthly_wind_1948_2014.nc')


    

In [ ]:

    

dset_u.close()
dset_v.close()


    

In [4]:

    

dset_hgt = xray.open_dataset(os.path.join(dpath, 'hgt/hgt.mon.mean.nc'))


    

In [5]:

    

dset_hgt


    

    
Out[5]:





<xray.Dataset>
Dimensions:  (lat: 73, level: 17, lon: 144, time: 812)
Coordinates:
  * level    (level) float32 1000.0 925.0 850.0 700.0 600.0 500.0 400.0 ...
  * lat      (lat) float32 90.0 87.5 85.0 82.5 80.0 77.5 75.0 72.5 70.0 67.5 ...
  * lon      (lon) float32 0.0 2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0 22.5 ...
  * time     (time) datetime64[ns] 1948-01-01 1948-02-01 1948-03-01 ...
Data variables:
    hgt      (time, level, lat, lon) float64 110.0 110.0 110.0 110.0 110.0 ...
Attributes:
    description:  Data from NCEP initialized reanalysis (4x/day).  These are interpolated to pressure surfaces from model (sigma) surfaces.
    platform: Model
    Conventions: COARDS
    references: http://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis.derived.html
    NCO: 20121012
    history: Created by NOAA-CIRES Climate Diagnostics Center (SAC) from the NCEP
reanalysis data set on 07/07/97 by calc.mon.mean.year.f using
/Datasets/nmc.reanalysis.derived/pressure/hgt.mon.mean.nc
from /Datasets/nmc.reanalysis/pressure/hgt.79.nc to hgt.95.nc
Converted to chunked, deflated non-packed NetCDF4 2014/09
    title: monthly mean hgt from the NCEP Reanalysis

In [6]:

    

lat = dset_hgt['lat'].data
lon = dset_hgt['lon'].data


    

In [7]:

    

dset_hgt = dset_hgt.sel(time=slice('1948','2014'))


    

In [8]:

    

dates = dset_hgt['time'].data


    

In [9]:

    

hgt_1000 = dset_hgt.sel(level=1000)['hgt'].data


    

In [10]:

    

hgt_850 = dset_hgt.sel(level=850)['hgt'].data


    

In [11]:

    

hgt_200 = dset_hgt.sel(level=200)['hgt'].data


    

In [12]:

    

d = {}
d['time'] = ('time',dates)
d['latitudes'] = ('latitudes',lat)
d['longitudes'] = ('longitudes', lon)
d['hgt_1000'] = (['time','latitudes','longitudes'], hgt_1000)
d['hgt_850'] = (['time','latitudes','longitudes'], hgt_850)
d['hgt_200'] = (['time','latitudes','longitudes'], hgt_200)


    

In [13]:

    

hgt = xray.Dataset(d)


    

In [14]:

    

hgt


    

    
Out[14]:





<xray.Dataset>
Dimensions:     (latitudes: 73, longitudes: 144, time: 804)
Coordinates:
  * time        (time) datetime64[ns] 1948-01-01 1948-02-01 1948-03-01 ...
  * longitudes  (longitudes) float32 0.0 2.5 5.0 7.5 10.0 12.5 15.0 17.5 ...
  * latitudes   (latitudes) float32 90.0 87.5 85.0 82.5 80.0 77.5 75.0 72.5 ...
Data variables:
    hgt_200     (time, latitudes, longitudes) float64 1.08e+04 1.08e+04 ...
    hgt_1000    (time, latitudes, longitudes) float64 110.0 110.0 110.0 ...
    hgt_850     (time, latitudes, longitudes) float64 1.279e+03 1.279e+03 ...

In [15]:

    

hgt.to_netcdf('./outputs/NCEP1_monthly_hgt_1948_2014.nc')


    

In [17]:

    

!cp ./outputs/NCEP1_monthly_hgt_1948_2014.nc /Users/nicolasf/research/NIWA/paleo/pict/data/NCEP1_monthly_hgt_1948_2014.nc


    

In [16]:

    

dset_hgt.close()


    

In [ ]:

    

dset_omega = xray.open_dataset(os.path.join(dpath, 'omega/omega.mon.mean.nc'))


    

In [ ]:

    

dset_omega


    

In [ ]:

    

dset_omega = dset_omega.sel(level=500, time=slice('1948','2014'))


    

In [ ]:

    

lat = dset_omega['lat'].data
lon = dset_omega['lon'].data


    

In [ ]:

    

dates = dset_omega['time'].data


    

In [ ]:

    

omega_500 = dset_omega['omega'].data


    

In [ ]:

    

d = {}
d['time'] = ('time',dates)
d['latitudes'] = ('latitudes',lat)
d['longitudes'] = ('longitudes', lon)
d['omega_500'] = (['time','latitudes','longitudes'], omega_500)


    

In [ ]:

    

omega = xray.Dataset(d)


    

In [ ]:

    

omega.to_netcdf('./outputs/NCEP1_monthly_omega_1948_2014.nc')


    

In [ ]:

    

omega.close()


    

In [ ]:

    

omega


    

In [ ]:

    




    

In [ ]: