acf.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import numpy\n",
    "import xarray as xr\n",
    "import matplotlib.pyplot as plot\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def normalize_axis(axis,ndim):\n",
    "    \n",
    "    # axis is either a scalar or iterable. \n",
    "    #\n",
    "    try:\n",
    "        iterator = iter(axis)\n",
    "    except TypeError:\n",
    "        # not iterable\n",
    "        axis = [axis]\n",
    "  \n",
    "    \n",
    "    axis = [x if x >= 0 else ndim+x for x in axis]\n",
    "    \n",
    "    return axis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def read_var(file, ftype, varname, ts, dims=[1,1,1]):\n",
    "    print('file = ', file)\n",
    "    if ftype == \"netcdf\":\n",
    "        ds = xr.open_dataset(file, decode_times=False)\n",
    "        data = ds.get(varname).isel(time=ts).values\n",
    "    else:\n",
    "        nelements = np.prod(dims)\n",
    "        data = np.fromfile(file, dtype=np.float32, count=nelements)\n",
    "        data = np.squeeze(data.reshape(dims))\n",
    "\n",
    "    return (data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def autocorr1(x,lags):\n",
    "    '''numpy.corrcoef, partial'''\n",
    "\n",
    "    corr=[1. if l==0 else numpy.corrcoef(x[l:],x[:-l])[0][1] for l in lags]\n",
    "    return numpy.array(corr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def autocorr2(x,lags):\n",
    "    '''manualy compute, non partial'''\n",
    "\n",
    "    mean=numpy.mean(x)\n",
    "    var=numpy.var(x)\n",
    "    xp=x-mean\n",
    "    corr=[1. if l==0 else numpy.sum(xp[l:]*xp[:-l])/len(x)/var for l in lags]\n",
    "\n",
    "    return numpy.array(corr)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "#\n",
    "# FFT-based autocorrelation with padding\n",
    "# \n",
    "def autocorr3(x,lags,axis=-1):\n",
    "    '''fft, pad 0s, non partial'''\n",
    "\n",
    "    n=x.shape[axis]\n",
    "    # pad 0s to 2n-1\n",
    "    ext_size=2*n-1\n",
    "    # nearest power of 2\n",
    "    fsize=2**numpy.ceil(numpy.log2(ext_size)).astype('int')\n",
    "\n",
    "    xp=x-numpy.mean(x,axis=axis,keepdims=True)\n",
    "    var=numpy.var(x,axis=axis,keepdims=True)\n",
    "\n",
    "    # do fft and ifft\n",
    "    cf=numpy.fft.fft(xp,fsize,axis=axis)\n",
    "    sf=cf.conjugate()*cf\n",
    "    corr=numpy.fft.ifft(sf,axis=axis).real\n",
    "    corr = np.divide(corr,var,where=(var!=0.0))\n",
    "    corr = corr / n\n",
    "   \n",
    "    naxis = tuple(set(np.arange(x.ndim).tolist()) - set(normalize_axis(axis,x.ndim)))\n",
    "    \n",
    "    return np.mean(corr,axis=naxis)[:len(lags)]\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def autocorr4(x,lags):\n",
    "    '''fft, don't pad 0s, non partial'''\n",
    "    mean=x.mean()\n",
    "    var=numpy.var(x)\n",
    "    xp=x-mean\n",
    "\n",
    "    cf=numpy.fft.fft(xp)\n",
    "    sf=cf.conjugate()*cf\n",
    "    corr=numpy.fft.ifft(sf).real/var/len(x)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def autocorr5(x,lags):\n",
    "    '''numpy.correlate, non partial'''\n",
    "    mean=x.mean()\n",
    "    var=numpy.var(x)\n",
    "    xp=x-mean\n",
    "    corr=numpy.correlate(xp,xp,'full')[len(x)-1:]/var/len(x)\n",
    "\n",
    "    return corr[:len(lags)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(64, 64)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fe8138e11c0>]"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD4CAYAAADlwTGnAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAodUlEQVR4nO3dd3hUddrG8e+TRu+EIh2pIVJjQEpioYuA2LAiCojS2WJZX3d9131FdzcURRFFFAtYQEGkBdxN6BB67wgBhAACSg/83j8y7sYYIULImWTuz3V5TeaUmfuAyT3nzOTBnHOIiEjgCfI6gIiIeEMFICISoFQAIiIBSgUgIhKgVAAiIgEqxOsAv0Xp0qVd1apVvY4hIpKrrFix4rBzLjzj8lxVAFWrViUpKcnrGCIiuYqZfZvZcl0CEhEJUCoAEZEApQIQEQlQKgARkQClAhARCVCXLQAze9fMDpnZ+l9Zb2Y2ysy2m9laM2ucbl17M9viW/dMuuUlzSzezLb5bktkz+GIiEhWZeUM4D2g/SXWdwBq+v7rA7wJYGbBwGjf+gjgfjOL8O3zDDDPOVcTmOe7LyIiOeiyBeCcSwSOXmKTLsAEl2YJUNzMygPRwHbn3E7n3Dlgkm/bn/Z53/f1+0DXK8yfJUt2HmHcgl1cuKjR1yIiP8mO9wAqAHvT3U/2Lfu15QBlnXMHAHy3ZX7twc2sj5klmVlSSkrKFQX8eu0B/jp9I3ePWcS2gz9c0WOIiOQ12VEAlskyd4nlv4lzbqxzLso5FxUe/ovfZM6S/+1SjxH3NWT34ZPcPmoBo+Zt41zqxSt6LBGRvCI7CiAZqJTufkVg/yWWAxz0XSbCd3soG3L8KjOja6MKxA+NpV1kOeLit9L59QWsTT52LZ9WRMSvZUcBTAMe8X0aqBlw3HdZZzlQ08yqmVkY0N237U/79PB93QOYmg05Lqt04Xy8dn8j3n4kiu9PnaPr6IW8PGMTp89dyImnFxHxK5cdBmdmE4GbgdJmlgz8GQgFcM6NAWYAHYHtwCmgp29dqpn1B2YDwcC7zrkNvocdBnxqZo8De4B7svGYLqtNRFmiq5Vk2MxNvJW4k9kbvmPYXfVpVr1UTsYQEfGU5aZ/FD4qKspl9zTQRdsP88yUdew5eooHm1bmmQ51KJI/NFufQ0TES2a2wjkXlXF5wP8mcPMapZk9OIZeLasxcdke2g5P5JvNB72OJSJyzQV8AQAUCAvm+U4RTH6yOUXyh/DYe0kMnrSKoyfPeR1NROSaUQGk06hyCaYPaMWg22ry9boDtI5LYNqa/eSmy2QiIlmlAsggLCSIIW1q8dWAllQqUYCBE1fRe8IKvjt+xutoIiLZSgXwK+qUK8qUp1rw/O11WbA9hTZxCUxctkdnAyKSZ6gALiE4yOjVqjqzB8cQWaEYz05ZxwNvL+XbIye9jiYictVUAFlQpVQhPu7dlJe73cD6fcdpNyKRd+bv1HA5EcnVVABZZGbcH12Z+KGxtKxRmpe+3kS3Nxex5TsNlxOR3EkF8BuVK5aftx+JYtT9jUg+eopOr81nxNytGi4nIrmOCuAKmBmdG1xH/NBYbr+hPCPmbuOO1xaweu8xr6OJiGSZCuAqlCwUxojujXj30ShOnDlPtzcW8tL0jRouJyK5ggogG9xapyxzhsRwf3Rl3lmwi3YjElm047DXsURELkkFkE2K5A/lb3fewKQ+zQgyeODtpTw7ZS0nzpz3OpqISKZUANmsWfVSzBwUwxMx1flk+V7axCUwd6OGy4mI/1EBXAMFwoJ5tmNdvuzXghIFw+g1IYkBE1dx5MezXkcTEfkPFcA1VL9icab1b8nQNrWYtT5tuNyXq/ZpnISI+AUVwDUWFhLEwNtq8vXAVlQpVYjBn6zm8feT2H/stNfRRCTAqQBySK2yRZj8ZHP+p1MEi3ccoe3wRD5c8i0XNU5CRDyiAshBwUHG4y2rMXtwDA0qFeP5L9dz/9tL2HVYw+VEJOepADxQuVRBPny8Ka/cdQMbD5yg/YhE3krYQeoFjZMQkZyjAvCImXHfjZWZOzSWmFrhvDxzM93eXMSmAye8jiYiAUIF4LGyRfMz9uEmvP5AI/YfO80dry0gbs4WzqZqnISIXFsqAD9gZnSqfx3xQ2Lp3OA6Rn2zndtHLWDFt997HU1E8jAVgB8pUSiMuPsaMr7njZw6m8rdYxbx4lcbOHUu1etoIpIHqQD80C21yzB7SAwPNa3C+IW7aTs8kQXbNFxORLKXCsBPFckfyl+7RvLpEzcRGhzEQ+OW8sfP13D8lIbLiUj2UAH4uehqJZk5qBVP3nw9k1fuo/XwBGat/87rWCKSB2SpAMysvZltMbPtZvZMJutLmNkXZrbWzJaZWWS6dYPMbL2ZbTCzwemWNzSzJWa22sySzCw6W44oD8ofGszT7eswtV8Lwgvno++HK+j30UpSftBwORG5cpctADMLBkYDHYAI4H4zi8iw2XPAaudcfeARYKRv30igNxANNAA6mVlN3z6vAi865xoCL/juyyVEVijG1P4t+EO72sRvPEjruAQmr0jWcDkRuSJZOQOIBrY753Y6584Bk4AuGbaJAOYBOOc2A1XNrCxQF1jinDvlnEsFEoA7ffs4oKjv62LA/qs6kgARGhxEv1tqMGNQK2qUKczvPlvDo+OXs0/D5UTkN8pKAVQA9qa7n+xblt4aoBuA71JOFaAisB6IMbNSZlYQ6AhU8u0zGPi7me0F/gE8m9mTm1kf3yWipJSUlCwdVCCoUaYwnz1xE3+5I4Llu4/SNi6BCYt3a7iciGRZVgrAMlmW8afMMKCEma0GBgCrgFTn3CbgFSAemEVaUfz0ofYngSHOuUrAEGBcZk/unBvrnItyzkWFh4dnIW7gCAoyHm2RNlyucZUSvDB1A/eNXcyOlB+9jiYiuUBWCiCZ/75qh7RX9j+7XOOcO+Gc6+m7nv8IEA7s8q0b55xr7JyLAY4C23y79QCm+L7+jLRLTXIFKpUsyITHovnHPQ3YevBHOoyczxv/3s55DZcTkUvISgEsB2qaWTUzCwO6A9PSb2BmxX3rAHoBic65E751ZXy3lUm7TDTRt91+INb39a38txjkCpgZdzepSPzQGG6rU4ZXZ22h6+iFrN933OtoIuKnQi63gXMu1cz6A7OBYOBd59wGM+vrWz+GtDd7J5jZBWAj8Hi6h5hsZqWA80A/59xPA256AyPNLAQ4A/TJroMKZGWK5OfNh5owc90B/mfqBrqMXkjf2OoMuLUm+UODvY4nIn7EctNHCKOiolxSUpLXMXKNY6fO8bevN/HZimSqhxfi1bvqE1W1pNexRCSHmdkK51xUxuX6TeA8rHjBMP5+TwMmPBbN2fMXueetxfxl2gZOntVwORFRAQSEmFrhzBkSQ4+bqvL+4rThcglb9ZFakUCnAggQhfKF8JfO9fjsiZvIHxpEj3eX8btP13Ds1Dmvo4mIR1QAASaqakm+HtiK/rfU4MvV+2gdl8jMdQe8jiUiHlABBKD8ocH8vl1tpvVvQdmi+Xjyo5X0/WAFh06c8TqaiOQgFUAAq3ddMab2a8HT7evwzZZDtI5L4LOkvRouJxIgVAABLiQ4iCdvvp6Zg1pRu1wR/vD5Wh55dxl7j57yOpqIXGMqAAHg+vDCfNLnJv7apR4rv/2ediMSGb9wFxc0XE4kz1IByH8EBRkP31SV2UNiuLFqSV78aiP3vrWY7Yd+8DqaiFwDKgD5hYolCvJezxuJu7cBO1J+pOPIBbz+zTYNlxPJY1QAkikzo1vjisQPiaVNvbL8Y85WOr+u4XIieYkKQC4pvEg+Rj/QmLcebsLhH8/SZfRChs3czJnzF7yOJiJXSQUgWdKuXjnmDonl7sYVGZOwg44j57Ns11GvY4nIVVABSJYVKxjKK3fX58PHm3LuwkXufWsx//Plen44c97raCJyBVQA8pu1rFmaOUNieKxFNT5c+i3thifyry2HvI4lIr+RCkCuSMGwEF64I4LP+zanYL4Qeo5fztBPVvP9SQ2XE8ktVAByVZpUKcHXA1sy8NYaTFuzn9ZxCUxfu1/jJERyARWAXLV8IcEMbVubrwa05LriBej/8Sqe+GAFBzVcTsSvqQAk29QtX5QvnmrOsx3qkLA1hdZxCXyyfI/OBkT8lApAslVIcBBPxF7PrMEx1C1flKcnr+OhcUvZc0TD5UT8jQpArolqpQsxqXczXuoayZq9x2k3IpFxCzRcTsSfqADkmgkKMh5qVoU5Q2JoVr0kf52+kbvHLGLbQQ2XE/EHKgC55q4rXoB3H72Rkd0bsvvwSW4ftYBR87ZxLlXD5US8pAKQHGFmdGlYgblDY2kXWY64+K10fn0Ba5OPeR1NJGCpACRHlSqcj9fub8Tbj0Tx/alzdB29kJdnbOL0OQ2XE8lpKgDxRJuIssQPjeW+GyvxVuJOOoxMZMnOI17HEgkoKgDxTNH8obzcrT4f92rKRQfdxy7hT1+s03A5kRySpQIws/ZmtsXMtpvZM5msL2FmX5jZWjNbZmaR6dYNMrP1ZrbBzAZn2G+A73E3mNmrV300kis1r1Ga2YNj6N2qGhOX7aHt8ES+2XzQ61gied5lC8DMgoHRQAcgArjfzCIybPYcsNo5Vx94BBjp2zcS6A1EAw2ATmZW07fuFqALUN85Vw/4R7YckeRKBcKC+dPtEUx5qgVF84fy2HtJDJq0iiM/nvU6mkielZUzgGhgu3Nup3PuHDCJtB/c6UUA8wCcc5uBqmZWFqgLLHHOnXLOpQIJwJ2+fZ4Ehjnnzvr20zxhoWGl4nw1oCWDW9dkxroDtBmeyLQ1Gi4nci1kpQAqAHvT3U/2LUtvDdANwMyigSpARWA9EGNmpcysINARqOTbpxbQysyWmlmCmd2Y2ZObWR8zSzKzpJSUlKwel+RiYSFBDG5di+kDWlGpZEEGTlxF7wlJfHdcw+VEslNWCsAyWZbx5dgwoISZrQYGAKuAVOfcJuAVIB6YRVpRpPr2CQFKAM2APwCfmtkvnss5N9Y5F+WciwoPD89CXMkrapcrwpQnm/P87XVZsP0wbeISmLhMw+VEsktWCiCZ/75qh7RX9vvTb+CcO+Gc6+mca0jaewDhwC7funHOucbOuRjgKLAt3eNOcWmWAReB0ldzMJL3BAcZvVpVZ/bgGCIrFOPZKet44O2lfHvkpNfRRHK9rBTAcqCmmVUzszCgOzAt/QZmVty3DqAXkOicO+FbV8Z3W5m0y0QTfdt9CdzqW1cLCAMOX9XRSJ5VpVQhPu7dlGHdbmD9vrThcm8n7tRwOZGrEHK5DZxzqWbWH5gNBAPvOuc2mFlf3/oxpL3ZO8HMLgAbgcfTPcRkMysFnAf6Oee+9y1/F3jXzNYD54AeTuf2cglmRvfoytxcuwzPf7mOv83YxPS1+3n17gbULlfE63giuY7lpp+5UVFRLikpyesY4gecc0xfe4C/TNvAiTPneermGvS7pQZhIfrdRpGMzGyFcy4q43J9t0iuZGbc0eA64ofGcvsN5Rk5bxudXpvPqj3fX35nEQFUAJLLlSwUxojujXj30Sh+OJPKXW8u4qXpGzVcTiQLVACSJ9xapyxzhsTwQNPKvLNgF+1GJLJouz5TIHIpKgDJM4rkD+WlrjcwqU8zggweeGcpz0xey/HTGi4nkhkVgOQ5zaqXYtbgGJ6Irc6nSXtpOzyB+I0aLieSkQpA8qT8ocE826EuX/ZrQYmCYfSekET/j1dyWMPlRP5DBSB5Wv2KxZnWvyW/a1OLORsO0iYugS9X7dM4CRFUABIAwkKCGHBbTb4e2JKqpQsx+JPVPP5+EvuPnfY6moinVAASMGqWLcLnfZvzQqcIFu84QtvhiXy45FsuapyEBCgVgASU4CDjsZbVmDMkhoaVivP8l+vp/vYSdh3WcDkJPCoACUiVShbkg8ejefWu+mw6cIL2IxIZk7CD1AsXvY4mkmNUABKwzIx7b6zE3KGxxNYKZ9jMzdz5xiI27j/hdTSRHKECkIBXtmh+3nq4CaMfaMyB46fp/PoC/jlnC2dTNU5C8jYVgAhpZwO31y9P/JBYOje8jte+2c7toxaw4lsNl5O8SwUgkk6JQmHE3duQ93reyOlzF7h7zCJe/GoDJ8+mXn5nkVxGBSCSiZtrl2H2kBgeblaF8Qt3025EIvO3pXgdSyRbqQBEfkXhfCH8b5dIPn3iJsKCg3h43DL++Pkajp/ScDnJG1QAIpcRXa0kMwa14smbr2fyyn20Hp7ArPXfeR1L5KqpAESyIH9oME+3r8PUfi0IL5yPvh+uoN9HK0n5QcPlJPdSAYj8BpEVijG1fwv+0K428ZsO0jougckrkjVcTnIlFYDIbxQaHES/W2owY2ArapQpzO8+W0OP8ctJ/v6U19FEfhMVgMgVqlGmMJ89cRMvdq5H0u6jtBueyITFuzVcTnINFYDIVQgKMno0r8rswTE0rlKCF6Zu4L6xi9mR8qPX0UQuSwUgkg0qlSzIhMei+cc9Ddh68Ec6jJzPG//eznkNlxM/pgIQySZmxt1NKhI/NIbWdcvw6qwtdB29kPX7jnsdTSRTKgCRbFamSH7eeLAJYx5qzMETZ+kyeiGvztrMmfMaLif+RQUgco20jyzPvKGxdGtUgTf+vYOOo+aTtPuo17FE/iNLBWBm7c1si5ltN7NnMllfwsy+MLO1ZrbMzCLTrRtkZuvNbIOZDc5k39+bmTOz0ld1JCJ+qFjBUP5+TwMmPBbN2fMXueetxfx56np+1HA58QOXLQAzCwZGAx2ACOB+M4vIsNlzwGrnXH3gEWCkb99IoDcQDTQAOplZzXSPXQloA+y5+kMR8V8xtcKZMySGHjdVZcKSb2k3PJGErRouJ97KyhlANLDdObfTOXcOmAR0ybBNBDAPwDm3GahqZmWBusAS59wp51wqkADcmW6/4cAfAX1wWvK8QvlC+Evnenze9ybyhwbR491l/O7TNRw7dc7raBKgslIAFYC96e4n+5altwboBmBm0UAVoCKwHogxs1JmVhDoCFTybdcZ2OecW3OpJzezPmaWZGZJKSl6xSS5X5MqJfl6YCv631KDqav30TougRnrDngdSwJQVgrAMlmW8RX7MKCEma0GBgCrgFTn3CbgFSAemEVaUaT6yuBPwAuXe3Ln3FjnXJRzLio8PDwLcUX8X/7QYH7frjZT+7egXLH8PPXRSvp+sIJDJ854HU0CSFYKIBnfq3afisD+9Bs4504453o65xqS9h5AOLDLt26cc66xcy4GOApsA64HqgFrzGy37zFXmlm5qzsckdyl3nXF+PKpFjzdvg7fbDlE67gEPk3aq+FykiOyUgDLgZpmVs3MwoDuwLT0G5hZcd86gF5AonPuhG9dGd9tZdIuE010zq1zzpVxzlV1zlUlrWQaO+c0ZF0CTkhwEE/efD2zBrWiTrmi/PHztTzy7jL2HtVwObm2LlsAvjdv+wOzgU3Ap865DWbW18z6+jarC2wws82kfVpoULqHmGxmG4GvgH7OOf0r2yKZqB5emEl9mvHXrpGs/PZ72o1IZPzCXVzQcDm5Riw3nWpGRUW5pKQkr2OIXHP7jp3mT1+s499bUmhcuTiv3l2fGmWKeB1LcikzW+Gci8q4XL8JLOKHKhQvwPhHb2T4fQ3YefgkHUcu4PVvtmm4nGQrFYCInzIz7mxUkblDY2lTryz/mLOVO15bwLpkDZeT7KECEPFzpQvnY/QDjXnr4SYcPXmOrm8sZNhMDZeTq6cCEMkl2tUrR/zQWO5uXJExCTvoMHI+S3ce8TqW5GIqAJFcpFiBUF65uz4f9WpK6sWL3Dd2Cc9/uY4fzpz3OprkQioAkVyoRY3SzB4cw2MtqvHR0j20G57Iv7Yc8jqW5DIqAJFcqmBYCC/cEcHkJ5tTKF8IPccvZ8gnqzl6UsPlJGtUACK5XOPKJZg+sCUDb6vJV2v20yYugelr92uchFyWCkAkD8gXEszQNrX4akBLKpQoQP+PV9HngxUc1HA5uQQVgEgeUrd8UaY82ZznOtYhcWsKreMS+GT5Hp0NSKZUACJ5TEhwEH1irmf24Bgiyhfl6cnrePCdpew5ouFy8nMqAJE8qmrpQkzs3Yy/3RnJ2uTjtBuRyLgFGi4n/6UCEMnDgoKMB5tWIX5oDDddX4q/Tt/IXW8uYuvBH7yOJn5ABSASAMoXK8C4HlGM7N6Qb4+c5PZR8xk1bxvnUjVcLpCpAEQChJnRpWEF5g6NpUNkeeLit9L59QWs2XvM62jiERWASIApVTgfo+5vxDuPRHHs1HnufGMhL8/YxOlzGi4XaFQAIgGqdURZ5gyN4b4bK/NW4k46jExkiYbLBRQVgEgAK5o/lJe73cDHvZvigO5jl/DcF+s4oeFyAUEFICI0v740swbF0KtlNSYt20PbuES+2XzQ61hyjakARASAAmHBPN8pgilPtaBYgVAeey+JQZNWceTHs15Hk2tEBSAiP9OwUnG+GtCSwa1rMmPdAdoMT2Tq6n0aJ5EHqQBE5BfCQoIY3LoW0we0olLJggyatJpe7ydx4Phpr6NJNlIBiMivql2uCFOebM7zt9dl4Y7DtI1L5OOle7iocRJ5ggpARC4pOMjo1ao6swfHEFmhGM99sY4H3lnC7sMnvY4mV0kFICJZUqVUIT7u3ZRh3W5gw74TtB+ZyNuJO0m9oHESuZUKQESyzMzoHl2Z+KGxtKwRzt9mbOKuNxex+bsTXkeTK6ACEJHfrFyx/Lz9SBNeu78Ryd+fptOoBcTFb+VsqsZJ5CZZKgAza29mW8xsu5k9k8n6Emb2hZmtNbNlZhaZbt0gM1tvZhvMbHC65X83s82+fb4ws+LZcUAikjPMjDsaXEf80FjuaHAdo+Zt447XFrBqz/deR5MsumwBmFkwMBroAEQA95tZRIbNngNWO+fqA48AI337RgK9gWigAdDJzGr69okHIn37bAWevfrDEZGcVrJQGMPva8j4R2/khzOpdHtzEX+dvpFT51K9jiaXkZUzgGhgu3Nup3PuHDAJ6JJhmwhgHoBzbjNQ1czKAnWBJc65U865VCABuNO33RzfMoAlQMWrPhoR8cwtdcowZ0gMDzatzLgFu2g3IpGF2w97HUsuISsFUAHYm+5+sm9ZemuAbgBmFg1UIe0H+nogxsxKmVlBoCNQKZPneAyYmdmTm1kfM0sys6SUlJQsxBURrxTJH8pLXW/gkz7NCAkK4sF3lvL052s5flrD5fxRVgrAMlmW8bdAhgElzGw1MABYBaQ65zYBr5B2uWcWaUXxs/NCM/uTb9lHmT25c26scy7KORcVHh6ehbgi4rWm1Usxc1Arnoitzucrk2kTl8CcDd95HUsyyEoBJPPzV+0Vgf3pN3DOnXDO9XTONSTtPYBwYJdv3TjnXGPnXAxwFNj2035m1gPoBDzoNGhEJE/JHxrMsx3q8uVTLShVOB99PlhB/49XcljD5fxGVgpgOVDTzKqZWRjQHZiWfgMzK+5bB9ALSHTOnfCtK+O7rUzaZaKJvvvtgaeBzs65U9lxMCLif26oWIxp/Vvw+7a1mLPhIK3jEvhiVbKGy/mByxaA743a/sBsYBPwqXNug5n1NbO+vs3qAhvMbDNpnxYalO4hJpvZRuAroJ9z7qfPiL0OFAHizWy1mY3JnkMSEX8TGhxE/1trMmNQS6qXLsSQT9bQ873l7Dum4XJestzUwlFRUS4pKcnrGCJyFS5cdExYvJtXZ20hyOCZjnV5MLoyQUGZvd0o2cHMVjjnojIu128Ci0iOCg4yeraoxpwhMTSqXIL/+XI93ccuYWfKj15HCzgqABHxRKWSBfng8Whevbs+m787QYeR8xmTsEPD5XKQCkBEPGNm3BtViblDY7m5djjDZm6m6xsL2bhfw+VyggpARDxXpmh+xjzUhDcebMx3x8/Q+fUF/HPOFg2Xu8ZUACLiF8yMjjeUJ35ILF0aVuC1b7Zz+6gFrPhWw+WuFRWAiPiVEoXC+Oe9DXj/sWhOn7vA3WMW8eJXGzh5VsPlspsKQET8UmytcGYPieGRZlUYv3A37UYkMn+b5oFlJxWAiPitwvlCeLFLJJ/1vYmwkCAeHreMP3y2huOnNFwuO6gARMTv3Vi1JDMGtuKpm69nyqp9tB6ewKz1Gi53tVQAIpIr5A8N5o/t6zC1XwvCC+ej74creOqjFRz64YzX0XItFYCI5CqRFYoxtX8L/tCuNnM3HaJNXCKTV2i43JVQAYhIrhMaHES/W2owY2ArapYpzO8+W0OP8ctJ/l6DhX8LFYCI5Fo1yhTm0ydu4sXO9UjafZS2wxN5f9FuLl7U2UBWqABEJFcLCjJ6NK/KnCExRFUtyZ+nbeDetxazQ8PlLksFICJ5QsUSBXm/5438854GbDv0Ix1Gzmf0v7ZzXsPlfpUKQETyDDPjriYVmTs0ltZ1y/D32Vvo8vpC1u877nU0v6QCEJE8J7xIPt54sAljHmpMyo9n6TJ6Ia/M2syZ8xoul54KQETyrPaR5Zk7JJZujSrw5r930HHkfJbvPup1LL+hAhCRPK1YwVD+fk8DPng8mnMXLnLPmMW8MHU9P2q4nApARAJDq5rhzB4cQ88WVflgybe0G57Iv7cc8jqWp1QAIhIwCuUL4c931OPzvs0pEBbMo+OXM/TT1Xx/8pzX0TyhAhCRgNOkSgm+HtiSgbfWYNrq/bQZnsCMdQcCbpyECkBEAlK+kGCGtq3NtP4tKV+sAE99tJK+H67g0InAGS6nAhCRgBZxXVG+eKo5z3aow7+3pNA6LoFPk/YGxNmACkBEAl5IcBBPxF7PzEGtqFO+KH/8fC0Pj1vG3qN5e7icCkBExKd6eGEm9W7GS10jWb33GG2HJ/Lugl1cyKPD5VQAIiLpBAUZDzWrwpwhMTStXpL/nb6Re8YsYtvBH7yOlu2yVABm1t7MtpjZdjN7JpP1JczsCzNba2bLzCwy3bpBZrbezDaY2eB0y0uaWbyZbfPdlsiWIxIRyQbXFS/A+EdvZMR9Ddl1+CS3j1rAa/O25anhcpctADMLBkYDHYAI4H4zi8iw2XPAaudcfeARYKRv30igNxANNAA6mVlN3z7PAPOcczWBeb77IiJ+w8zo2qgC8UNjaVuvLP+M38odry1gXXLeGC6XlTOAaGC7c26nc+4cMAnokmGbCNJ+iOOc2wxUNbOyQF1giXPulHMuFUgA7vTt0wV43/f1+0DXqzkQEZFrpXThfLz+QGPGPtyEoyfP0WX0Al6euSnXD5fLSgFUAPamu5/sW5beGqAbgJlFA1WAisB6IMbMSplZQaAjUMm3T1nn3AEA322ZzJ7czPqYWZKZJaWkpGTtqEREroG29coRPzSWe6Mq8VbCTjqMnM/SnUe8jnXFslIAlsmyjG+JDwNKmNlqYACwCkh1zm0CXgHigVmkFcVvmsDknBvrnItyzkWFh4f/ll1FRLJdsQKhDLurPh/1asqFi477xi7h+S/X8cOZ815H+82yUgDJ/PdVO6S9st+ffgPn3AnnXE/nXEPS3gMIB3b51o1zzjV2zsUAR4Ftvt0Omll5AN9tYE9lEpFcpUWN0swa3IpeLavx8dI9tBueyL82564fY1kpgOVATTOrZmZhQHdgWvoNzKy4bx1ALyDROXfCt66M77YyaZeJJvq2mwb08H3dA5h6NQciIpLTCoaF8HynCCY/2ZxC+ULo+d5yhnyymqO5ZLjcZQvA9+Ztf2A2sAn41Dm3wcz6mllf32Z1gQ1mtpm0TwsNSvcQk81sI/AV0M85971v+TCgjZltA9r47ouI5DqNKpdg+sCWDLqtJl+t2U+buAS+WrPf78dJmL8HTC8qKsolJSV5HUNE5Fdt/u4ET3++ljXJx2ldtyx/uzOSskXze5rJzFY456IyLtdvAouIZKM65Yoy+cnm/KljXeZvSxsuN2nZHr88G1ABiIhks5DgIHrHVGf24BgiyhflmSnrePCdpew54l/D5VQAIiLXSNXShZjYuxn/d+cNrEs+TtsRCbwzf6ffDJdTAYiIXENBQcYDTSszZ2gMLa4vzUtfb6Lbm4vY8p33w+VUACIiOaB8sQK80yOKkd0bsvfoKTq9Np8Rc7dyLtW74XIqABGRHGJmdGlYgfghMXS8oTwj5m7jjtcWsGbvMU/yqABERHJYqcL5GNm9EeN6RHH89HnufGMhf/t6I6fP5exwORWAiIhHbqtbljlDY+geXZm35++i/chEFu04nGPPrwIQEfFQ0fyh/N+dNzCxdzMAHnh7Kc9OWceJHBgupwIQEfEDN11filmDYugTU51Plu+hTVwCczcevKbPqQIQEfETBcKCea5jXb54qgUlCobRa0ISAyeu4siPZ6/J86kARET8TINKxZnWvyVD29Ri5voDtI5LYPGO7P+HZ1QAIiJ+KCwkiIG31eTrga2IrFCMqqULZvtzhGT7I4qISLapVbYIHzze9Jo8ts4AREQClApARCRAqQBERAKUCkBEJECpAEREApQKQEQkQKkAREQClApARCRAmT/+S/W/xsxSgG+vcPfSQM7NWb0y/p7R3/OB/2f093ygjNnB3/JVcc6FZ1yYqwrgaphZknMuyuscl+LvGf09H/h/Rn/PB8qYHfw93090CUhEJECpAEREAlQgFcBYrwNkgb9n9Pd84P8Z/T0fKGN28Pd8QAC9ByAiIj8XSGcAIiKSjgpARCRABUQBmFl7M9tiZtvN7Bmv86RnZpXM7F9mtsnMNpjZIK8z/RozCzazVWY23essGZlZcTP73Mw2+/4sb/I6U0ZmNsT3d7zezCaaWX4/yPSumR0ys/XplpU0s3gz2+a7LeFn+f7u+3tea2ZfmFlxr/L58vwiY7p1vzczZ2alvch2OXm+AMwsGBgNdAAigPvNLMLbVD+TCvzOOVcXaAb087N86Q0CNnkd4leMBGY55+oADfCznGZWARgIRDnnIoFgoLu3qQB4D2ifYdkzwDznXE1gnu++V97jl/nigUjnXH1gK/BsTofK4D1+mREzqwS0AfbkdKCsyvMFAEQD251zO51z54BJQBePM/2Hc+6Ac26l7+sfSPvBVcHbVL9kZhWB24F3vM6SkZkVBWKAcQDOuXPOuWOehspcCFDAzEKAgsB+j/PgnEsEjmZY3AV43/f1+0DXnMyUXmb5nHNznHOpvrtLgIo5HuzneTL7MwQYDvwR8NtP2gRCAVQA9qa7n4wf/oAFMLOqQCNgqcdRMjOCtP+ZL3qcIzPVgRRgvO8S1TtmVsjrUOk55/YB/yDt1eAB4Lhzbo63qX5VWefcAUh7gQKU8TjPpTwGzPQ6REZm1hnY55xb43WWSwmEArBMlvldI5tZYWAyMNg5d8LrPOmZWSfgkHNuhddZfkUI0Bh40znXCDiJt5ctfsF3Hb0LUA24DihkZg95myp3M7M/kXYJ9SOvs6RnZgWBPwEveJ3lcgKhAJKBSunuV8QPTr3TM7NQ0n74f+Scm+J1nky0ADqb2W7SLqHdamYfehvpZ5KBZOfcT2dOn5NWCP6kNbDLOZfinDsPTAGae5zp1xw0s/IAvttDHuf5BTPrAXQCHnT+98tM15NW9Gt83zMVgZVmVs7TVJkIhAJYDtQ0s2pmFkbaG2/TPM70H2ZmpF273uSci/M6T2acc8865yo656qS9uf3jXPOb169Oue+A/aaWW3fotuAjR5GysweoJmZFfT9nd+Gn71Rnc40oIfv6x7AVA+z/IKZtQeeBjo75055nScj59w651wZ51xV3/dMMtDY9/+pX8nzBeB7s6g/MJu0b7hPnXMbvE31My2Ah0l7Vb3a919Hr0PlQgOAj8xsLdAQ+D9v4/yc7+zkc2AlsI607z3PxwWY2URgMVDbzJLN7HFgGNDGzLaR9imWYX6W73WgCBDv+34Z41W+S2TMFTQKQkQkQOX5MwAREcmcCkBEJECpAEREApQKQEQkQKkAREQClApARCRAqQBERALU/wPvujQ85b/TKQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plot\n",
    "\n",
    "t = np.arange(0, 64*64, 1,dtype=np.float32) / 64.0 * np.pi * 25.0\n",
    "s = np.sin(t).reshape(64,64)\n",
    "\n",
    "print (s.shape)\n",
    "n = 16\n",
    "lags = range(n)\n",
    "\n",
    "plot.plot(autocorr3(t,lags,axis=-1))\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CLDTOT   0.0   1.0\n",
      "         dim = 0   [0.99999996 0.87487303 0.76533522 0.68190394]\n",
      "         dim = 1   [0.99869905 0.89372246 0.79388989 0.71467519]\n",
      "FLDS   94.57982   459.7665\n",
      "         dim = 0   [0.99999999 0.98678879 0.97468588 0.96366206]\n",
      "         dim = 1   [1.         0.95852748 0.92165327 0.88889622]\n",
      "FLDSC   94.16072   431.92297\n",
      "         dim = 0   [1.00000001 0.99632354 0.99179718 0.98696969]\n",
      "         dim = 1   [1.         0.99268171 0.9794218  0.96418709]\n",
      "FLNS   -80.99102   435.33313\n",
      "         dim = 0   [1.00000002 0.95611181 0.91626015 0.8821644 ]\n",
      "         dim = 1   [1.         0.94786605 0.89823612 0.85454549]\n",
      "FLUT   80.489204   374.09515\n",
      "         dim = 0   [1.00000002 0.97721298 0.94316017 0.90783344]\n",
      "         dim = 1   [1.         0.96967847 0.92413242 0.87539718]\n",
      "FLUTC   173.52068   376.0187\n",
      "         dim = 0   [1.         0.99505005 0.9885495  0.9818061 ]\n",
      "         dim = 1   [1.         0.9890874  0.97100861 0.95109461]\n",
      "FSDS   0.0   1104.7214\n",
      "         dim = 0   [0.99999993 0.98134859 0.96216939 0.94413644]\n",
      "         dim = 1   [1.         0.98877617 0.97521367 0.96174705]\n",
      "FSDSC   0.0   1105.4778\n",
      "         dim = 0   [1.00000018 0.99327633 0.98648944 0.97966932]\n",
      "         dim = 1   [1.         0.99832446 0.99609686 0.99370695]\n",
      "FSNS   0.0   1074.9822\n",
      "         dim = 0   [1.00000001 0.9724596  0.94392309 0.91828166]\n",
      "         dim = 1   [1.         0.98820293 0.9739142  0.95976989]\n",
      "FSNSC   0.0   1075.5436\n",
      "         dim = 0   [0.99999988 0.99051507 0.97924684 0.96798234]\n",
      "         dim = 1   [1.         0.99805455 0.99538615 0.99256676]\n",
      "FSNTOA   0.0   1327.1697\n",
      "         dim = 0   [0.99999996 0.9866243  0.97302563 0.96044673]\n",
      "         dim = 1   [1.         0.99392303 0.98674246 0.97960542]\n",
      "FSNTOAC   0.0   1327.1697\n",
      "         dim = 0   [1.00000018 0.99408107 0.98744635 0.9807201 ]\n",
      "         dim = 1   [1.         0.99816487 0.99568002 0.99304404]\n",
      "LHFLX   -230.00731   924.3833\n",
      "         dim = 0   [1.         0.97973442 0.95018745 0.92065756]\n",
      "         dim = 1   [1.         0.97395217 0.93938567 0.90549762]\n",
      "OMEGA500   -8.563181   3.4370985\n",
      "         dim = 0   [1.00000001 0.8678716  0.67351946 0.53348116]\n",
      "         dim = 1   [1.         0.89743527 0.73335588 0.59413851]\n",
      "OMEGA850   -7.199399   6.226275\n",
      "         dim = 0   [1.00000002 0.81509618 0.54035252 0.34442633]\n",
      "         dim = 1   [1.         0.87596378 0.65476877 0.46963822]\n",
      "PRECC   0.0   1.0735267e-06\n",
      "         dim = 0   [0.99999969 0.88528274 0.74131623 0.61874711]\n",
      "         dim = 1   [0.8972247  0.77028356 0.62993582 0.51055571]\n",
      "PRECL   -1.6854496e-20   1.8852295e-05\n",
      "         dim = 0   [1.         0.7392424  0.51367172 0.3682822 ]\n",
      "         dim = 1   [0.99869904 0.8000277  0.59652904 0.46127777]\n",
      "PRECSC   0.0   1.1438748e-07\n",
      "         dim = 0   [0.44480601 0.19084303 0.08797585 0.0512994 ]\n",
      "         dim = 1   [0.37263698 0.26433051 0.19939193 0.15347832]\n",
      "PRECSL   -1.3576364e-19   3.5453468e-06\n",
      "         dim = 0   [0.99999996 0.793574   0.61343179 0.48802743]\n",
      "         dim = 1   [0.98699054 0.54681812 0.44280655 0.3759205 ]\n",
      "PRECT   -1.6854496e-20   1.8852295e-05\n",
      "         dim = 0   [0.99999998 0.82871157 0.63064182 0.48510033]\n",
      "         dim = 1   [0.99869905 0.86140082 0.67924896 0.53685144]\n",
      "PREC_pcw   -1.6854496e-20   1.8852295e-05\n",
      "         dim = 0   [1.         0.7392424  0.51367172 0.3682822 ]\n",
      "         dim = 1   [0.99869904 0.8000277  0.59652904 0.46127777]\n",
      "PREC_zmc   0.0   1.0735267e-06\n",
      "         dim = 0   [0.99999969 0.88528274 0.74131623 0.61874711]\n",
      "         dim = 1   [0.8972247  0.77088139 0.63352111 0.51791933]\n",
      "PS   49771.2   103770.83\n",
      "         dim = 0   [1.00000001 0.98627163 0.9662682  0.9437164 ]\n",
      "         dim = 1   [0.99869905 0.98582282 0.95626873 0.92171148]\n",
      "PSL   95243.63   106444.81\n",
      "         dim = 0   [0.99999997 0.99785502 0.99479433 0.99115704]\n",
      "         dim = 1   [0.99869905 0.99441121 0.98802246 0.98137674]\n",
      "QREFHT   3.3239467e-05   0.023403471\n",
      "         dim = 0   [0.99999998 0.99666233 0.99218381 0.98747995]\n",
      "         dim = 1   [0.99869905 0.99005478 0.97714318 0.96336531]\n",
      "RHREFHT   14.376854   2262.7522\n",
      "         dim = 0   [1.         0.96803656 0.93121324 0.89890837]\n",
      "         dim = 1   [1.         0.96673324 0.9192534  0.87463613]\n",
      "SHFLX   -152.55318   2514.837\n",
      "         dim = 0   [1.00000001 0.96033859 0.91817421 0.88401355]\n",
      "         dim = 1   [1.         0.96789838 0.92942892 0.8953728 ]\n",
      "SOLIN   0.0   1408.253\n",
      "         dim = 0   [0.99999979 0.99421916 0.98840343 0.98255361]\n",
      "         dim = 1   [0.99869905 0.99324417 0.98638276 0.97998823]\n",
      "T500   226.54243   271.30182\n",
      "         dim = 0   [0.99999999 0.99843636 0.99662813 0.99459704]\n",
      "         dim = 1   [0.99869904 0.99305351 0.98575248 0.97804533]\n",
      "T850   227.19751   302.54214\n",
      "         dim = 0   [1.00000001 0.99618395 0.99167138 0.98697688]\n",
      "         dim = 1   [1.         0.98758612 0.96776959 0.94829757]\n",
      "TGCLDIWP   0.0   0.7850221\n",
      "         dim = 0   [0.9999999  0.90455323 0.78394074 0.67756305]\n",
      "         dim = 1   [0.99609714 0.91950408 0.81731186 0.7203355 ]\n",
      "TGCLDLWP   0.0   2.2505817\n",
      "         dim = 0   [0.99999991 0.88084232 0.74242739 0.62581999]\n",
      "         dim = 1   [0.98811915 0.86923931 0.74987952 0.6486541 ]\n",
      "TMQ   0.20576538   75.97875\n",
      "         dim = 0   [1.00000001 0.996643   0.99052245 0.98306493]\n",
      "         dim = 1   [1.         0.99254494 0.97721561 0.95847421]\n",
      "TREFHT   215.13908   303.04626\n",
      "         dim = 0   [1.         0.99483072 0.9887749  0.98264136]\n",
      "         dim = 1   [1.         0.99091293 0.97659491 0.96161925]\n",
      "TREFHTMN   215.13908   303.04626\n",
      "         dim = 0   [1.         0.99483072 0.9887749  0.98264136]\n",
      "         dim = 1   [1.         0.99091293 0.97659491 0.96161925]\n",
      "TREFHTMX   215.13908   303.04626\n",
      "         dim = 0   [1.         0.99483072 0.9887749  0.98264136]\n",
      "         dim = 1   [1.         0.99091293 0.97659491 0.96161925]\n",
      "TREFMNAV   215.13908   303.04626\n",
      "         dim = 0   [1.         0.99483072 0.9887749  0.98264136]\n",
      "         dim = 1   [1.         0.99091293 0.97659491 0.96161925]\n",
      "TREFMXAV   215.13908   303.04626\n",
      "         dim = 0   [1.         0.99483072 0.9887749  0.98264136]\n",
      "         dim = 1   [1.         0.99091293 0.97659491 0.96161925]\n"
     ]
    }
   ],
   "source": [
    "file = '/Users/clyne/Data/CAM/AMIP_0.25degree/cam5_1_amip_run2.cam2.h1.1979-01-01-00000.nc'\n",
    "\n",
    "\n",
    "\n",
    "ds = xr.open_dataset(file, decode_times=False)\n",
    "for varname in ds.data_vars:\n",
    "\n",
    "    if ds[varname].ndim == 3:\n",
    "        data = ds.get(varname).isel(time=0).values\n",
    "        \n",
    "        lags=range(4)\n",
    "\n",
    "        # print 4-lag autocorrelation. Note lag=0 should be ~1.0\n",
    "        print (varname, ' ', np.min(data), ' ', np.max(data))\n",
    "        for d in (0,1):\n",
    "            print('        ', 'dim =', d, ' ', autocorr3(data,lags,axis=d))\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.8.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}