test_ERAInterim_Interpolated.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using matplotlib backend: Qt4Agg\n"
     ]
    }
   ],
   "source": [
    "%run basics\n",
    "%matplotlib\n",
    "import pysolar\n",
    "import pytz\n",
    "from scipy.interpolate import InterpolatedUnivariateSpline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "do_plot = True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "site_name = \"HowardSprings\"\n",
    "tower_name = \"../Sites/\"+site_name+\"/Data/Portal/\"+site_name+\"_2014_L3.nc\"\n",
    "ds_tower = qcio.nc_read_series(tower_name)\n",
    "site_timezone = ds_tower.globalattributes[\"time_zone\"]\n",
    "site_latitude = float(ds_tower.globalattributes[\"latitude\"])\n",
    "site_longitude = float(ds_tower.globalattributes[\"longitude\"])\n",
    "tower_timestep = int(ds_tower.globalattributes[\"time_step\"])\n",
    "dt_tower = ds_tower.series[\"DateTime\"][\"Data\"]\n",
    "Fsd_tower,f,a = qcutils.GetSeriesasMA(ds_tower,\"Fsd\")\n",
    "Fsu_tower,f,a = qcutils.GetSeriesasMA(ds_tower,\"Fsu\")\n",
    "Fn_sw_tower = Fsd_tower - Fsu_tower\n",
    "Fld_tower,f,a = qcutils.GetSeriesasMA(ds_tower,\"Fld\")\n",
    "Flu_tower,f,a = qcutils.GetSeriesasMA(ds_tower,\"Flu\")\n",
    "Fn_lw_tower = Fld_tower - Flu_tower\n",
    "Fn_tower,f,a = qcutils.GetSeriesasMA(ds_tower,\"Fn\")\n",
    "Fh_tower,f,a = qcutils.GetSeriesasMA(ds_tower,\"Fh\")\n",
    "Fe_tower,f,a = qcutils.GetSeriesasMA(ds_tower,\"Fe\")\n",
    "Fg_tower,f,a = qcutils.GetSeriesasMA(ds_tower,\"Fg\")\n",
    "Fa_tower,f,a = qcutils.GetSeriesasMA(ds_tower,\"Fa\")\n",
    "ps_tower,f,a = qcutils.GetSeriesasMA(ds_tower,\"ps\")\n",
    "Ta_tower,f,a = qcutils.GetSeriesasMA(ds_tower,\"Ta\")\n",
    "Ah_tower,f,a = qcutils.GetSeriesasMA(ds_tower,\"Ah\")\n",
    "Precip_tower,f,a = qcutils.GetSeriesasMA(ds_tower,\"Precip\")\n",
    "Sws_tower,f,a = qcutils.GetSeriesasMA(ds_tower,\"Sws\")\n",
    "Ts_tower,f,a = qcutils.GetSeriesasMA(ds_tower,\"Ts\")\n",
    "Ws_tower,f,a = qcutils.GetSeriesasMA(ds_tower,\"Ws\")\n",
    "Wd_tower,f,a = qcutils.GetSeriesasMA(ds_tower,\"Wd\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "erai_name = \"../ECMWF/ECMWF_2014.nc\"\n",
    "erai_file = netCDF4.Dataset(erai_name)\n",
    "latitude = erai_file.variables[\"latitude\"][:]\n",
    "longitude = erai_file.variables[\"longitude\"][:]\n",
    "lat_resolution = abs(latitude[-1]-latitude[0])/(len(latitude)-1)\n",
    "lon_resolution = abs(longitude[-1]-longitude[0])/(len(longitude)-1)\n",
    "# index of the site in latitude dimension\n",
    "site_lat_index = int(((latitude[0]-site_latitude)/lat_resolution)+0.5)\n",
    "erai_latitude = latitude[site_lat_index]\n",
    "# index of the site in longitude dimension\n",
    "site_lon_index = int(((site_longitude-longitude[0])/lon_resolution)+0.5)\n",
    "erai_longitude = longitude[site_lon_index]\n",
    "erai_timestep = 180"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# get the time and convert to Python datetime object\n",
    "erai_time = erai_file.variables[\"time\"][:]\n",
    "time_units = getattr(erai_file.variables[\"time\"],\"units\")\n",
    "dt_erai =  netCDF4.num2date(erai_time,time_units)\n",
    "hour_utc = numpy.array([dt.hour for dt in dt_erai])\n",
    "# get the UTC and local datetime series\n",
    "site_tz = pytz.timezone(site_timezone)\n",
    "# make utc_dt timezone aware\n",
    "dt_erai_utc = [x.replace(tzinfo=pytz.utc) for x in dt_erai]\n",
    "# get local time from UTC\n",
    "dt_erai_loc = [x.astimezone(site_tz) for x in dt_erai_utc]\n",
    "#dt_erai_loc = [x.astimezone(site_tz) for x in dt_erai_utc]\n",
    "# NOTE: will have to disable daylight saving at some stage, towers stay on Standard Time\n",
    "# PRI hopes that the following line will do this ...\n",
    "#dt_ecmwf_loc = [x-x.dst() for x in dt_ecmwf_loc]\n",
    "# make local time timezone naive to match datetimes in OzFluxQC\n",
    "dt_erai_loc_ntz = [x.replace(tzinfo=None) for x in dt_erai_loc]\n",
    "dt_erai_utc_ntz = [x.replace(tzinfo=None) for x in dt_erai_utc]\n",
    "# get the datetime in the middle of the accumulation period\n",
    "erai_offset = datetime.timedelta(minutes=float(erai_timestep)/2)\n",
    "# subtract this from the ERA-I datetime to align ERA-I 3 hourly points with the middle\n",
    "# of the accumulation period\n",
    "dt_erai_loc_cor = [x - erai_offset for x in dt_erai_loc_ntz]\n",
    "dt_erai_utc_cor = [x - erai_offset for x in dt_erai_utc_ntz]\n",
    "# now we get the datetime series at the tower time step\n",
    "tdts = datetime.timedelta(minutes=tower_timestep)\n",
    "# local time for plotting\n",
    "start_date = dt_erai_loc_cor[0]\n",
    "end_date = dt_erai_loc_cor[-1]\n",
    "dt_erai_loc_tts = [x for x in qcutils.perdelta(start_date,end_date,tdts)]\n",
    "# UTC for calculating solar altitude\n",
    "start_date = dt_erai_utc_cor[0]\n",
    "end_date = dt_erai_utc_cor[-1]\n",
    "dt_erai_utc_tts = [x for x in qcutils.perdelta(start_date,end_date,tdts)]\n",
    "# UTC at tower time step as number for interpolation\n",
    "erai_time_tts = netCDF4.date2num(dt_erai_utc_tts,time_units)\n",
    "erai_time_3hr = netCDF4.date2num(dt_erai_utc_cor,time_units)\n",
    "# get the solar altitude, we will use this later to interpolate the ERA Interim\n",
    "# data from the ERA-I 3 hour time step to the tower time step.\n",
    "# NOTE: alt_solar is in degrees\n",
    "alt_solar_3hr = numpy.array([pysolar.GetAltitude(erai_latitude,erai_longitude,dt) for dt in dt_erai_utc_cor])\n",
    "# get the solar altitude at the tower time step\n",
    "alt_solar_tts = numpy.array([pysolar.GetAltitude(erai_latitude,erai_longitude,dt) for dt in dt_erai_utc_tts])\n",
    "idx = numpy.where(alt_solar_tts<=0)[0]\n",
    "alt_solar_tts[idx] = float(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Interpolate the 3 hourly accumulated downwelling shortwave to the tower time step\n",
    "# NOTE: ERA-I variables are dimensioned [time,latitude,longitude]\n",
    "Fsd_3d = erai_file.variables[\"ssrd\"][:,:,:]\n",
    "Fsd_accum = Fsd_3d[:,site_lat_index,site_lon_index]\n",
    "# Downwelling shortwave in ERA-I is a cummulative value that is reset to 0 at 0300 and 1500 UTC.\n",
    "# Here we convert the cummulative values to 3 hourly values.\n",
    "Fsd_erai_3hr = numpy.ediff1d(Fsd_accum,to_begin=0)\n",
    "idx = numpy.where((hour_utc==3)|(hour_utc==15))[0]\n",
    "Fsd_erai_3hr[idx] = Fsd_accum[idx]\n",
    "Fsd_erai_3hr = Fsd_erai_3hr/(erai_timestep*60)\n",
    "# normalise the ERA-I downwelling shortwave by the solar altitude\n",
    "coef_3hr = Fsd_erai_3hr/numpy.sin(numpy.deg2rad(alt_solar_3hr))\n",
    "# get the spline interpolation function\n",
    "s = InterpolatedUnivariateSpline(erai_time_3hr, coef_3hr, k=1)\n",
    "# get the coefficient at the tower time step\n",
    "coef_tts = s(erai_time_tts)\n",
    "# get the downwelling solar radiation at the tower time step\n",
    "Fsd_erai_tts = coef_tts*numpy.sin(numpy.deg2rad(alt_solar_tts))\n",
    "if do_plot==True:\n",
    "    fig = plt.figure()\n",
    "    plt.plot(dt_tower,Fsd_tower,'b-')\n",
    "    plt.plot(dt_erai_loc_cor,Fsd_erai_3hr,'r+')\n",
    "    plt.plot(dt_erai_loc_tts,Fsd_erai_tts,'g^')\n",
    "    plt.tight_layout()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Interpolate the 3 hourly accumulated net shortwave to the tower time step\n",
    "# NOTE: ERA-I variables are dimensioned [time,latitude,longitude]\n",
    "Fn_sw_3d = erai_file.variables[\"ssr\"][:,:,:]\n",
    "Fn_sw_accum = Fn_sw_3d[:,site_lat_index,site_lon_index]\n",
    "# Net shortwave in ERA-I is a cummulative value that is reset to 0 at 0300 and 1500 UTC.\n",
    "# Here we convert the cummulative values to 3 hourly values.\n",
    "Fn_sw_erai_3hr = numpy.ediff1d(Fn_sw_accum,to_begin=0)\n",
    "# deal with the reset times at 0300 and 1500\n",
    "idx = numpy.where((hour_utc==3)|(hour_utc==15))[0]\n",
    "Fn_sw_erai_3hr[idx] = Fn_sw_accum[idx]\n",
    "# get the average value over the 3 hourly period\n",
    "Fn_sw_erai_3hr = Fn_sw_erai_3hr/(erai_timestep*60)\n",
    "# normalise the ERA-I et shortwave by the solar altitude\n",
    "coef_3hr = Fn_sw_erai_3hr/numpy.sin(numpy.deg2rad(alt_solar_3hr))\n",
    "# get the spline interpolation function\n",
    "s = InterpolatedUnivariateSpline(erai_time_3hr, coef_3hr, k=1)\n",
    "# get the coefficient at the tower time step\n",
    "coef_tts = s(erai_time_tts)\n",
    "# get the downwelling solar radiation at the tower time step\n",
    "Fn_sw_erai_tts = coef_tts*numpy.sin(numpy.deg2rad(alt_solar_tts))\n",
    "if do_plot==True:\n",
    "    fig = plt.figure()\n",
    "    plt.plot(dt_tower,Fn_sw_tower,'b-')\n",
    "    plt.plot(dt_erai_loc_cor,Fn_sw_erai_3hr,'r+')\n",
    "    plt.plot(dt_erai_loc_tts,Fn_sw_erai_tts,'g^')\n",
    "    plt.tight_layout()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Fsu_erai_tts = Fsd_erai_tts - Fn_sw_erai_tts\n",
    "if do_plot==True:\n",
    "    fig = plt.figure()\n",
    "    plt.plot(dt_tower,Fsu_tower,'b-')\n",
    "    plt.plot(dt_erai_loc_tts,Fsu_erai_tts,'g^')\n",
    "    plt.tight_layout()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Interpolate the 3 hourly accumulated downwelling longwave to the tower time step\n",
    "# NOTE: ERA-I variables are dimensioned [time,latitude,longitude]\n",
    "Fld_3d = erai_file.variables[\"strd\"][:,:,:]\n",
    "Fld_accum = Fld_3d[:,site_lat_index,site_lon_index]\n",
    "# Downwelling longwave in ERA-I is a cummulative value that is reset to 0 at 0300 and 1500 UTC.\n",
    "# Here we convert the cummulative values to 3 hourly values.\n",
    "Fld_erai_3hr = numpy.ediff1d(Fld_accum,to_begin=0)\n",
    "idx = numpy.where((hour_utc==3)|(hour_utc==15))[0]\n",
    "Fld_erai_3hr[idx] = Fld_accum[idx]\n",
    "Fld_erai_3hr = Fld_erai_3hr/(erai_timestep*60)\n",
    "# get the spline interpolation function\n",
    "s = InterpolatedUnivariateSpline(erai_time_3hr, Fld_erai_3hr, k=1)\n",
    "# get the downwelling longwave at the tower time step\n",
    "Fld_erai_tts = s(erai_time_tts)\n",
    "if do_plot==True:\n",
    "    fig = plt.figure()\n",
    "    plt.plot(dt_tower,Fld_tower,'b-')\n",
    "    plt.plot(dt_erai_loc_cor,Fld_erai_3hr,'r+')\n",
    "    plt.plot(dt_erai_loc_tts,Fld_erai_tts,'g^')\n",
    "    plt.tight_layout()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Interpolate the 3 hourly accumulated net longwave to the tower time step\n",
    "# NOTE: ERA-I variables are dimensioned [time,latitude,longitude]\n",
    "Fn_lw_3d = erai_file.variables[\"str\"][:,:,:]\n",
    "Fn_lw_accum = Fn_lw_3d[:,site_lat_index,site_lon_index]\n",
    "# Net longwave in ERA-I is a cummulative value that is reset at 0300 and 1500 UTC.\n",
    "# Here we convert the cummulative values to 3 hourly values.\n",
    "Fn_lw_erai_3hr = numpy.ediff1d(Fn_lw_accum,to_begin=0)\n",
    "# deal with the reset times at 0300 and 1500\n",
    "idx = numpy.where((hour_utc==3)|(hour_utc==15))[0]\n",
    "Fn_lw_erai_3hr[idx] = Fn_lw_accum[idx]\n",
    "# get the average value over the 3 hourly period\n",
    "Fn_lw_erai_3hr = Fn_lw_erai_3hr/(erai_timestep*60)\n",
    "# get the spline interpolation function\n",
    "s = InterpolatedUnivariateSpline(erai_time_3hr, Fn_lw_erai_3hr, k=1)\n",
    "# get the net longwave at the tower time step\n",
    "Fn_lw_erai_tts = s(erai_time_tts)\n",
    "if do_plot==True:\n",
    "    fig = plt.figure()\n",
    "    plt.plot(dt_tower,Fn_lw_tower,'b-')\n",
    "    plt.plot(dt_erai_loc_cor,Fn_lw_erai_3hr,'r+')\n",
    "    plt.plot(dt_erai_loc_tts,Fn_lw_erai_tts,'g^')\n",
    "    plt.tight_layout()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Flu_erai_tts = Fld_erai_tts - Fn_lw_erai_tts\n",
    "if do_plot==True:\n",
    "    fig = plt.figure()\n",
    "    plt.plot(dt_tower,Flu_tower,'b-')\n",
    "    plt.plot(dt_erai_loc_tts,Flu_erai_tts,'g^')\n",
    "    plt.tight_layout()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Fn_erai_tts = (Fsd_erai_tts - Fsu_erai_tts) + (Fld_erai_tts - Flu_erai_tts)\n",
    "if do_plot==True:\n",
    "    fig = plt.figure()\n",
    "    plt.plot(dt_tower,Fn_tower,'b-')\n",
    "    plt.plot(dt_erai_loc_tts,Fn_erai_tts,'g^')\n",
    "    plt.tight_layout()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Interpolate the 3 hourly accumulated sensible heat flux to the tower time step\n",
    "# NOTE: ERA-I variables are dimensioned [time,latitude,longitude]\n",
    "Fh_3d = erai_file.variables[\"sshf\"][:,:,:]\n",
    "Fh_accum = float(-1)*Fh_3d[:,site_lat_index,site_lon_index]\n",
    "# Sensible heat flux in ERA-I is a cummulative value that is reset at 0300 and 1500 UTC.\n",
    "# Here we convert the cummulative values to 3 hourly values.\n",
    "Fh_erai_3hr = numpy.ediff1d(Fh_accum,to_begin=0)\n",
    "# deal with the reset times at 0300 and 1500\n",
    "idx = numpy.where((hour_utc==3)|(hour_utc==15))[0]\n",
    "Fh_erai_3hr[idx] = Fh_accum[idx]\n",
    "# get the average value over the 3 hourly period\n",
    "Fh_erai_3hr = Fh_erai_3hr/(erai_timestep*60)\n",
    "# get the spline interpolation function\n",
    "s = InterpolatedUnivariateSpline(erai_time_3hr, Fh_erai_3hr, k=1)\n",
    "# get the net longwave at the tower time step\n",
    "Fh_erai_tts = s(erai_time_tts)\n",
    "if do_plot==True:\n",
    "    fig = plt.figure()\n",
    "    plt.plot(dt_tower,Fh_tower,'b-')\n",
    "    plt.plot(dt_erai_loc_cor,Fh_erai_3hr,'r+')\n",
    "    plt.plot(dt_erai_loc_tts,Fh_erai_tts,'g^')\n",
    "    plt.tight_layout()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Interpolate the 3 hourly accumulated latent heat flux to the tower time step\n",
    "# NOTE: ERA-I variables are dimensioned [time,latitude,longitude]\n",
    "Fe_3d = erai_file.variables[\"slhf\"][:,:,:]\n",
    "Fe_accum = float(-1)*Fe_3d[:,site_lat_index,site_lon_index]\n",
    "# Latent heat flux in ERA-I is a cummulative value that is reset at 0300 and 1500 UTC.\n",
    "# Here we convert the cummulative values to 3 hourly values.\n",
    "Fe_erai_3hr = numpy.ediff1d(Fe_accum,to_begin=0)\n",
    "# deal with the reset times at 0300 and 1500\n",
    "idx = numpy.where((hour_utc==3)|(hour_utc==15))[0]\n",
    "Fe_erai_3hr[idx] = Fe_accum[idx]\n",
    "# get the average value over the 3 hourly period\n",
    "Fe_erai_3hr = Fe_erai_3hr/(erai_timestep*60)\n",
    "# get the spline interpolation function\n",
    "s = InterpolatedUnivariateSpline(erai_time_3hr, Fe_erai_3hr, k=1)\n",
    "# get the net longwave at the tower time step\n",
    "Fe_erai_tts = s(erai_time_tts)\n",
    "if do_plot==True:\n",
    "    fig = plt.figure()\n",
    "    plt.plot(dt_tower,Fe_tower,'b-')\n",
    "    plt.plot(dt_erai_loc_cor,Fe_erai_3hr,'r+')\n",
    "    plt.plot(dt_erai_loc_tts,Fe_erai_tts,'g^')\n",
    "    plt.tight_layout()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# get Fg as residual\n",
    "Fg_erai_tts = Fn_erai_tts - Fh_erai_tts - Fe_erai_tts\n",
    "# and then Fa\n",
    "Fa_erai_tts = Fn_erai_tts - Fg_erai_tts\n",
    "if do_plot==True:\n",
    "    fig = plt.figure()\n",
    "    ax1 = plt.subplot(211)\n",
    "    plt.plot(dt_tower,Fg_tower,'b-')\n",
    "    plt.plot(dt_erai_loc_tts,Fg_erai_tts,'g^')\n",
    "    ax2 = plt.subplot(212,sharex=ax1)\n",
    "    plt.plot(dt_tower,Fa_tower,'b-')\n",
    "    plt.plot(dt_erai_loc_tts,Fa_erai_tts,'g^')\n",
    "    plt.tight_layout()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Interpolate the 3 hourly air pressure to the tower time step\n",
    "# NOTE: ERA-I variables are dimensioned [time,latitude,longitude]\n",
    "ps_3d = erai_file.variables[\"sp\"][:,:,:]\n",
    "ps_erai_3hr = ps_3d[:,site_lat_index,site_lon_index]/float(1000)\n",
    "# get the spline interpolation function\n",
    "s = InterpolatedUnivariateSpline(erai_time_3hr, ps_erai_3hr, k=1)\n",
    "# get the air pressure at the tower time step\n",
    "ps_erai_tts = s(erai_time_tts)\n",
    "if do_plot==True:\n",
    "    fig = plt.figure()\n",
    "    plt.plot(dt_tower,ps_tower,'b-')\n",
    "    plt.plot(dt_erai_loc_cor,ps_erai_3hr,'r+')\n",
    "    plt.plot(dt_erai_loc_tts,ps_erai_tts,'g^')\n",
    "    plt.tight_layout()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Interpolate the 3 hourly air temperature to the tower time step\n",
    "# NOTE: ERA-I variables are dimensioned [time,latitude,longitude]\n",
    "Ta_3d = erai_file.variables[\"t2m\"][:,:,:]\n",
    "Ta_erai_3hr = Ta_3d[:,site_lat_index,site_lon_index] - 273.15\n",
    "# get the spline interpolation function\n",
    "s = InterpolatedUnivariateSpline(erai_time_3hr, Ta_erai_3hr, k=1)\n",
    "# get the air temperature at the tower time step\n",
    "Ta_erai_tts = s(erai_time_tts)\n",
    "if do_plot==True:\n",
    "    fig = plt.figure()\n",
    "    plt.plot(dt_tower,Ta_tower,'b-')\n",
    "    plt.plot(dt_erai_loc_cor,Ta_erai_3hr,'r+')\n",
    "    plt.plot(dt_erai_loc_tts,Ta_erai_tts,'g^')\n",
    "    plt.tight_layout()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Interpolate the 3 hourly dew point temperature to the tower time step\n",
    "# and convert to Ah, RH and q\n",
    "# NOTE: ERA-I variables are dimensioned [time,latitude,longitude]\n",
    "Td_3d = erai_file.variables[\"d2m\"][:,:,:]\n",
    "Td_erai_3hr = Td_3d[:,site_lat_index,site_lon_index] - 273.15\n",
    "# get the spline interpolation function\n",
    "s = InterpolatedUnivariateSpline(erai_time_3hr, Td_erai_3hr, k=1)\n",
    "# get the dew point temperature at the tower time step\n",
    "Td_erai_tts = s(erai_time_tts)\n",
    "# get the relative humidity\n",
    "es_erai_tts = mf.es(Ta_erai_tts)\n",
    "e_erai_tts = mf.es(Td_erai_tts)\n",
    "RH_erai_tts = float(100)*e_erai_tts/es_erai_tts\n",
    "# get the absolute humidity\n",
    "Ah_erai_tts = mf.absolutehumidityfromRH(Ta_erai_tts,RH_erai_tts)\n",
    "# get the specific humidity\n",
    "q_erai_tts = mf.specifichumidityfromRH(RH_erai_tts,Ta_erai_tts,ps_erai_tts)\n",
    "if do_plot==True:\n",
    "    fig = plt.figure()\n",
    "    plt.plot(dt_tower,Ah_tower,'b-')\n",
    "    plt.plot(dt_erai_loc_tts,Ah_erai_tts,'g^')\n",
    "    plt.tight_layout()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Interpolate the 3 hourly boundary layer height to the tower time step\n",
    "# NOTE: ERA-I variables are dimensioned [time,latitude,longitude]\n",
    "Habl_3d = erai_file.variables[\"blh\"][:,:,:]\n",
    "Habl_erai_3hr = Habl_3d[:,site_lat_index,site_lon_index]\n",
    "# get the spline interpolation function\n",
    "s = InterpolatedUnivariateSpline(erai_time_3hr, Habl_erai_3hr, k=1)\n",
    "# get the boundary layer height at the tower time step\n",
    "Habl_erai_tts = s(erai_time_tts)\n",
    "if do_plot==True:\n",
    "    fig = plt.figure()\n",
    "    plt.plot(dt_erai_loc_cor,Habl_erai_3hr,'r+')\n",
    "    plt.plot(dt_erai_loc_tts,Habl_erai_tts,'g^')\n",
    "    plt.tight_layout()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Spread the 3 hourly accumulated precipitation to the tower time step\n",
    "# NOTE: ERA-I variables are dimensioned [time,latitude,longitude]\n",
    "Precip_3d = erai_file.variables[\"tp\"][:,:,:]\n",
    "Precip_accum = Precip_3d[:,site_lat_index,site_lon_index]\n",
    "\n",
    "Precip_erai_3hr = numpy.ediff1d(Precip_accum,to_begin=0)\n",
    "idx = numpy.where((hour_utc==3)|(hour_utc==15))[0]\n",
    "Precip_erai_3hr[idx] = Precip_accum[idx]\n",
    "Precip_erai_3hr = Precip_erai_3hr*float(1000)\n",
    "\n",
    "Precip_erai_tts = numpy.zeros(len(dt_erai_loc_tts))\n",
    "idx = qcutils.FindIndicesOfBInA(dt_erai_loc_cor,dt_erai_loc_tts)\n",
    "Precip_erai_tts[idx] = Precip_erai_3hr\n",
    "\n",
    "if do_plot:\n",
    "    fig = plt.figure()\n",
    "    plt.plot(dt_tower,Precip_tower,'b-')\n",
    "    plt.plot(dt_erai_loc_cor,Precip_erai_3hr,'r+')\n",
    "    plt.plot(dt_erai_loc_tts,Precip_erai_tts,'g^')\n",
    "    plt.show()\n",
    "\n",
    "#idx = qcutils.FindIndicesOfBInA(dt_erai_loc_tts,dt_tower)\n",
    "\n",
    "#Precip_erai_tts_ma = numpy.ma.array(Precip_erai_tts,mask=Precip_tower.mask)\n",
    "#Precip_erai_tts_cs = numpy.ma.cumsum(Precip_erai_tts)\n",
    "#Precip_tower_cs = numpy.ma.cumsum(Precip_tower)\n",
    "\n",
    "#if do_plot:\n",
    "#    fig = plt.figure()\n",
    "#    plt.plot(dt_tower,Precip_tower_cs,'b-')\n",
    "#    plt.plot(dt_erai_loc_tts,Precip_erai_tts_cs,'g^')\n",
    "#    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Interpolate the 3 hourly soil moisture to the tower time step\n",
    "# NOTE: ERA-I variables are dimensioned [time,latitude,longitude]\n",
    "Sws_3d = erai_file.variables[\"swvl1\"][:,:,:]\n",
    "Sws_erai_3hr = Sws_3d[:,site_lat_index,site_lon_index]\n",
    "# get the spline interpolation function\n",
    "s = InterpolatedUnivariateSpline(erai_time_3hr, Sws_erai_3hr, k=1)\n",
    "# get the soil moisture at the tower time step\n",
    "Sws_erai_tts = s(erai_time_tts)\n",
    "if do_plot==True:\n",
    "    fig = plt.figure()\n",
    "    plt.plot(dt_tower,Sws_tower,'b-')\n",
    "    plt.plot(dt_erai_loc_cor,Sws_erai_3hr,'r+')\n",
    "    plt.plot(dt_erai_loc_tts,Sws_erai_tts,'g^')\n",
    "    plt.tight_layout()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Interpolate the 3 hourly soil temperature to the tower time step\n",
    "# NOTE: ERA-I variables are dimensioned [time,latitude,longitude]\n",
    "Ts_3d = erai_file.variables[\"stl1\"][:,:,:]\n",
    "Ts_erai_3hr = Ts_3d[:,site_lat_index,site_lon_index] - 273.15\n",
    "# get the spline interpolation function\n",
    "s = InterpolatedUnivariateSpline(erai_time_3hr, Ts_erai_3hr, k=1)\n",
    "# get the soil moisture at the tower time step\n",
    "Ts_erai_tts = s(erai_time_tts)\n",
    "if do_plot==True:\n",
    "    fig = plt.figure()\n",
    "    plt.plot(dt_tower,Ts_tower,'b-')\n",
    "    plt.plot(dt_erai_loc_cor,Ts_erai_3hr,'r+')\n",
    "    plt.plot(dt_erai_loc_tts,Ts_erai_tts,'g^')\n",
    "    plt.tight_layout()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Interpolate the 3 hourly U and V components to the tower time step\n",
    "# NOTE: ERA-I variables are dimensioned [time,latitude,longitude]\n",
    "# U first ...\n",
    "U_3d = erai_file.variables[\"u10\"][:,:,:]\n",
    "U_erai_3hr = U_3d[:,site_lat_index,site_lon_index]\n",
    "# get the spline interpolation function\n",
    "s = InterpolatedUnivariateSpline(erai_time_3hr, U_erai_3hr, k=1)\n",
    "# get the soil moisture at the tower time step\n",
    "U_erai_tts = s(erai_time_tts)\n",
    "# ... then V\n",
    "V_3d = erai_file.variables[\"v10\"][:,:,:]\n",
    "V_erai_3hr = V_3d[:,site_lat_index,site_lon_index]\n",
    "# get the spline interpolation function\n",
    "s = InterpolatedUnivariateSpline(erai_time_3hr, V_erai_3hr, k=1)\n",
    "# get the soil moisture at the tower time step\n",
    "V_erai_tts = s(erai_time_tts)\n",
    "# now get the wind speed and direction\n",
    "Ws_erai_tts = numpy.sqrt(U_erai_tts*U_erai_tts + V_erai_tts*V_erai_tts)\n",
    "Wd_erai_tts = float(270) - numpy.arctan2(V_erai_tts,U_erai_tts)*float(180)/numpy.pi\n",
    "idx = numpy.where(Wd_erai_tts>360)[0]\n",
    "if len(idx)>0: Wd_erai_tts[idx] = Wd_erai_tts[idx] - float(360)\n",
    "\n",
    "if do_plot==True:\n",
    "    fig = plt.figure()\n",
    "    ax1 = plt.subplot(211)\n",
    "    plt.plot(dt_tower,Ws_tower,'b-')\n",
    "    plt.plot(dt_erai_loc_tts,Ws_erai_tts,'g^')\n",
    "    ax2 = plt.subplot(212,sharex=ax1)\n",
    "    plt.plot(dt_tower,Wd_tower,'b-')\n",
    "    plt.plot(dt_erai_loc_tts,Wd_erai_tts,'g^')\n",
    "    plt.tight_layout()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}