RHLRC_NoD.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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def residuals_RHLRC_NoD(params,NEE,Fsd,weighted,weight_type):\n",
    "    alpha = params[0]\n",
    "    beta = params[1]\n",
    "    gamma = params[2]\n",
    "    resid = NEE - NEE_RHLRC_NoD(Fsd,alpha,beta,gamma)\n",
    "    if weighted:\n",
    "        weights = CalculateWeights(resid,weight_type=weight_type)\n",
    "        return resid*weights\n",
    "    else:\n",
    "        return resid\n",
    "\n",
    "def NEE_RHLRC_NoD(Fsd,alpha,beta,gamma):\n",
    "    NEE = -1*GPP_RHLRC_NoD(Fsd,alpha,beta) + gamma\n",
    "    return NEE\n",
    "\n",
    "def GPP_RHLRC_NoD(Fsd,alpha,beta):\n",
    "    GPP = alpha*beta*Fsd/(alpha*Fsd+beta)\n",
    "    return GPP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def irls_leastsq(func,p0,args=(),maxfev=3,weight_type='Huber',mode=\"quiet\"):\n",
    "    weighted = False\n",
    "    popt,pcov = scipy.optimize.leastsq(func,p0,args=args+(weighted,weight_type))\n",
    "    if mode!=\"quiet\":\n",
    "        print 'After non-weighted call to leastsq: ',popt\n",
    "    n=1\n",
    "    weighted = True\n",
    "    while n<=maxfev:\n",
    "        popt,pcov = scipy.optimize.leastsq(func,popt,args=args+(weighted,weight_type))\n",
    "        if mode!=\"quiet\":\n",
    "            print 'Weighted call '+str(n)+' to leastsq: ',popt\n",
    "        n = n + 1\n",
    "    return popt,pcov\n",
    "\n",
    "def CalculateWeights(resid,weight_type='Huber'):\n",
    "    if weight_type.lower()=='tukey':\n",
    "        weights = TukeyBiweight(MAD_scale(resid))\n",
    "    elif weight_type.lower()=='huber':\n",
    "        weights = HuberTWeight(MAD_scale(resid))\n",
    "    elif weight_type.lower()=='ols':\n",
    "        weights = OLSWeight(MAD_scale(resid))\n",
    "    else:\n",
    "        print \"CalculateWeights: Unknown weight type, used Huber's t weight\"\n",
    "        weights = HuberTWeight(MAD_scale(resid))\n",
    "    return weights\n",
    "\n",
    "def TukeyBiweight(x):\n",
    "    return ((numpy.abs(x)<1.0)*(1.0 - x**2)**2)\n",
    "\n",
    "def HuberTWeight(x):\n",
    "    weights = numpy.ones(len(x))\n",
    "    idx = numpy.where(numpy.abs(x)>1.345)[0]\n",
    "    weights[idx] = 1.345/numpy.abs(x[idx])\n",
    "    return weights\n",
    "\n",
    "def OLSWeight(x):\n",
    "    weights = numpy.ones(len(x))\n",
    "    return weights\n",
    "\n",
    "def MAD_scale(resid):\n",
    "    scale = numpy.median(numpy.abs(resid - numpy.median(resid))) / 0.6745\n",
    "    return resid/scale"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# get the netCDF file name\n",
    "ncname = \"../Sites/Whroo/Data/Portal/Whroo_L4.nc\"\n",
    "# read the contents into a data structure\n",
    "ds = qcio.nc_read_series(ncname)\n",
    "# get the datetime series\n",
    "dt = ds.series[\"DateTime\"][\"Data\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# open an Excel workbook for the results\n",
    "xlname = ncname.replace(\".nc\",\"_LRC.xls\")\n",
    "xlfile = qcio.xl_open_write(xlname)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/peter/anaconda/lib/python2.7/site-packages/matplotlib/pyplot.py:424: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).\n",
      "  max_open_warning, RuntimeWarning)\n",
      "/home/peter/anaconda/lib/python2.7/site-packages/numpy/ma/core.py:3900: UserWarning: Warning: converting a masked element to nan.\n",
      "  warnings.warn(\"Warning: converting a masked element to nan.\")\n"
     ]
    }
   ],
   "source": [
    "# rectangular hyperbolic LRC, no D function\n",
    "# set the downwelling shortwave radiation and friction velocity thresholds\n",
    "Fsd_threshold = 10\n",
    "ustar_threshold = 0.40\n",
    "# add a sheet to the Excel workbook\n",
    "xlsheet = xlfile.add_sheet(\"RHLRC_NoD\")\n",
    "# get the index of the first period at the start of the first whole month\n",
    "si = qcutils.GetDateIndex(dt,str(dt[0]),match=\"startnextmonth\")\n",
    "# get the start datetime\n",
    "start_date = dt[si]\n",
    "# create a dictionary for the results\n",
    "RHLRC_NoD = {\"DateTime\":{\"data\":[],\"units\":\"-\",\"format\":\"dd/mm/yyyy hh:mm\"},\n",
    "             \"alpha\":{\"data\":[],\"units\":\"-\",\"format\":\"0.0000\"},\n",
    "             \"beta\":{\"data\":[],\"units\":\"-\",\"format\":\"0.00\"},\n",
    "             \"gamma\":{\"data\":[],\"units\":\"-\",\"format\":\"0.00\"},\n",
    "             \"ER\":{\"data\":[],\"units\":\"-\",\"format\":\"0.00\"}}\n",
    "# loop over periods until done\n",
    "while start_date<dt[-1]:\n",
    "    # get the end datetime for this pass\n",
    "    end_date = start_date+dateutil.relativedelta.relativedelta(months=1)\n",
    "    # indices of start and end datetimes\n",
    "    si = qcutils.GetDateIndex(dt,str(start_date),default=0)\n",
    "    ei = qcutils.GetDateIndex(dt,str(end_date),default=len(dt))\n",
    "    # minimum number of points is 10%\n",
    "    min_points = int(0.1*(ei-si))\n",
    "    mi = int((si+ei)/2)\n",
    "    # get the downwelling shortwave, CO2 flux and friction velocity\n",
    "    Fsd,f,a = qcutils.GetSeriesasMA(ds,\"Fsd\",si=si,ei=ei)\n",
    "    Fc,f,a = qcutils.GetSeriesasMA(ds,\"Fc\",si=si,ei=ei)\n",
    "    ustar,f,a = qcutils.GetSeriesasMA(ds,\"ustar\",si=si,ei=ei)\n",
    "    # mask out the night-time conditions\n",
    "    Fc_day = numpy.ma.masked_where(Fsd<Fsd_threshold,Fc)\n",
    "    # impose the same mask on Fsd and Fc\n",
    "    mask = numpy.ma.mask_or(numpy.ma.getmaskarray(Fsd),\n",
    "                            numpy.ma.getmaskarray(Fc_day),\n",
    "                            copy=True)\n",
    "    # remove masked elements\n",
    "    NEE_day = numpy.ma.compressed(numpy.ma.masked_where(mask,Fc_day))\n",
    "    Fsd_day = numpy.ma.compressed(numpy.ma.masked_where(mask,Fsd))\n",
    "    # add the mid-period datetime to the results dictionary\n",
    "    RHLRC_NoD[\"DateTime\"][\"data\"].append(dt[mi])\n",
    "    # call the optimisation if more than the minimum number of\n",
    "    # points available, otherwise set results to missing value\n",
    "    if len(Fsd_day)>=min_points:\n",
    "        popt,pcov = irls_leastsq(residuals_RHLRC_NoD,[0.1,100,1],\n",
    "                                 args=(NEE_day,Fsd_day),maxfev=3,\n",
    "                                 weight_type='Huber')\n",
    "        RHLRC_NoD[\"alpha\"][\"data\"].append(popt[0])\n",
    "        RHLRC_NoD[\"beta\"][\"data\"].append(popt[1])\n",
    "        RHLRC_NoD[\"gamma\"][\"data\"].append(popt[2])\n",
    "    else:\n",
    "        RHLRC_NoD[\"alpha\"][\"data\"].append(-9999)\n",
    "        RHLRC_NoD[\"beta\"][\"data\"].append(-9999)\n",
    "        RHLRC_NoD[\"gamma\"][\"data\"].append(-9999)\n",
    "    # plot the data\n",
    "    Fsd_fit = numpy.arange(0,numpy.ma.max(Fsd_day),1)\n",
    "    NEE_fit = NEE_RHLRC_NoD(Fsd_fit,popt[0],popt[1],popt[2])\n",
    "    fig=plt.figure()\n",
    "    plt.plot(Fsd_day,NEE_day,'b.')\n",
    "    plt.plot(Fsd_fit,NEE_fit,'r-')\n",
    "    plt.show()\n",
    "    # get the average of the u*-filtered, nocturnal Fc\n",
    "    Fc_night = numpy.ma.masked_where(Fsd>=Fsd_threshold,Fc)\n",
    "    NEE_night = numpy.ma.masked_where(ustar<=ustar_threshold,Fc_night)\n",
    "    if len(NEE_night)>=min_points:\n",
    "        RHLRC_NoD[\"ER\"][\"data\"].append(numpy.ma.mean(NEE_night))\n",
    "    else:\n",
    "        RHLRC_NoD[\"ER\"][\"data\"].append(-9999)\n",
    "    start_date = end_date\n",
    "# write the results to the Excel workbook\n",
    "qcio.xl_write_data(xlsheet,RHLRC_NoD,xlCol=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# plot a time series of gamma and observed ER\n",
    "ER_modelled = numpy.ma.array(RHLRC_NoD[\"gamma\"][\"data\"])\n",
    "ER_modelled = numpy.ma.masked_where(ER_modelled==-9999,ER_modelled)\n",
    "ER_observed = numpy.ma.array(RHLRC_NoD[\"ER\"][\"data\"])\n",
    "ER_observed = numpy.ma.masked_where(ER_observed==-9999,ER_observed)\n",
    "fig = plt.figure()\n",
    "plt.plot(RHLRC_NoD[\"DateTime\"][\"data\"],ER_modelled,'bo')\n",
    "plt.plot(RHLRC_NoD[\"DateTime\"][\"data\"],ER_observed,'ro')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "xlfile.save(xlname)"
   ]
  }
 ],
 "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}