RDA_frequency_autocorrelation.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#quick script to find the autocorrelation in time of front frequency for each 1 degree latitude bin\n",
    "#without accounting for autocorrelation, changes in frequency will be found to be more significant\n",
    "#than they actually are\n",
    "\n",
    "#a bunch of this code is copied over from RDA_hov_periods.nc and slightly adapted."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import time\n",
    "import os\n",
    "import netCDF4 as nc\n",
    "import datetime\n",
    "import matplotlib.pyplot as plt\n",
    "from autocorr import autocorr\n",
    "%matplotlib notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#Access NetCDF files (currently hosted in Ferret directory)\n",
    "RDA_path_1 = \"/Users/Siwen/Desktop/ferret/bin/meiyu_clean.nc\"\n",
    "RDA_path_2 = \"/Users/Siwen/Desktop/ferret/bin/meiyu_2_clean.nc\"\n",
    "RDA_1 = nc.Dataset(RDA_path_1, 'r')\n",
    "RDA_2 = nc.Dataset(RDA_path_2, 'r')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#load data from NetCDF files to notebook\n",
    "lat_1_all =  RDA_1.variables['lat_115'][:]\n",
    "lat_2_all =  RDA_2.variables['lat_115'][:]\n",
    "intensity_1_all = RDA_1.variables['intensity'][:] \n",
    "intensity_2_all = RDA_2.variables['intensity'][:]\n",
    "\n",
    "RDA_1.close()\n",
    "RDA_2.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#Assign a calendar date to each time point\n",
    "startday = datetime.datetime(1951,1,1)\n",
    "date_list = np.array([datetime.timedelta(days=x) + startday for x in range(0, 20819)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#quick function that bins everything in a 2D array and puts floor and ceiling on array values.\n",
    "def findbins(lat):\n",
    "    index = np.ndarray.round(lat-19.5)\n",
    "    index[index < 0] = 0\n",
    "    index[index > 21] = 21\n",
    "    return index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "//anaconda/lib/python3.5/site-packages/ipykernel/__main__.py:4: RuntimeWarning: invalid value encountered in less\n",
      "//anaconda/lib/python3.5/site-packages/ipykernel/__main__.py:5: RuntimeWarning: invalid value encountered in greater\n"
     ]
    }
   ],
   "source": [
    "#bin latitudes\n",
    "lat_1_indices = findbins(lat_1_all)\n",
    "lat_2_indices = findbins(lat_2_all)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "//anaconda/lib/python3.5/site-packages/ipykernel/__main__.py:7: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future\n",
      "//anaconda/lib/python3.5/site-packages/ipykernel/__main__.py:11: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future\n"
     ]
    }
   ],
   "source": [
    "#create a time series of whether a rainband exists for each latitude bin\n",
    "occupancy_1 = np.zeros([20819,22])\n",
    "occupancy_2 = np.zeros([20819,22])\n",
    "\n",
    "for i, index in enumerate(lat_1_indices):\n",
    "    if ~np.isnan(index):\n",
    "        occupancy_1[i,index] = 1\n",
    "        \n",
    "for i, index in enumerate(lat_2_indices):\n",
    "    if ~np.isnan(index):\n",
    "        occupancy_2[i,index] = 1\n",
    "\n",
    "#also create combined occupancy matrix of both primary and secondary fronts\n",
    "occupancy = occupancy_1 + occupancy_2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def smooth(P,dayrange,latrange):\n",
    "    Psmooth=np.zeros(P.shape)\n",
    "    \n",
    "    ll=int((dayrange-1)/2)\n",
    "    yy=int((latrange-1)/2)\n",
    "    \n",
    "    days = P.shape[0]\n",
    "    lats = P.shape[1]\n",
    "    \n",
    "    for d in range(days):\n",
    "        for j in range(lats):\n",
    "            \n",
    "            Psample = P.take(range(d-ll,d+ll+1), mode='wrap', axis=0)\n",
    "            Ps = Psample[:,max(0,j-yy):min(days,j+yy+1)]\n",
    "            Psmooth[d,j] = np.mean(np.mean(Ps))\n",
    "    \n",
    "    return Psmooth"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#switching it to binary - hence the *5. Wanted to keep same def of smooth function as elsewhere\n",
    "#NOTE: in some cases, values can reach 2 (if both fronts are within small latitude range on a given day)\n",
    "#should be a minimal effect - we test both keeping the 2 and reducing to 1.\n",
    "occupancy_1_smth = smooth(occupancy_1,1,5)*5\n",
    "occupancy_2_smth = smooth(occupancy_2,1,5)*5\n",
    "occupancy_smth = smooth(occupancy,1,5)*5\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#as constructed, occupancy_smth can have values of 2 if fronts are close enough - this fixes that.\n",
    "occupancy_smth_binary = np.copy(occupancy_smth) #without using copy, still points to same array.\n",
    "occupancy_smth_binary[occupancy_smth_binary>1] = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.00048033  0.00086459  0.00384264  0.02877179  0.02656227  0.03227821\n",
      "  0.03443969  0.03818627  0.04303761  0.03669725  0.03247034  0.02987656\n",
      "  0.01902109  0.01998175  0.01565877  0.01508238  0.01196023  0.00979874\n",
      "  0.00667659  0.00614823  0.00403478  0.00134493]\n",
      "[ 0.00864595  0.04244921  0.06052164  0.09231952  0.12589462  0.16023824\n",
      "  0.17450406  0.18463903  0.18483116  0.18026802  0.16110284  0.13804698\n",
      "  0.1170085   0.09962054  0.08170421  0.07248187  0.05917671  0.04966617\n",
      "  0.03861857  0.02800327  0.02275566  0.01921322]\n"
     ]
    }
   ],
   "source": [
    "#quick comparison of how much the smoothing in latitude changes mean front occupancy\n",
    "print(np.mean(occupancy, axis=0))\n",
    "print(np.mean(occupancy_smth, axis=0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ac_1 = np.array([ autocorr(x)[1] for x in occupancy_1.T ])\n",
    "ac_2 = np.array([ autocorr(x)[1] for x in occupancy_2.T ])\n",
    "ac = np.array([ autocorr(x)[1] for x in occupancy.T ])\n",
    "\n",
    "ac_smth_1 = np.array([ autocorr(x)[1] for x in occupancy_1_smth.T ])\n",
    "ac_smth_2 = np.array([ autocorr(x)[1] for x in occupancy_2_smth.T ])\n",
    "ac_smth = np.array([ autocorr(x)[1] for x in occupancy_smth.T ])\n",
    "\n",
    "ac_smth_binary = np.array([ autocorr(x)[1] for x in occupancy_smth_binary.T ])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 1.26800427  2.06563188  2.19680495  2.41971398  2.63872406  2.91905467\n",
      "  3.13065542  3.14236761  3.09594696  3.02263046  2.7621429   2.7015617\n",
      "  2.78921635  3.09993799  3.48683284  3.55399918  3.32081316  3.06894202\n",
      "  2.65484529  2.4286421   2.09128261  1.93877854]\n"
     ]
    }
   ],
   "source": [
    "print(ac_smth_binary)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 365\n",
      "[ 1.26786473  2.06815904  2.19667512  2.41851382  2.63731245  2.91657888\n",
      "  3.12864151  3.14013304  3.0929943   3.02005989  2.75926658  2.69866058\n",
      "  2.78678746  3.09818152  3.48510832  3.55233922  3.31952351  3.067708\n",
      "  2.65380782  2.42789364  2.0907967   1.93847172]\n",
      "121 303\n",
      "[ 1.25910986  1.7751358   1.77983667  1.94169548  2.24122216  2.55674112\n",
      "  3.13345049  3.18741965  3.04275965  3.07778672  2.65459302  2.55384506\n",
      "  2.59913389  2.80532258  3.00901351  2.93672505  2.70481988  2.43793793\n",
      "  2.09884922  2.01092165  1.80647995  1.75284482]\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/Siwen/RDA/autocorr.py:21: RuntimeWarning: invalid value encountered in true_divide\n",
      "  yy_norm = (yy - mn) / stdv\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "304 120\n",
      "[ 0.97307575  1.30752568  1.83982645  2.31805864  2.54438857  2.91290018\n",
      "  2.94216576  2.9485238   2.98618978  2.76847747  2.57404321  2.36503662\n",
      "  2.10530157  1.75479571  1.70420055  1.52759612  1.31999537  1.19659115\n",
      "  1.13164305  1.09312178  1.10815138  0.98879426]\n",
      "61 120\n",
      "[ 0.93190695  1.20765554  1.61502509  1.97290452  1.87763004  1.93885051\n",
      "  1.73141673  1.71222539  1.8471951   1.95880991  1.90398223  1.85987348\n",
      "  1.84289566  1.55381197  1.66999157  1.5074714   1.15597072  1.07008494\n",
      "  1.08501374  1.0341087   0.95000157  0.97677   ]\n",
      "121 160\n",
      "[ 1.1175325   1.15846921  1.2845494   1.31815761  1.56142744  1.61892266\n",
      "  1.7742786   1.82963536  1.3958954   1.61416402  1.54258473  1.61677349\n",
      "  1.82639008  1.6085412   1.36110305  1.33357804  1.14662994  1.15135891\n",
      "  1.14879411  1.14182157  1.27146641  1.17023637]\n",
      "161 200\n",
      "[ 1.08979805  1.3084742   1.42795773  1.51423239  1.80624523  2.25555689\n",
      "  2.80670439  3.08561437  3.22276015  3.08813804  2.67759279  2.4756915\n",
      "  2.63995025  2.8925691   3.12244753  2.92989564  2.8554656   2.62409303\n",
      "  2.19509008  1.83795318  1.52364984  1.37011803]\n",
      "201 273\n",
      "[ 1.41455676  1.78289896  1.81158746  1.83063489  1.84718891  1.85113519\n",
      "  1.75079816  1.74320927  1.89871866  2.27421508  2.28560598  2.0207719\n",
      "  1.9709001   2.01801737  2.22604411  2.23898853  2.12298377  1.86656824\n",
      "  1.59036905  1.7579992   1.66926665  1.75005726]\n",
      "274 320\n",
      "[ 1.2696393   1.5608861   1.43832242  1.31893855  1.44495032  1.67643644\n",
      "  2.12005921  2.0448747   2.04140934  1.96120433  1.68654988  1.85655765\n",
      "  1.96641987  1.93473357  2.02219429  2.10715746  1.67711595  1.64839005\n",
      "  1.47377921  1.37167232  1.28704678  1.39090205]\n",
      "321 60\n",
      "[ 0.99212767  1.09661944  1.65620644  2.06909904  2.20224112  2.56558637\n",
      "  2.55681926  2.37002794  2.40524803  2.25336107  2.1493343   2.07527497\n",
      "  1.79373969  1.48817364  1.37638818  1.27046372  1.37264037  1.16100152\n",
      "  0.97031696  0.98743158         nan         nan]\n"
     ]
    }
   ],
   "source": [
    "## FIND AUTOCORRELATION FOR ARBITRARY TIME PERIODS ##\n",
    "#much like the RDA_freq_diff script, now taking the next step and letting autocorrelation be calculated over arbitrary\n",
    "#time periods.\n",
    "\n",
    "#it's reasonable for autocorrelation to change a lot across the year, so worth checking.\n",
    "day_beg = [1,121,304,61,121,161,201,274,321]\n",
    "day_end = [365,303,120,120,160,200,273,320,60]\n",
    "\n",
    "#stores autocorrelation time scales for each latitude box.\n",
    "tau_1 = [] #primary front only\n",
    "tau_2 = [] #secondary front only\n",
    "tau = [] #both primary and secondary\n",
    "\n",
    "for db, de in zip(day_beg, day_end):\n",
    "    \n",
    "    #date_list defined above as list of dates between January 1, 1951 and Dec 31, 2007\n",
    "    #mysterious .timetuple().tm_yday function below returns day of year given date.\n",
    "    if db <= de:\n",
    "        mymap = map(lambda date: (date.timetuple().tm_yday >= db) & (date.timetuple().tm_yday <= de), date_list)\n",
    "    else: #if de is earlier in the year than db, then we have to wrap around the new year.\n",
    "        mymap = map(lambda date: not((date.timetuple().tm_yday > de) & (date.timetuple().tm_yday < db)), date_list)\n",
    "        \n",
    "    extract = np.array(list(mymap))\n",
    "    \n",
    "    ## primary fronts\n",
    "    oc_1 = np.copy(occupancy_1_smth)\n",
    "    oc_1_mean = np.mean(occupancy_1_smth[extract],0) #this is the mean during time period of interest\n",
    "    oc_1[~extract] = oc_1_mean #change all values outside of time period of interest to mean\n",
    "    #therefore, won't contribute to autocorrelation\n",
    "    \n",
    "    ## secondary fronts\n",
    "    oc_2 = np.copy(occupancy_2_smth)\n",
    "    oc_2_mean = np.mean(occupancy_2_smth[extract],0) #this is the mean during time period of interest\n",
    "    oc_2[~extract] = oc_2_mean #change all values outside of time period of interest to mean\n",
    "    #therefore, won't contribute to autocorrelation\n",
    "    \n",
    "    ## all fronts\n",
    "    oc = np.copy(occupancy_smth_binary)\n",
    "    oc_mean = np.mean(occupancy_smth_binary[extract],0) #this is the mean during time period of interest\n",
    "    oc[~extract] = oc_mean #change all values outside of time period of interest to mean\n",
    "    #therefore, won't contribute to autocorrelation\n",
    "    \n",
    "    tau_1.append(np.array([ autocorr(x)[1] for x in oc_1.T ]))\n",
    "    tau_2.append(np.array([ autocorr(x)[1] for x in oc_2.T ]))\n",
    "    tau.append(np.array([ autocorr(x)[1] for x in oc.T ]))\n",
    "    \n",
    "    print(db,de)\n",
    "    #print(tau_1[-1])\n",
    "    #print(tau_2[-1])\n",
    "    print(tau[-1])\n",
    "    \n",
    "    #small debugging module to verify that the actual standard deviation of our hacked time series is the same as\n",
    "    #the theoretical expectation for a Bernoulli time-series (p*(1-p))**.5\n",
    "    #print((oc_mean*(1-oc_mean))**.5)\n",
    "    #print(np.std(oc,0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "## clear from the above that autocorrelation of rainbands changes dramatically between seasons.\n",
    "# the yearly behavior is dominated by what goes on during Meiyu season, when the autocorrelation is super-high,\n",
    "# but then the decorrelation time scale is significantly longer outside of it.\n",
    "\n",
    "#Therefore, let's experiment with both - use a single yearly autocorrelation time scale value, and another, highly\n",
    "#seasonalized value, and see what it does to the resulting frequency difference graphs.\n",
    "\n",
    "## UPDATE - the seasonalized value seems to let changes come through more clearly. Otherwise, losing information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 2.75926658  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658\n",
      "  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658\n",
      "  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658\n",
      "  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658\n",
      "  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658\n",
      "  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658\n",
      "  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658\n",
      "  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658\n",
      "  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658\n",
      "  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658\n",
      "  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658\n",
      "  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658\n",
      "  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658\n",
      "  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658\n",
      "  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658\n",
      "  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658\n",
      "  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658\n",
      "  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658\n",
      "  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658\n",
      "  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658\n",
      "  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658\n",
      "  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658\n",
      "  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658\n",
      "  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658\n",
      "  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658\n",
      "  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658\n",
      "  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658\n",
      "  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658\n",
      "  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658\n",
      "  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658\n",
      "  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658\n",
      "  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658\n",
      "  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658\n",
      "  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658\n",
      "  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658\n",
      "  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658\n",
      "  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658\n",
      "  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658\n",
      "  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658\n",
      "  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658\n",
      "  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658\n",
      "  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658\n",
      "  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658\n",
      "  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658\n",
      "  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658\n",
      "  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658\n",
      "  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658\n",
      "  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658\n",
      "  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658\n",
      "  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658\n",
      "  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658\n",
      "  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658\n",
      "  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658\n",
      "  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658\n",
      "  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658\n",
      "  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658\n",
      "  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658\n",
      "  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658\n",
      "  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658\n",
      "  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658\n",
      "  2.75926658  2.75926658  2.75926658  2.75926658  2.75926658]\n",
      "[ 2.57404321  2.57404321  2.57404321  2.57404321  2.57404321  2.57404321\n",
      "  2.57404321  2.57404321  2.57404321  2.57404321  2.57404321  2.57404321\n",
      "  2.57404321  2.57404321  2.57404321  2.57404321  2.57404321  2.57404321\n",
      "  2.57404321  2.57404321  2.57404321  2.57404321  2.57404321  2.57404321\n",
      "  2.57404321  2.57404321  2.57404321  2.57404321  2.57404321  2.57404321\n",
      "  2.57404321  2.57404321  2.57404321  2.57404321  2.57404321  2.57404321\n",
      "  2.57404321  2.57404321  2.57404321  2.57404321  2.57404321  2.57404321\n",
      "  2.57404321  2.57404321  2.57404321  2.57404321  2.57404321  2.57404321\n",
      "  2.57404321  2.57404321  2.57404321  2.57404321  2.57404321  2.57404321\n",
      "  2.57404321  2.57404321  2.57404321  2.57404321  2.57404321  2.57404321\n",
      "  2.57404321  2.57404321  2.57404321  2.57404321  2.57404321  2.57404321\n",
      "  2.57404321  2.57404321  2.57404321  2.57404321  2.57404321  2.57404321\n",
      "  2.57404321  2.57404321  2.57404321  2.57404321  2.57404321  2.57404321\n",
      "  2.57404321  2.57404321  2.57404321  2.57404321  2.57404321  2.57404321\n",
      "  2.57404321  2.57404321  2.57404321  2.57404321  2.57404321  2.57404321\n",
      "  2.57404321  2.57404321  2.57404321  2.57404321  2.57404321  2.57404321\n",
      "  2.57404321  2.57404321  2.57404321  2.57404321  2.57404321  2.57404321\n",
      "  2.57404321  2.57404321  2.57404321  2.57404321  2.57404321  2.57404321\n",
      "  2.57404321  2.57404321  2.57404321  2.57404321  2.57404321  2.57404321\n",
      "  2.57404321  2.57404321  2.57404321  2.57404321  2.57404321  2.57404321\n",
      "  2.65459302  2.65459302  2.65459302  2.65459302  2.65459302  2.65459302\n",
      "  2.65459302  2.65459302  2.65459302  2.65459302  2.65459302  2.65459302\n",
      "  2.65459302  2.65459302  2.65459302  2.65459302  2.65459302  2.65459302\n",
      "  2.65459302  2.65459302  2.65459302  2.65459302  2.65459302  2.65459302\n",
      "  2.65459302  2.65459302  2.65459302  2.65459302  2.65459302  2.65459302\n",
      "  2.65459302  2.65459302  2.65459302  2.65459302  2.65459302  2.65459302\n",
      "  2.65459302  2.65459302  2.65459302  2.65459302  2.65459302  2.65459302\n",
      "  2.65459302  2.65459302  2.65459302  2.65459302  2.65459302  2.65459302\n",
      "  2.65459302  2.65459302  2.65459302  2.65459302  2.65459302  2.65459302\n",
      "  2.65459302  2.65459302  2.65459302  2.65459302  2.65459302  2.65459302\n",
      "  2.65459302  2.65459302  2.65459302  2.65459302  2.65459302  2.65459302\n",
      "  2.65459302  2.65459302  2.65459302  2.65459302  2.65459302  2.65459302\n",
      "  2.65459302  2.65459302  2.65459302  2.65459302  2.65459302  2.65459302\n",
      "  2.65459302  2.65459302  2.65459302  2.65459302  2.65459302  2.65459302\n",
      "  2.65459302  2.65459302  2.65459302  2.65459302  2.65459302  2.65459302\n",
      "  2.65459302  2.65459302  2.65459302  2.65459302  2.65459302  2.65459302\n",
      "  2.65459302  2.65459302  2.65459302  2.65459302  2.65459302  2.65459302\n",
      "  2.65459302  2.65459302  2.65459302  2.65459302  2.65459302  2.65459302\n",
      "  2.65459302  2.65459302  2.65459302  2.65459302  2.65459302  2.65459302\n",
      "  2.65459302  2.65459302  2.65459302  2.65459302  2.65459302  2.65459302\n",
      "  2.65459302  2.65459302  2.65459302  2.65459302  2.65459302  2.65459302\n",
      "  2.65459302  2.65459302  2.65459302  2.65459302  2.65459302  2.65459302\n",
      "  2.65459302  2.65459302  2.65459302  2.65459302  2.65459302  2.65459302\n",
      "  2.65459302  2.65459302  2.65459302  2.65459302  2.65459302  2.65459302\n",
      "  2.65459302  2.65459302  2.65459302  2.65459302  2.65459302  2.65459302\n",
      "  2.65459302  2.65459302  2.65459302  2.65459302  2.65459302  2.65459302\n",
      "  2.65459302  2.65459302  2.65459302  2.65459302  2.65459302  2.65459302\n",
      "  2.65459302  2.65459302  2.65459302  2.65459302  2.65459302  2.65459302\n",
      "  2.65459302  2.65459302  2.65459302  2.65459302  2.65459302  2.65459302\n",
      "  2.65459302  2.65459302  2.65459302  2.65459302  2.65459302  2.65459302\n",
      "  2.65459302  2.65459302  2.65459302  2.57404321  2.57404321  2.57404321\n",
      "  2.57404321  2.57404321  2.57404321  2.57404321  2.57404321  2.57404321\n",
      "  2.57404321  2.57404321  2.57404321  2.57404321  2.57404321  2.57404321\n",
      "  2.57404321  2.57404321  2.57404321  2.57404321  2.57404321  2.57404321\n",
      "  2.57404321  2.57404321  2.57404321  2.57404321  2.57404321  2.57404321\n",
      "  2.57404321  2.57404321  2.57404321  2.57404321  2.57404321  2.57404321\n",
      "  2.57404321  2.57404321  2.57404321  2.57404321  2.57404321  2.57404321\n",
      "  2.57404321  2.57404321  2.57404321  2.57404321  2.57404321  2.57404321\n",
      "  2.57404321  2.57404321  2.57404321  2.57404321  2.57404321  2.57404321\n",
      "  2.57404321  2.57404321  2.57404321  2.57404321  2.57404321  2.57404321\n",
      "  2.57404321  2.57404321  2.57404321  2.57404321  2.57404321]\n",
      "[ 2.1493343   2.1493343   2.1493343   2.1493343   2.1493343   2.1493343\n",
      "  2.1493343   2.1493343   2.1493343   2.1493343   2.1493343   2.1493343\n",
      "  2.1493343   2.1493343   2.1493343   2.1493343   2.1493343   2.1493343\n",
      "  2.1493343   2.1493343   2.1493343   2.1493343   2.1493343   2.1493343\n",
      "  2.1493343   2.1493343   2.1493343   2.1493343   2.1493343   2.1493343\n",
      "  2.1493343   2.1493343   2.1493343   2.1493343   2.1493343   2.1493343\n",
      "  2.1493343   2.1493343   2.1493343   2.1493343   2.1493343   2.1493343\n",
      "  2.1493343   2.1493343   2.1493343   2.1493343   2.1493343   2.1493343\n",
      "  2.1493343   2.1493343   2.1493343   2.1493343   2.1493343   2.1493343\n",
      "  2.1493343   2.1493343   2.1493343   2.1493343   2.1493343   2.1493343\n",
      "  1.90398223  1.90398223  1.90398223  1.90398223  1.90398223  1.90398223\n",
      "  1.90398223  1.90398223  1.90398223  1.90398223  1.90398223  1.90398223\n",
      "  1.90398223  1.90398223  1.90398223  1.90398223  1.90398223  1.90398223\n",
      "  1.90398223  1.90398223  1.90398223  1.90398223  1.90398223  1.90398223\n",
      "  1.90398223  1.90398223  1.90398223  1.90398223  1.90398223  1.90398223\n",
      "  1.90398223  1.90398223  1.90398223  1.90398223  1.90398223  1.90398223\n",
      "  1.90398223  1.90398223  1.90398223  1.90398223  1.90398223  1.90398223\n",
      "  1.90398223  1.90398223  1.90398223  1.90398223  1.90398223  1.90398223\n",
      "  1.90398223  1.90398223  1.90398223  1.90398223  1.90398223  1.90398223\n",
      "  1.90398223  1.90398223  1.90398223  1.90398223  1.90398223  1.90398223\n",
      "  1.54258473  1.54258473  1.54258473  1.54258473  1.54258473  1.54258473\n",
      "  1.54258473  1.54258473  1.54258473  1.54258473  1.54258473  1.54258473\n",
      "  1.54258473  1.54258473  1.54258473  1.54258473  1.54258473  1.54258473\n",
      "  1.54258473  1.54258473  1.54258473  1.54258473  1.54258473  1.54258473\n",
      "  1.54258473  1.54258473  1.54258473  1.54258473  1.54258473  1.54258473\n",
      "  1.54258473  1.54258473  1.54258473  1.54258473  1.54258473  1.54258473\n",
      "  1.54258473  1.54258473  1.54258473  1.54258473  2.67759279  2.67759279\n",
      "  2.67759279  2.67759279  2.67759279  2.67759279  2.67759279  2.67759279\n",
      "  2.67759279  2.67759279  2.67759279  2.67759279  2.67759279  2.67759279\n",
      "  2.67759279  2.67759279  2.67759279  2.67759279  2.67759279  2.67759279\n",
      "  2.67759279  2.67759279  2.67759279  2.67759279  2.67759279  2.67759279\n",
      "  2.67759279  2.67759279  2.67759279  2.67759279  2.67759279  2.67759279\n",
      "  2.67759279  2.67759279  2.67759279  2.67759279  2.67759279  2.67759279\n",
      "  2.67759279  2.67759279  2.28560598  2.28560598  2.28560598  2.28560598\n",
      "  2.28560598  2.28560598  2.28560598  2.28560598  2.28560598  2.28560598\n",
      "  2.28560598  2.28560598  2.28560598  2.28560598  2.28560598  2.28560598\n",
      "  2.28560598  2.28560598  2.28560598  2.28560598  2.28560598  2.28560598\n",
      "  2.28560598  2.28560598  2.28560598  2.28560598  2.28560598  2.28560598\n",
      "  2.28560598  2.28560598  2.28560598  2.28560598  2.28560598  2.28560598\n",
      "  2.28560598  2.28560598  2.28560598  2.28560598  2.28560598  2.28560598\n",
      "  2.28560598  2.28560598  2.28560598  2.28560598  2.28560598  2.28560598\n",
      "  2.28560598  2.28560598  2.28560598  2.28560598  2.28560598  2.28560598\n",
      "  2.28560598  2.28560598  2.28560598  2.28560598  2.28560598  2.28560598\n",
      "  2.28560598  2.28560598  2.28560598  2.28560598  2.28560598  2.28560598\n",
      "  2.28560598  2.28560598  2.28560598  2.28560598  2.28560598  2.28560598\n",
      "  2.28560598  2.28560598  2.28560598  1.68654988  1.68654988  1.68654988\n",
      "  1.68654988  1.68654988  1.68654988  1.68654988  1.68654988  1.68654988\n",
      "  1.68654988  1.68654988  1.68654988  1.68654988  1.68654988  1.68654988\n",
      "  1.68654988  1.68654988  1.68654988  1.68654988  1.68654988  1.68654988\n",
      "  1.68654988  1.68654988  1.68654988  1.68654988  1.68654988  1.68654988\n",
      "  1.68654988  1.68654988  1.68654988  1.68654988  1.68654988  1.68654988\n",
      "  1.68654988  1.68654988  1.68654988  1.68654988  1.68654988  1.68654988\n",
      "  1.68654988  1.68654988  1.68654988  1.68654988  1.68654988  1.68654988\n",
      "  1.68654988  1.68654988  2.1493343   2.1493343   2.1493343   2.1493343\n",
      "  2.1493343   2.1493343   2.1493343   2.1493343   2.1493343   2.1493343\n",
      "  2.1493343   2.1493343   2.1493343   2.1493343   2.1493343   2.1493343\n",
      "  2.1493343   2.1493343   2.1493343   2.1493343   2.1493343   2.1493343\n",
      "  2.1493343   2.1493343   2.1493343   2.1493343   2.1493343   2.1493343\n",
      "  2.1493343   2.1493343   2.1493343   2.1493343   2.1493343   2.1493343\n",
      "  2.1493343   2.1493343   2.1493343   2.1493343   2.1493343   2.1493343\n",
      "  2.1493343   2.1493343   2.1493343   2.1493343   2.1493343 ]\n"
     ]
    }
   ],
   "source": [
    "## based on autocorrelation estimates above, create 3 sets of seasonally changing decorrelation time scales:\n",
    "# 1) constant\n",
    "# 2) half-years\n",
    "# 3) seasonal, as defined in paper.\n",
    "\n",
    "#each is defined for each latitude bin.\n",
    "\n",
    "tau_constant = np.zeros([365,22])\n",
    "tau_halfyear = np.zeros([365,22])\n",
    "tau_seasonal = np.zeros([365,22])\n",
    "\n",
    "#relies on the list of days defined earlier in this script:\n",
    "#day_beg = [1,121,304,61,121,161,201,274,321]\n",
    "#day_end = [365,303,120,120,160,200,273,320,60]\n",
    "\n",
    "\n",
    "## CONSTANT YEARLY VALUE FOR TAU\n",
    "for dd in range(365):\n",
    "    tau_constant[dd,:] = tau[0]\n",
    "    \n",
    "print(tau_constant[:,10])\n",
    "            \n",
    "\n",
    "## HALF-YEAR TAU VALUES\n",
    "for mytau, db, de in zip(tau[1:3],day_beg[1:3],day_end[1:3]):\n",
    "    \n",
    "    if db <= de:\n",
    "        for dd in range(db-1,de):\n",
    "            tau_halfyear[dd,:] = mytau\n",
    "        \n",
    "    else:\n",
    "        for dd in range(db-1,de+365):\n",
    "            tau_halfyear[dd%365,:] = mytau\n",
    "            \n",
    "print(tau_halfyear[:,10])            \n",
    "\n",
    "\n",
    "## SEASONAL TAU VALUES\n",
    "for mytau, db, de in zip(tau[3:],day_beg[3:],day_end[3:]):\n",
    "    \n",
    "    if db <= de:\n",
    "        for dd in range(db-1,de):\n",
    "            tau_seasonal[dd,:] = mytau\n",
    "        \n",
    "    else:\n",
    "        for dd in range(db-1,de+365):\n",
    "            tau_seasonal[dd%365,:] = mytau\n",
    "            \n",
    "print(tau_seasonal[:,10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## SAVE AUTOCORRELATIONS BY LATITUDE BIN AND SEASON AS NETCDF\n",
    "#after considering differences between primary and primary+secondary, just doing the latter\n",
    "outfilename = \"RDA_freq_tau.nc\"\n",
    "ferretpath = \"/Users/Siwen/Desktop/Ferret/bin/\"\n",
    "\n",
    "#netCDF output file formatting formatting\n",
    "fullname = ferretpath + outfilename\n",
    "\n",
    "##the following is only to make sure that we don't exceed the bounds of the number of rows of actual data.\n",
    "if os.path.isfile(fullname):\n",
    "    os.remove(fullname)\n",
    "    \n",
    "outfile = nc.Dataset(fullname, 'w')\n",
    "\n",
    "#create dimensions in NetCDF file\n",
    "days_dim = outfile.createDimension('time', None)\n",
    "lat_dim = outfile.createDimension('lat', 22)\n",
    "dim_tuple = (\"time\",\"lat\")\n",
    "\n",
    "tau_const = outfile.createVariable(\"tau_constant\", 'f8', dim_tuple)\n",
    "tau_hy = outfile.createVariable(\"tau_halfyear\", 'f8', dim_tuple)\n",
    "tau_seas = outfile.createVariable(\"tau_seasonal\", 'f8', dim_tuple)\n",
    "tau_isband = outfile.createVariable(\"tau_band\", 'f8', dim_tuple)\n",
    "\n",
    "tau_const[:] = tau_constant\n",
    "tau_hy[:] = tau_halfyear\n",
    "tau_seas[:] = tau_seasonal\n",
    "tau_isband[:] = tau_band\n",
    "\n",
    "\n",
    "outfile.close()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}