erai2nc.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": [
    "cf = qcio.load_controlfile(path=\"controlfiles/\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "erai_name = os.path.join(cf[\"Files\"][\"erai_path\"],cf[\"Files\"][\"erai_file\"])\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",
    "# 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",
    "erai_timestep = 180"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# get a list of sites\n",
    "site_list = cf[\"Sites\"].keys()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# now loop over the sies\n",
    "for site in site_list:\n",
    "    # get the metadata from the control file\n",
    "    out_filename = os.path.join(cf[\"Sites\"][site][\"out_filepath\"],\n",
    "                               cf[\"Sites\"][site][\"out_filename\"])\n",
    "    site_name = cf[\"Sites\"][site][\"site_name\"]\n",
    "    site_timezone = cf[\"Sites\"][site][\"site_timezone\"]\n",
    "    site_latitude = cf[\"Sites\"][site][\"site_latitude\"]\n",
    "    site_longitude = cf[\"Sites\"][site][\"site_longitude\"]\n",
    "    site_timestep = cf[\"Sites\"][site][\"site_timestep\"]\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",
    "\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)\n",
    "    \n",
    "    # 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"
   ]
  }
 ],
 "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
}