matt-long/etopo_grid_file.ipynb

## etopo_grid_file.ipynb
{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Demonstrate grid file generation in SCRIP format"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import xarray as xr\n",
    "from datetime import datetime\n",
    "\n",
    "from itertools import product\n",
    "\n",
    "import esmlab\n",
    "\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define a function to generate regular lat x lon grids."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def latlon_to_scrip(nx, ny, lon0=-180., grid_imask=None, file_out=None):\n",
    "    \"\"\"Generate a SCRIP grid file for a regular lat x lon grid.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    \n",
    "    nx : int\n",
    "       Number of points in x (longitude).\n",
    "    ny : int\n",
    "       Number of points in y (latitude).\n",
    "    lon0 : float, optional [default=-180]\n",
    "       Longitude on lefthand grid boundary.\n",
    "    grid_imask : array-like, optional [default=None]       \n",
    "       If the value is set to 0 for a grid point, then that point is\n",
    "       considered masked out and won't be used in the weights \n",
    "       generated by the application. \n",
    "    file_out : string, optional [default=None]\n",
    "       File to which to write the grid.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    \n",
    "    ds : xarray.Dataset\n",
    "       The grid file dataset.       \n",
    "    \"\"\"\n",
    "    \n",
    "    # compute coordinates of regular grid\n",
    "    dx = 360. / nx\n",
    "    dy = 180. / ny\n",
    "    lat = np.arange(-90. + dy / 2., 90., dy)\n",
    "    lon = np.arange(lon0 + dx / 2., lon0 + 360., dx)\n",
    "\n",
    "    # make 2D\n",
    "    y_center = np.broadcast_to(lat[:, None], (ny, nx))\n",
    "    x_center = np.broadcast_to(lon[None, :], (ny, nx))\n",
    "\n",
    "    # compute corner points: must be counterclockwise\n",
    "    y_corner = np.stack((y_center - dy / 2.,  # SW\n",
    "                         y_center - dy / 2.,  # SE\n",
    "                         y_center + dy / 2.,  # NE\n",
    "                         y_center + dy / 2.), # NW\n",
    "                        axis=2)\n",
    "\n",
    "    x_corner = np.stack((x_center - dx / 2.,  # SW\n",
    "                         x_center + dx / 2.,  # SE\n",
    "                         x_center + dx / 2.,  # NE\n",
    "                         x_center - dx / 2.), # NW\n",
    "                        axis=2)\n",
    "\n",
    "    # compute area\n",
    "    y0 = np.sin(y_corner[:, :, 0] * np.pi / 180.) # south\n",
    "    y1 = np.sin(y_corner[:, :, 3] * np.pi / 180.) # north\n",
    "    x0 = x_corner[:, :, 0] * np.pi / 180.         # west\n",
    "    x1 = x_corner[:, :, 1] * np.pi / 180.         # east\n",
    "    grid_area = (y1 - y0) * (x1 - x0)\n",
    "    \n",
    "    # sum of area should be equal to area of sphere\n",
    "    np.testing.assert_allclose(grid_area.sum(), 4.*np.pi)\n",
    "    \n",
    "    # construct mask\n",
    "    if grid_imask is None:\n",
    "        grid_imask = np.ones((ny, nx), dtype=np.int32)\n",
    "    \n",
    "    # generate output dataset\n",
    "    dso = xr.Dataset()    \n",
    "    dso['grid_dims'] = xr.DataArray(np.array([nx, ny], dtype=np.int32), \n",
    "                                    dims=('grid_rank',)) \n",
    "    dso.grid_dims.encoding = {'dtype': np.int32}\n",
    "\n",
    "    dso['grid_center_lat'] = xr.DataArray(y_center.reshape((-1,)), \n",
    "                                          dims=('grid_size'),\n",
    "                                          attrs={'units': 'degrees'})\n",
    "\n",
    "    dso['grid_center_lon'] = xr.DataArray(x_center.reshape((-1,)), \n",
    "                                          dims=('grid_size'),\n",
    "                                          attrs={'units': 'degrees'})\n",
    "    \n",
    "    dso['grid_corner_lat'] = xr.DataArray(y_corner.reshape((-1, 4)), \n",
    "                                          dims=('grid_size', 'grid_corners'), \n",
    "                                          attrs={'units': 'degrees'})\n",
    "    dso['grid_corner_lon'] = xr.DataArray(x_corner.reshape((-1, 4)), \n",
    "                                      dims=('grid_size', 'grid_corners'), \n",
    "                                      attrs={'units': 'degrees'})    \n",
    "\n",
    "    dso['grid_imask'] = xr.DataArray(grid_imask.reshape((-1,)), \n",
    "                                     dims=('grid_size'),\n",
    "                                     attrs={'units': 'unitless'})\n",
    "    dso.grid_imask.encoding = {'dtype': np.int32}\n",
    "    \n",
    "    dso['grid_area'] = xr.DataArray(grid_area.reshape((-1,)), \n",
    "                                     dims=('grid_size'),\n",
    "                                     attrs={'units': 'radians^2',\n",
    "                                            'long_name': 'area weights'})\n",
    "    \n",
    "    # force no '_FillValue' if not specified\n",
    "    for v in dso.variables:\n",
    "        if '_FillValue' not in dso[v].encoding:\n",
    "            dso[v].encoding['_FillValue'] = None\n",
    "\n",
    "    dso.attrs = {'title': f'{dy} x {dx} (lat x lon) grid',\n",
    "                 'created_by': 'latlon_to_scrip',\n",
    "                 'date_created': f'{datetime.now()}',\n",
    "                 'conventions': 'SCRIP',\n",
    "                }\n",
    "            \n",
    "    # write output file\n",
    "    if file_out is not None:\n",
    "        print(f'writing {file_out}')\n",
    "        dso.to_netcdf(file_out)\n",
    "        \n",
    "    return dso\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Test grid file generation with very simple, low resolution grid. Plot grid vertice (.) and grid corners (x)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "xarray.Dataset {\n",
      "dimensions:\n",
      "\tgrid_corners = 4 ;\n",
      "\tgrid_rank = 2 ;\n",
      "\tgrid_size = 54 ;\n",
      "\n",
      "variables:\n",
      "\tint32 grid_dims(grid_rank) ;\n",
      "\tfloat64 grid_center_lat(grid_size) ;\n",
      "\t\tgrid_center_lat:units = degrees ;\n",
      "\tfloat64 grid_center_lon(grid_size) ;\n",
      "\t\tgrid_center_lon:units = degrees ;\n",
      "\tfloat64 grid_corner_lat(grid_size, grid_corners) ;\n",
      "\t\tgrid_corner_lat:units = degrees ;\n",
      "\tfloat64 grid_corner_lon(grid_size, grid_corners) ;\n",
      "\t\tgrid_corner_lon:units = degrees ;\n",
      "\tint32 grid_imask(grid_size) ;\n",
      "\t\tgrid_imask:units = unitless ;\n",
      "\tfloat64 grid_area(grid_size) ;\n",
      "\t\tgrid_area:units = radians^2 ;\n",
      "\t\tgrid_area:long_name = area weights ;\n",
      "\n",
      "// global attributes:\n",
      "\t:title = 20.0 x 60.0 (lat x lon) grid ;\n",
      "\t:created_by = latlon_to_scrip ;\n",
      "\t:date_created = 2019-04-07 11:48:07.006968 ;\n",
      "\t:conventions = SCRIP ;\n",
      "}"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3X1cVPeZ///XZ87MgCAgKooigmiMBrxBkBATW5NGVxsTY7Zuk26qxp+5YXeTml/9tk2b1k3v0huzSc23tVHWu5omadIYExNdtYktiUEEQcR4UyGDyI2CIiAIw8z5fP84jGBMV2cGcjjyeT4eeczMkcN5G+U6x+vcXEJKiaIoinL9s5kdQFEURfliqIKvKIrSR6iCryiK0keogq8oitJHqIKvKIrSR6iCryiK0keogq8oitJHqIKvKIrSR6iCryiK0kfYzQ7Q1eDBg2ViYqLZMRRFUSyloKCgTkoZc7Wv61UFPzExkfz8fLNjKIqiWIoQovxavk61dBRFUfoIVfAVRVH6CFXwFUVR+ghV8BVFUfoIVfAVRVH6CMsX/PLylzhX//Fly87Vf0x5+UsmJQrMtm3PcbB422XLDhZvY9u250xKFLjiDTm4Pjxy2TLXh0co3pBjUqLA5bzw35zYufeyZSd27iXnhf82KVFgtj6/jQM7Lr8C7sCOfLY+v+0frNF7FX9nJa7Xd1+2zPX6boq/s9KkRIE7m51Nc+6+y5Y15+7jbHZ2j2zP8gU/InIiJSVPXCr65+o/pqTkCSIiJ5qczD/xI29k2zsfXSr6B4u3se2dj4gfeaPJyfwXOWYwnnerLxV914dH8LxbTeSYwSYn89+wm8azY8MLl4r+iZ172bHhBYbdNN7kZP6JHx/Lx1tOXyr6B3bk8/GW08SPjzU5mf8ip06m8SdPXSr6rtd30/iTp4icOtnkZP4LTZlA5ZNPXir6zbn7qHzySUJTJvTI9kRvGnGYnp4uA7kO/1z9xxQXP4LdHo3bfZqwsDE4HFE9kLBn1dToHCgYjdMpaG21MWhQJOHhg8yOFZDw8xpTT8fjdYCz3QYDnYRGhZsdKyBVNcf565GXCQkJp6W1kYGxcfQbaL2/XxcaQrnYnIlm19E9DiIHO+g/MNLsWAGxlxWTuOcFZL8ItOZ67AmjCI2JNjtWQLyNjbSVlhJ+6620FhcT9/zzhGfe7Nf3EEIUSCnTr/Z1lj/CBxgYfQt2ezRtbZU4nUMtWewBYmNt2O1uLl60ExamWbbYAzQP8NKKh9B2O95QYdliDzA8dixOWz+aW88TFh5lyWIP0D+qFWQrujcUR6i0bLEH8CRNxGMPx3GhDqIGWbbYA2iRkWCz0fzXvxL9wP1+F3t/9Ko7bQN1rv5j3O7ThITEoesXGTXqcQZG32J2LL8dLN6G251HSOhF3G4bU9JimDRxrtmxAuL68Aji2BmabW1obTZakh0k3matNojPiZ17ufj+BcIcUXha20i75R7GzJpmdiy/HdiRT+2WeuAC7a0aCck6U2Zf9aCwV3K9vptmdwNt/QZia2lB3rOQxAV3mh0rIM25+zi5dCn24cOpf+VVwjJu7rGib/kjfF/PPixsDP36jSAlZdVlPX2r8PXsw/tfpH+4ZO7dt17W07cSX8/+gqONthAd+13DLuvpW4mvZx8VFkNE+EBmL152WU/fKnw9e5vWhN3p5pb5Qy/r6VuJr2ffFjUc74ChRP7w2ct6+lbi69mHjB6NMy6OuOefv6yn390sX/CbGotJSVl1qY0zMPoWUlJW0dRYbHIy/1ScPMbcu28lxCkAmDRxLnPvvpWKk8dMTua/xhN12O8ahrfj34+Jt43HftcwGk/UmRssANWfHGH24mWEOMMAGDNrGrMXL6P6E2vtvCqO1HDL/KFodh2AKbPTuWX+UCqO1JiczH+N+4uI/OGzyH79AUhccKdR9PcXmZzMf60lh4h7/nmjrQOEZ95M3PPP01pyqEe2d12ctAUoOPANANKm/LE7I33hXnzxOwA8/vivTE4SvOJndgIwccUsk5MEb/OjjwPw4EsvmpwkOC898QYAj676mslJgldw+30ApH3wpslJglf+zYUAJPxhU0Dr96mTtoqiKMrVqYKvKIrSR6iCryiK0keogq8oitJHqILvh4aGA7hcq2loOGB2lKBVVFSQk5NDRUWF2VGC0lbeSOMHFbSVN5odJWhVx4+wb8ufqDpurSuAPk9NWQMFO1zUlDWYHSVoLYWF1L20hpbCQrOjBO26uPHqi9DQcIADhd9E193YbE6mpP6BqKgpZscKSEVFBRs3bsTr9aJpGosWLSI+Pt7sWH5rK2+kLvsQ0qMj7DYGL51ASII17x6tOn6E13/yA7weD5rdzoIf/ozhY615o1pNWQNbny/E69HR7DbmPZlKbJI1705uKSzk5ENLkG43wulk5Pp1hKWmmh0rYOoI/xrV1+9D192Ajq63U1/fMzdGfBFcLhderxcpJV6vF5fLZXakgLSVNSA9OkiQHp02Cx9NVhw+hNfjQeo6Xo+HisM9cx32F6HyeD1ej46U4PXqVB6vNztSwFry9iPdbtB1ZHs7LXn7zY4UFFXwr1F09M3YbE5Aw2ZzEB3dc8+76GmJiYlomoYQAk3TSExMNDtSQEKSohB2GwgQdhshFj2KBIhPnoBmtyNsNjS7nfjknnla4hchbmw0mt2GsIGm2Ygba93n3IRlTEU4naBpCIeDsIypZkcKStAtHSHEjcBrXRYlAT8CBgAPA7Udy78vpXwv2O2ZJSpqClNS/0B9/T6io2+2bDsHID4+nkWLFuFyuUhMTLRkOwcgJCGSwUsn0FbWQEhSlGXbOQDDx45nwQ9/RsXhQ8QnT7BsOwcgNimKeU+mUnm8nrix0ZZt5wCEpaYycv06WvL2E5Yx1dLtHOiGgi+lPAZMBhBCaEAlsAV4CHheSmm9qQT/QFTUFEsX+q7i4+MtW+i7CkmItHSh72r42PGWLvRdxSZFWbrQdxWWmmr5Qu/T3S2drwClUsrybv6+iqIoSpC6u+DfD7zS5fN/CCGKhRDrhBDWbeQpiqJcB7qt4AshnMA9wOsdi1YDozHaPdXA5w5nFUI8IoTIF0Lk19bWft6XKIqiKN2gO4/w5wAHpJSnAaSUp6WUXimlDqwFMj5vJSnlGillupQyPSYmphvjKIqiKF11541XD9ClnSOEGCalrO74OB8o6cZtXfL+B/OJHnD59KHCwl9Tf34vd9y+pSc22SNW/+JZhg+/fIe3dVM2VVW1ZH3vKZNSBebgU28hhoUA2qVlxau2I6vbmPTsveYFC8Dvl/w7sYljL1v21o9/Q43rOI+t+61Jqfy3celmBseFAOLSsndXvE5dZRuLsh80L1gA/jJzOlpyChFdlu1ZloX3cAlf2ZVjWq5A7Pjpo8RMmcbgLssK3ttI7YG9zH76pW7fXrcc4QshwoCZQNcHU/9KCHFICFEM3A482R3b+qzoAdM4e+73NDWdBIxif/bc76/YCfR2w4fHUFh6ipZGB2AU+8LSU1fsBKxADAshujIc50Vj1kLxqu1EV4Z37ASsJTZxLKWHd9Fw4RxgFPvSw7uu2An0doPjQnDVDES/aOyE313xOq6agR07AWvRklOI2bEHvbESMIp9zI49aMkpJifzX8yUaXif/iVN54xBNAXvbcT79C+JmdIz9eu6GIBSWPhrztX/HgApob09El0f0N3xety5svGU1g0xPkjoL0KJdPY3N1SA4hpCmCTGASClpFlvwutoNzlVYI7UfkJ5Y+dIQ5sWiTNsoImJAuNsGYI78jbjg5SEttcTprnNDRUgUbudSUeNWiGBugEaFweFmRsqQI4WN7HVbTRFhyDa3Gg//S5pX13k1/foUwNQUlP/z6X3uu6wZLEHGJjU+dAsh7BZttgDVEa1XXrfLt2WLfYA42Nu6vLJacliD+AOO3Ppvc3batliDyBj5lx6f9GJZYs9QHuYk7pIiKxvo252mt/F3h/XxcPTCgt/jZTQ3Dyc8PAqhg/76mU7AavYuikb5CmG0J8zXGDoiFDmLVxqdqyAFK/ajmyQ1LVVMThkOBdiQ5j4xJyrr9gLvfXj3wBg04aje6uIG3kD9/7oWyan8p/RxtGJaHLRFDGKyNgW7npmgdmxArJnWRY6cGK4jTFVOiE3pDHjhdVmxwpIwXsb8fzgF+yYEcmtOwoouG1jjxV9yx/h+3r2tWdmcO7sEgYNfIyz535PYeGvzY7mF1/PfpIcwW22iaSOHkFh6SljJ2Axvp79J437OHzxI+rjmomuDKd41Xazo/nN17N3hmUSNngWo5NnUnp416WdgFX4evbxp3aRXP0eibHncNUM5N0Vr1995V7G17N/e1oor94TSe3sGcTs2MOeZVlmR/Obr2f/+jcTKL73JrSffhfv07+k4L2NPbI9yxf8+vN7GTTwMVpbjd5kaur/YdDAx6g/v/cqa/YuVVW1pI4ewWhtJADzFi4ldfQIqqqsd2+CrG6jPq6ZWt04kT7xiTnUxzUjq9uusmbvU+M6zujkmTgjjPMR9/7oW4xOnkmN67jJyfxTV9lGYuw54lr+DsBdzywgMfYcdZXW+zPxHi6hdvYMiieGAjDjhdXUzp6B93CPXAjYo2oP7EX76Xc5e9MwANK+ugjtp9+l9kDP1K/r4qQtwPr16wF46KGHujPSF674mZ0ATFwxy+Qkwdv86OMAPPjSiyYnCd5LT7wBwKOrvmZykuAU3H4fAGkfvHmVr+z9HnjhFgBeWfaxyUmC99AOo26tn70+oPX71ElbRVEU5epUwVcURekjVMFXFEXpI1TBVxRF6SNUwfdDRUUFOTk5VFRUmB0laG3ljTR+UEFbeaPZUYJSdfwI+7b8iarjR67+xb1cTVkDBTtc1Fh4Nq9PS2EhdS+toaWw0OwoQSs6U0T2oWyKzhSZHSVo18WNV1+EiooKNm7ciNfrRdM0Fi1aZNmJUW3ljdRlH0J6dITdxuClEyw5Narq+BFe/8kP8Ho8aHY7C374M8tOjKopa2Dr84V4PTqa3ca8J1MtOzGqpbCQkw8tQbrdCKeTkevXWXZiVNGZIh7e+TBurxun5mTtrLVMHjLZ7FgBU0f418jlcuH1epFS4vV6cblcZkcKWFtZA9KjgwTp0Wmz6BFlxeFDeD0epK7j9XioOHzI7EgBqzxej9ejIyV4vTqVx+vNjhSwlrz9SLcbdB3Z3k5L3n6zIwUs/3Q+bq8bHZ12vZ3804FdNt5bqIJ/jRITE9E0DSEEmqaRmJhodqSAhSRFIew2ECDsNkIseiQZnzwBzW5H2GxodjvxyRPMjhSwuLHRaHYbwgaaZiNurHUHxIVlTEU4naBpCIeDsIypZkcKWPrQdJyaE01oOGwO0ode9VL3Xk21dK5RfHw8ixYtwuVykZiYaNl2DhiDvwcvnUBbWQMhSVGWbOeAMfR7wQ9/RsXhQ8QnT7BsOweMod/znkyl8ng9cWOjLdvOAWPo98j162jJ209YxlTLtnMAJg+ZzNpZa8k/nU/60HRLt3NAFXy/xMfHW7rQdxWSEGnZQt/V8LHjLV3ou4pNirJ0oe8qLDXV0oW+q8lDJlu+0Puolo6iKEofoQq+oihKH6EKvqIoSh+hCr6iKEof0W0nbYUQLqAJ8AIeKWW6EGIg8BqQCLiAf5FSdusFxhtWZZM0atRly/72zl8o+/RTFj9hnWlRB596q2PIt3ZpWfGq7cjqNiY9e695wQLw+yX/fsWQ77d+/BtqXMd5bN1vTUoVmI1LN3cM+haXlr274nXqKttYlP2gecH89JeZ09GSU4josmzPsiy8h0v4yq4c03IFYuHqVCZG3HjZspWb76e46Ribsqx1Z++6d5aQMuLygeV5hdmUnNrLkrvXdfv2uvsI/3Yp5eQuz2X+HvAXKeUNwF86PnerpFGjeL8gh/Nnjf3I3975C+8X5FyxE+jtxLAQoivDcV405hP4pkYZOwFriU0cS+nhXTRcOAd0To367E7ACgbHheCqGYh+0dgR+yZHGTsB69CSU4jZsQe9sRLonBqlJaeYnMx/EyNuZJOnhFrHBcAo9ps8JVfsBKwgZcQ0lhe9QOOFasAo9suLXrhiJ9Bdum0ASscRfrqUsq7LsmPADClltRBiGLBHSvkP/1QCHYDiK/IASIiwhRHhsN5Q47iGECYJY7KSlJJmvcmyw7+P1H5CeWPn1B6bFmnZ4d/OliG4I42JakhJaHu9JQeAi9rtTDpq/HxJoG6AZtnh3wXjm3hlQkeDQkrivDBMCzU3VIAa9Xb+bvMSKxy06u2snLyMjFT/uhNmDECRwE4hRIEQ4pGOZUOllNUAHa9DPifoI0KIfCFEfm1tYOP8vnT3Vy69dwjNksUeoDKqc9xcu3RbttgDjI+5qcsnp2WLPYA77Myl9zZvqyWLPYCM6Rwif9GJZYs9QNqRzuZUuC4tW+wBIm0OYr061Xj4lwEpfhd7f3TnjVe3SimrhBBDgF1CiKPXspKUcg2wBowj/EA2/Ld3/gISYm3R1Mh6xk0Yf9lOwCqKV21HNkjq2qoYHDKcC7EhTHxiztVX7IV8Q75t2nB0bxVxI2/g3h99y+RUgTHaODoRTS6aIkYRGdvCXc8sMDuW3/Ysy0IHTgy3MaZKJ+SGNGa8sNrsWAFZufl+hKeEia1uikOdJNvHsPzBV82OFZC8wmyWF/4Xj16U/EkvIaMwu8eKfrcd4UspqzpezwBbgAzgdEcrh47XM//4OwTG1865pV8y98Xezh1p03m/IMfYCViIr2f/SeM+Dl/8iPq4ZqIrwyletd3saH7z9eydYZmEDZ7F6OSZlB7edWknYCW+nn38qV0kV79HYuw5XDUDeXfF62ZH84uvZ//2tFBevSeS2tkziNmxhz3LssyO5jdfz/7h+maebvay0J7CJk8JKzffb3Y0v/l69ivbI/iPkARWTl7G8qIXyCvM7pHtdUvBF0KECyEifO+BWUAJ8DawqOPLFgFbu2N7XZV9+il3pE0nNdo4NfClu7/CHWnTKfv00+7eVI+S1W3UxzVTq58EYOITc6iPa0ZWt11lzd6nxnWc0ckzcUYY5yPu/dG3GJ08kxrXcZOT+a+uso3E2HPEtfwdgLueWUBi7DnqKq315+I9XELt7BkUTzRaHzNeWE3t7Bl4D5eYnMx/xU3HWGhPYaZulK/lD77KQnsKxU3HTE7mv5JTe42evXMwABmpS1k5eRklp/ZeZc3AdMtJWyFEEsZRPRhtoj9KKX8mhBgE/AkYCZwEFkgpz/2j7xPoSVuAMy8VAzDk0YkBrd9bbH70cQAefOlFk5ME76Un3gDg0VVfMzlJ8Apuvw+AtA/eNDlJcB544RYAXln2sclJgnf0udEAjPt2qclJusH6u4zXh94NaPVrPWnbLT18KWUZMOlzlp8FrNdMVxRFuQ6pO20VRVH6CFXwFUVR+ghV8BVFUfoIVfD90FbeSOMHFbSVN5odJWhVx4+wb8ufqDp+xOwoQakpa6Bgh4sai87l7aqlsJC6l9bQUmit58F8nqIzRWQfyqboTJHZUYJXkQc5zxmvFqcmXl2jtvJG6rIPIT06wm5j8NIJlp0YVXX8CK//5Ad4PR40u50FP/yZJadG1ZQ1sPX5QrweHc1uY96TqZadGNVSWMjJh5Yg3W6E08nI9essOzGq6EwRD+98GLfXjVNzsnbWWutOjKrIg433gNcNmhMWvQ3xGWanCpg6wr9GbWUNSI8OEqRHp83CR5QVhw/h9XiQuo7X46Hi8CGzIwWk8ng9Xo+OlOD16lQe79YHsX6hWvL2I91u0HVkezstefvNjhSw/NP5uL1udHTa9XbyTwd2qXWv4Moxir30Gq8uaz1Z9LNUwb9GIUlRCLsNBAi7jRCLHkkCxCdPQLPbETYbmt1OfPIEsyMFJG5sNJrdhrCBptmIGxttdqSAhWVMRTidoGkIh4OwjKlmRwpY+tB0nJoTTWg4bA7Sh1718vDeK3G6cWQvNOM1cbrZiYKiWjrXKCQhksFLJ9BW1kBIUpRl2zlgDP5e8MOfUXH4EPHJEyzZzgFj6Pe8J1OpPF5P3Nhoy7ZzwBj6PXL9Olry9hOWMdWy7Rwwhn6vnbWW/NP5pA9Nt247B4z2zaK3jSP7xOmWbueAKvh+CUmItHSh72r42PGWLfRdxSZFWbrQdxWWmmrpQt/V5CGTrV3ou4rPsHyh91EtHUVRlD5CFXxFUZQ+QhV8RVGUPkIVfEVRlD7C8idtT/xiN86kKELpHCp98k/7cZc1MOZ7d5qYzD+/X/LvVwz5fuvHv6HGdZzH1v3WpFSB2bh0c8eQb3Fp2bsrXqeuso1F2Q+aFywAf5k5HS05hYguy/Ysy8J7uISv7LLONdkLV6deMeR75eb7KW46xqYsa93ZW/zzVFpiJzGgy7LcdYsJqznIxO9b6/eSu+mH9B+dQddR8iUfvcOF0jwyF/6k27dn+SN8Z1IUouAibWeMxx2c/NN+RMFFnBa7ciM2cSylh3fRcMEYF+CbGvXZnYAVDI4LwVUzEP2iBnROjTJ2AtaiJacQs2MPemMl0Dk5SktOucqavcvEiBvZ5Cmh1nEB6Jwa9dmdgBW0xE4io3wLES0tgFHsM8q30BJ7xRPae73+ozOI2/VvXGisA4xiH7fr3+g/umeuCuqWASjdJdABKCf/tB/bgVYApJRc5ALtFhw0faT2E8obOyfd2LRIyw7/drYMwR15m/FBSkLb6y07/FvUbmfSUePvpQTqBmiWHABeML6JVyZ0/KNeSuK8WHb4981nJI+1nABASqgUMTQ4Y01OFZhQTxOjvC7O2YegeVupnPk7Um6926/vca0DUCx/hA8w8l8670r04LZksQcYH3NTl09OyxZ7AHdY5/him7fVssUeQMZ0DpK/6MSSxR4g7UhnYypcl5Yt9gD7hnS2C5sItWyxB2i1R1DFYAZ7z3BsxAK/i70/LN/Dh442jpSc89YwUIslYnLcZTsBq/AN+bZpw9G9VcSNvIF7f/Qtk1MFxmjj6EQ0uWiKGEVkbAt3PbPA7FgB2bMsCx04MdzGmCqdkBvSmPHCarNj+W3l5vsRnhImtropDnWSbB/D8gdfNTtWQHLXLUaXJzggk5giSmmKHUPmkg1mxwpIyUfvELYzi/XaAu499TolH93eY0U/6CN8IUS8EOIDIcQRIcRhIcS3Opb/pxCiUghR1PHfV4OPeyVfz/7vsoijEQeQaf0QBRc5+SdrPXzK17N3hmUSNngWo5NnUnp416WdgJX4evbxp3aRXP0eibHncNUM5N0Vr5sdzW++nv3b00J59Z5IamfPIGbHHvYsyzI7ml98PfuH65t5utnLQnsKmzwlrNx8v9nR/Obr2f9Ov5tfa4+RlzCfjPIt5K5bbHY0v/l69r/o/z12DPn/qJz5O+J2/RslH73TI9vrjpaOB/i2lHI8kAn8uxDC15t4Xko5ueO/97phW1dwlzUg0/pRG14FGO0dmdYPt8WeZlnjOs7o5Jk4I8YBcO+PvsXo5JnUuI6bnMx/dZVtJMaeI67l7wDc9cwCEmPPUVfZZnIy/3kPl1A7ewbFE432x4wXVlM7ewbewyUmJ/NPcdMxFtpTmKkbP/LLH3yVhfYUipuOmZzMf2E1B8lLmE+O3XiQWeaSDeQlzCes5qDJyfx3oTSPypm/ozwyDYCUW++mcubvuFDaM8/e7/aTtkKIrcD/BW4FLkgpV17ruoGetAV47ZnvAfD1Fb8IaP3e4qUn3gDg0VVfMzlJ8Apuvw+AtA/eNDlJ8B544RYAXln2sclJgnP0udEAjPt2qclJgvf1Z9YA8NqKR0xOEryvv2T8vXrt0VsCWt+Uk7ZCiEQgFdjXseg/hBDFQoh1QgjrPrtWURTlOtBtBV8I0R/4M7BMStkIrAZGA5OBauC5f7DeI0KIfCFEfm1tbXfFURRFUT6jWwq+EMKBUexfllK+CSClPC2l9EopdWAt8Ll3Ekgp10gp06WU6TExMd0RR1EURfkc3XGVjgD+GzgipfyvLsuHdfmy+YC1znIpiqJcZ7rjOvxbgW8Ch4QQvhH13wceEEJMxrg50QU82g3bMlXV8SOWnxLlU1PWcF1MimopLLwupkSBMfz7upgSBcbw7+tkSlRBeT25ZWfJTBpEWoK1T0UGXfCllB/S9SlZnXrkMkyzVB0/wus/+QFejwfNbmfBD39m2aJfU9bA1ucL8Xp0NLuNeU+mWrLotxQWcvKhJUi3G+F0MnL9OssW/aIzRTy882HcXjdOzcnaWWutW/Qr8mDjPcbQb81pjAi0aNEvKK/nX7NzcXt0nHYbLy/NtHTRvy4erfBFqDh8CK/Hg9R1vB4PFYcPmR0pYJXH6/F6dKQEr1en8ni92ZEC0pK3H+l2g64j29tpybPWzXZd5Z/Ox+11o6PTrreTfzqwy5N7BVeOUeyl13h1Weepop+VW3YWt0dHl9Du0cktO2t2pKCogn+N4pMnoNntCJsNzW4nPnmC2ZECFjc2Gs1uQ9hA02zEjbXmEUtYxlSE0wmahnA4CMuw3uM0fNKHpuPUnGhCw2FzkD70qpdU916J040je6EZr4nTzU4UsMykQTjtNjQBDruNzKRBZkcKynXxLJ0vwvCx41nww59dFz382KQo5j2ZavkeflhqKiPXr7sueviTh0xm7ay110cPPz7DaONcBz38tIRoXl6aqXr4fdHwseMtXei7ik2Ksmyh7yosNdXShb6ryUMmW7vQdxWfYelC31VaQrTlC72PaukoiqL0EargK4qi9BGq4CuKovQRquAriqL0EZY/abv+//8u8SmXDy/eve6PVJQc5KH/+qVJqfy3cenmjiHfnfewvbvideoq21iU/aB5wQLwl5nT0ZJTiOiybM+yLLyHS/jKLmtdk71wdeoVg75Xbr6f4qZjbMoqNCmV/4p/nkpL7CQGdFmWu24xYTUHmfh96/w+ADKf3ULqsMvHM2Zt2E5hdSu5T803KVVgFr5VxG2Jl1/quaaogg9dZ9l0b/efwLf8EX58yiQO/s8fqa8+DRjF/uD//PGKnUBvNzguBFfNQPSLGtA5NcrYCViLlpxCzI496I2VQOfUKC05xeRk/psYcSObPCXUOi5LUDZMAAAgAElEQVQAnZOjPrsT6O1aYieRUb6FiJYWoHNqVEustX5OAFKHhbL9qE6pOxIwiv32o/oVOwEruC1xED9/s4SqC8ZwoDVFFfz8zZIrdgLdpdsHoAQj0AEoviLvY3NE4wy13iWHzpYhuCNvMz5ISWh7vWWHf4va7Uw6avxZSqBugGbZ4d8F45t4ZULHP4alJM6LJQeA33xG8ljLCQCkhEoRY9nh36vlnWxrMm60k4AWaiMswnoHRwBtbR7cje2EhNlxe3S+f18Kj0yO9+t7mDIAxSx3LvlG5wfhtGSxB3CHnbn03uZttWyxB5Axcy69v+jEssUeIO1IZ3MqXJeWLPYA+4Z0tgubCLVssQfIErsvvRcali32ACEhdgi14W7xkJYyxO9i7w/L9/DBOMIHsIeOwNN6ihunTbt8J2ARRhtHJ6LJRVPEKCJjW7jrmQVmxwrInmVZ6MCJ4TbGVOmE3JDGjBdWmx0rICs334/wlDCx1U1xqJNk+xiWP/iq2bH8lrtuMbo8wQGZxBRRSlPsGDKXbDA7VkCyNmxHNurYIm3IRp3pg9pYvXjO1VfshdYUVfCzP5fQb2wUBSVnWDO2oseKvuWP8H3tnMihMxk58REm/dM3OPg/f7y0E7AKX88+/tQukqvfIzH2HK6agby74nWzo/nN17N/e1oor94TSe3sGcTs2MOeZVlmR/Obr2f/cH0zTzd7WWhPYZOnhJWb7zc7ml98Pfvf6Xfza+0x8hLmk1G+hdx1i82O5jdfz94xKhTnpCHMGWdj+1GdrA3bzY7mN1/PfnhmLONSh/L9+1L4+ZslrCmq6JHtWb7gV5QcZNI/fYPBCcYDmu5c8g0m/dM3qCix1gT7uso2EmPPEdfydwDuemYBibHnqKtsMzmZ/7yHS6idPYPiiUbrY8YLq6mdPQPvYevNwCluOsZCewozdeNHZfmDr7LQnkJx0zGTk/knrOYgeQnzybEbPyeZSzaQlzCfsBpr/ZwAFFa3MmecDW2Ecc3R6sVzmDPORmF1q8nJ/Peh6yzfvy+F2OH9AXhkcjzfvy+FD10981TO6+KkLcCW5w4AMP/bU7oz0heu4Pb7AEj74E2TkwTvgRduAeCVZR+bnCR4R58bDcC4b5eanCQ4X39mDQCvrXjE5CTBu/HdAgCO3ZVmcpLgzS80DvS2pN4Q0Pp96qStoiiKcnWq4CuKovQRPV7whRCzhRDHhBAnhBDf6+ntKYqiKJ+vRwu+EEIDfgvMAW7CGGx+U09uU1EURfl8PX2EnwGckFKWSSndwKvAvB7eZo+pKWugYIeLmrIGs6MEraWwkLqX1tBSaK3nqHxW0Zkisg9lU3SmyOwowavIg5znjFeLKyiv57cfnKCg3JrzkrvKb2hmVflp8huazY4StJ6+8SoO6HpB6Sng5h7eZo+oKWtg6/OFeD06mt3GvCdTLTsxqqWwkJMPLUG63Qink5Hr11lyalTRmSIe3vkwbq8bp+Zk7ay11p0YVZEHG+8xhn5rTmNEoEUnRhWU1/Ov2bm4PTpOu42Xl2ZadmJUfkMzXys6QbsucdgEb0weQ3pUuNmxAtbTR/jic5Zddh2oEOIRIUS+ECK/tra2h+MErvJ4PV6PjpTg9epUHrfukUtL3n6k2w26jmxvpyVvv9mRApJ/Oh+3142OTrveTv7pwC7p7RVcOUaxl17j1WWtp4p2lVt2FrdHR5fQ7tHJLeuZa8q/CHvPX6Bdl3iBdl2y9/wFsyMFpacL/img6z3CI4Cqrl8gpVwjpUyXUqbHxMT0cJzAxY2NRrPbEDbQNBtxY615xAIQljEV4XSCpiEcDsIyppodKSDpQ9Nxak40oeGwOUgfetXLkHuvxOnGkb3QjNfE6WYnClhm0iCcdhuaAIfdRmZSzzz58YswbUB/HDaBBjhsgmkD+psdKSg93dLZD9wghBgFVAL3A9Z7yA3G0O95T6ZSebyeuLHRlm3ngDH4e+T6dbTk7ScsY6ol2zlgDP1eO2st+afzSR+abt12Dhjtm0VvG0f2idMt284BY+j3y0szyS07S2bSIMu2cwDSo8J5Y/IY9p6/wLQB/S3dzoEeLvhSSo8Q4j+A/wE0YJ2U8nBPbrMnxSZFWbrQdxWWmmrZQt/V5CGTrV3ou4rPsHSh7yotIdrShb6r9Khwyxd6nx5/WqaU8j3gvZ7ejqIoivK/U3faKoqi9BGq4CuKovQRquAriqL0EZafePXmd94hbuwAoPOkyr7sHCqPn+e+X91tXjA//WXmdLTkFCK6LNuzLAvv4RK+ssta12QvXJ16xZDvlZvvp7jpGJuyrHVnb/HPU2mJncSALsty1y0mrOYgE79vnd9L5rNbrhjynbVhO4XVreQ+Nd+kVIHJ3L2DlFAn0HkBxdIP36ek1U3unbPNCxaAnx58m7ToocDAS8u2n9xHQf1pnp50T7dvz/JH+HFjB5C/301j9XnAKPb5+90dOwHr0JJTiNmxB72xEuicGqUlp5iczH8TI25kk6eEWodxk4pvatRndwJW0BI7iYzyLUS0tACdk6NaYieZnMw/qcNC2X5Up9QdCXROjfrsTsAKUkKdbHMPoM1uzHxe+uH7bHMP6NgJWEta9FCeKHVzrrURMIr9E6Xujp1A97suBqDsy84hP7/d+CAl/TwN9LNZb/qNqN3OpKPG718CdQM0yw7/LhjfxCsTOv4BKSVxXiw7/PvmM5LHWk4AICVUihhLDgBfLe9kW5Nxk50EtFCbZYd/hw6r59QwYygNUhIj6ojVLpobKkBNuh2XHMYQ0cRFKVg12smckf49gaZPDUC5eWnnXYma3mbJYg8gYzqHMF90YtliD5B2pLM5Fa5LyxZ7gH1DOp8Q0kSoJYs9QJbYfem90LBssQdore68xj+Ui5Yt9gARNg+DqeOMjOSfo876Xez9YfkePhhH+EjJgIsVnO8XT9LU+Mt2AlaxZ1kWOnBiuI0xVTohN6Qx44XVZscKyMrN9yM8JUxsdVMc6iTZPoblD75qdqyA5K5bjC5PcEAmMUWU0hQ7hswlG8yO5besDduRjTq2SBuyUWf6oDZWL55z9RV7oaUfvs+7bp1R8gRlYgwJtv5k33aH2bECYrRx2viabQd/briVL5/c12NF3/JH+L6e/bjGHG6T75M+1Un+frexE7AQX8/+7WmhvHpPJLWzZxCzYw97lmWZHc1vvp79w/XNPN3sZaE9hU2eElZuvt/saH7z9ex/p9/Nr7XHyEuYT0b5FnLXLTY7ml98PXvHqFCck4YwZ5yN7Ud1sjZsNzua33w9+7n62zzm+CNznefZ5h7A0g/fNzua33w9+++E/ImHIopZNdrJE6Vutp/c1yPbs3zBrzx+nvSpTsbYywCjvZM+1Unl8fMmJ/OP93AJtbNnUDzRaH3MeGE1tbNn4D1cYnIy/xU3HWOhPYWZuvHXa/mDr7LQnkJx0zGTk/kvrOYgeQnzybEb/2LMXLKBvIT5hNUcNDmZfwqrW5kzzoY2wriYYfXiOcwZZ6Ow2nrtz5JWN3Od55nuNGYgZN92B3Od5ylpdZuczH8F9adZNdrJlNBzAMwZeTOrRjspqD/dI9u7Lk7aApR/cyEACX/Y1J2RvnAPvHALAK8s+9jkJME7+pxxUm3ct0tNThK8rz+zBoDXVjxicpLg3PhuAQDH7kozOUnwNv3tXgAWfuktk5MEr+CA8UzJtCl/DGj9PnXSVlEURbk6VfAVRVH6CFXwFUVR+ghV8BVFUfoIVfD90FJYSN1La2gptM4zVP6RojNFZB/KpuhMkdlRglORBznPGa8WV1Bez28/OEFBuXXnJfvkNzSzqvw0+Q3NZkcJWkPDAVyu1TQ0HDA7StCuixuvvggthYWcfGgJ0u1GOJ2MXL/OshOjis4U8fDOh3F73Tg1J2tnrbXm1KiKPNh4jzH0W3MaIwItOjGqoLyef83Oxe3RcdptvLw007ITo/Ibmvla0QnadYnDJnhj8hjLToxqaDjAgcJvoutubDYnU1L/QFTUFLNjBUwd4V+jlrz9SLcbdB3Z3k5L3n6zIwUs/3Q+bq8bHZ12vZ3804FdCms6V45R7KXXeHVZ62a7rnLLzuL26OgS2j06uWVnzY4UsL3nL9CuS7xAuy7Ze/6C2ZECVl+/D113Azq63k59fc/cEPVFUQX/GoVlTEU4naBpCIeDsIypZkcKWPrQdJyaE01oOGwO0ode9fLd3ilxunFkLzTjNdF6j9PwyUwahNNuQxPgsNvITBpkdqSATRvQH4dNoAEOm2DagP5mRwpYdPTN2GxOQMNmcxAd3XPPufkiBNXSEUL8GrgbcAOlwENSyvNCiETgCOC7tTJXSvlYMNsyW1hqKiPXr6Mlbz9hGVMt284BY/D32llryT+dT/rQdGu2c8Bo3yx62ziyT5xu2XYOGEO/X16aSW7ZWTKTBlm2nQPG0O83Jo9h7/kLTBvQ37LtHICoqClMSf0D9fX7iI6+2dLtHAi+h78LeEpK6RFC/BJ4Cvhux6+VSiktWkk+X1hqqqULfVeTh0y2bqHvKj7D0oW+q7SEaEsX+q7So8ItXei7ioqaYvlC7xNUS0dKuVNK6en4mAuMCD6SoiiK0hO6s4e/BOj66L1RQohCIcRfhRDWba4qiqJcJ67a0hFC7AY+b+LDD6SUWzu+5geAB3i549eqgZFSyrNCiDTgLSFEspSy8XO+/yPAIwAjR44M7HehKIqiXNVVC76U8s7/7deFEIuAucBXZMejN6WUbUBbx/sCIUQpMBa44vo/KeUaYA0YT8v09zegKIqiXJtgr9KZjXGS9stSypYuy2OAc1JKrxAiCbgBKAsq6T+w94Gv0i8zs8vMdyj8zY+5mJvLtFfe64lN9oiFq1OvGPK9cvP9FDcdY1OWte7sLf55Ki2xk+g6Rj533WLCag4y8fvW+r1kPrvlikHfWRu2U1jdSu5T801K5b/M3Ts6hnxHXVq29MP3KWl1k3vnbPOCBeD1HTNxhl8+RH5rznLczQdZMHuXSakCs+evTzEk5vKLDj75ZAtnavOY8eVnu317wfbw/y8QAewSQhQJIX7fsfxLQLEQ4iDwBvCYlPJckNv6XP0yM3GufoWmShdgFHvn6lfol5nZE5vrMRMjbmSTp4Rah3GTim9q1Gd3AlbQEjuJjPItRLQYxwC+qVEtsZOusmbvkzoslO1HdUrdkUDn5KjP7gR6u5RQJ9vcA2izG0NCfFOjjJ2AtTjDJxHu3kJ4WzVgFPtw95YrdgJWMCQmg5MVT9PcXAsYxf5kxdNX7AS6y3UxAKXwNz8mdPUrAEjg3EAHzdHW+oEEKBjfxCsTOv7RJSVxXiw7/PvmM5LHWk4AICVUihjLDv9eLe9kW5Nxo50EtFCbJQeAhw6r59QwYygNUhIj6iw7/Pvmpp3c3s+4RkRKaHNHghxwlbV6qxZCQurQvVHoso2R8T/lppv8+9djnxqAkvqtH116fzFEWLLYA6Qdibj0PlyXli32APuGiEvvmwi1bLEHyBK7L70XGpYs9gCt1Z3X+Idy0bLFHmBfxKxL771eh4WLPUAYbW3haPYGHI47/S72/rguHp5W+Jsf4wTKEpwklbsZMmvuZTsBq1i5+X6Ep4SJrW6KQ50k28ew/MFXzY4VkNx1i9HlCQ7IJKaIUppix5C5ZIPZsQKStWE7slHHFmlDNupMH9TG6sVzzI7lt6Ufvs+7bp1R8gRlYgwJtv5k33aH2bECsjVnOVJCQ9MQoiLO0NovjXnTV5odKyBGG+cHnK6ZysBBu/nkky09VvQtf4Tv69nvmDuUN5+cgjvrAZyrX6HwNz82O5pffD37h+ubebrZy0J7Cps8JazcfL/Z0fzm69n/Tr+bX2uPkZcwn4zyLeSuW2x2NL/5evaOUaE4Jw1hzjgb24/qZG3YfvWVexFfz36u/jaPOf7IXOd5trkHsPTD982O5jdfz951KpWqiuk0O+cT7t7C1pzlZkfzm69nX1M9j7a2f2Jk/E85WfE0n3yypUe2Z/mCfzE3F3fWA5TcOQow2jvurAe4mJtrcjL/FDcdY6E9hZm68Uey/MFXWWhPobjp2FXW7H3Cag6SlzCfHLtxv13mkg3kJcwnrOagycn8V1jdypxxNrQRRstg9eI5zBlno7C61eRk/ilpdTPXeZ7pTmP+QfZtdzDXeZ6SVrfJyfznbj5Is3M+FxuTAJg3fSXNzvm4m6339+tMbR4j43+KlGMBuOmm+YyM/ylnantmvsN1cdIW4KEdDwGwfvb67oz0hTv6nHFSbdy3S01OEryvP7MGgNdWPGJykuDd+G4BAMfuSjM5SXA2/e1eABZ+6S2TkwTvxRe/A8Djj//K5CTBW7/eqFsPPfRQQOv3qZO2iqIoytWpgq8oitJHqIKvKIrSR6iCryiK0keogu+HojNFZB/KpuhMkdlRgleRBznPGa8WVlBez28/OEFBeb3ZUYKW39DMqvLT5Dc0mx0laA0NB3C5VtPQcMDsKEGrqKggJyeHiooKs6ME7bq48eqLUHSmiId3Pozb68apOVk7a611J0ZV5MHGe4zB35rTGBNowalRBeX1/Gt2Lm6PjtNu4+WlmZadGJXf0MzXik7QrkscNsEbk8dYdmJUQ8MBDhR+E113Y7M5mZL6B8tOjKqoqGDjxo14vV40TWPRokXEx8ebHStg6gj/GuWfzsftdaOj0663k386sMtHewVXjlHspdd4deWYnSgguWVncXt0dAntHp3csrNmRwrY3vMXaNclXqBdl+w9f8HsSAGrr9+HrrsBHV1vp75+n9mRAuZyufB6vUgp8Xq9uFwusyMFRRX8a5Q+NB2n5kQTGg6bg/ShV73ktfdKnG4c2QvNeE205kCyzKRBOO02NAEOu43MpEFmRwrYtAH9cdgEGuCwCaYN6G92pIBFR9+MzeYENGw2B9HRN5sdKWCJiYlomoYQAk3TSExMNDtSUFRL5xpNHjKZtbPWkn86n/Sh6dZt54DRvln0tnFknzjdku0cMIZ+v7w0k9yys2QmDbJsOweMod9vTB7D3vMXmDagv2XbOWAM/Z6S+gfq6/cRHX2zZds5APHx8SxatAiXy0ViYqKl2zmgCr5fJg+ZbO1C31V8hmULfVdpCdGWLvRdpUeFW7rQdxUVNcXShb6r+Ph4yxd6H9XSURRF6SNUwVcURekjVMFXFEXpI4Iq+EKI/xRCVHbMsy0SQny1y689JYQ4IYQ4JoT4p+CjKoqiKMHojpO2z0spLxs1I4S4CbgfSAaGA7uFEGOllN5u2N5lsv40h1uGXX7ycVPOCj6uzmP1v1hnSEXxz1NpiZ1E10FtuesWE1ZzkInfLzQtVyAyn91yxZDvrA3bKaxuJfepnhvf1hMyd+/oGPQddWnZ0g/fp6TVTe6ds80L5qfXd8y8Ysj31pzluJsPsmD2LpNSBWb1L55l+PCYy5Zt3ZRNVVUtWd97yqRUgfnLpncZMXrkZcuOfXSIU6Un+crCu7p9ez3V0pkHvCqlbJNSfgqcAHrkkpBbhmWwsvTP1DS4AKPYryz98xU7gd6uJXYSGeVbiGhpATqnRrXETrrKmr1P6rBQth/VKXVHAp1Toz67E7CClFAn29wDaLMbg0J8k6OMnYB1OMMnEe7eQnhbNdA5NeqzOwErGD48hsLSU7Q0OgCj2BeWnrpiJ2AFI0aPZMuut2lpNB6nceyjQ2zZ9fYVO4HuEtQAFCHEfwKLgUYgH/i2lLJeCPF/gVwp5eaOr/tvYLuU8o3/7fsFOgBlU84Kfl32pvFBSkboNmJt1vqBBLj5jOSxlhMASAmVIsayw79XyzvZ1jQVAAlooTbLDv8OHVbPqWHGYBqkJEbUWXIA+M1NO7m9n/GvXimhzR1p2eHf58rGU1o3xPggob8IJdJpzZvVWj1uznkb6W/vh8frYf7Me7jx1gl+fY9uG4AihNgthCj5nP/mAauB0cBkoBp4zrfa53yrz92zCCEeEULkCyHya2trrxbncy2c/syl9/0lliz2APuGdP5vayLUssUeIEvsvvReaFi22AO0Vnde5x/KRUsWe4B9EbMuvfd6HZYt9gADk45ceu8QNssWe4BQu5NwQrjgvcikETf5Xez9cdUevpTyzmv5RkKItcC2jo+ngK53KowAqv7B918DrAHjCP9atvVZm3JWIKQk1QOFdrh99F2X7QSsInfdYnR5ggMyiSmilKbYMWQu2WB2rIBkbdiObNSxRdqQjTrTB7WxevEcs2MFZOmH7/OuW2eUPEGZGEOCrT/Zt91hdiy/bc1ZjpTQ0DSEqIgztPZLY970lVdfsRfauikb5CmG0J8zXGDoiFDmLVxqdqyAHPvoEFt2vs0UbQwHT31C0kdjeqzoB3uVzrAuH+cDJR3v3wbuF0KECCFGATcAPfIcXl/Pfnl7KBu1BJaP/mdWlv6ZTTkremJzPcbXs/+dfje/1h4jL2E+GeVbyF232OxofvP17B2jQnFOGsKccTa2H9XJ2mCdk+g+vp79XP1tHnP8kbnO82xzD2Dph++bHc0vvp6961QqVRXTaXbOJ9y9ha05y82O5jdfz36SHMFttomkjh5BYekpYydgMb6e/cz+U8kcMpH5M+9hy663OfbRoR7ZXrAnbX8lhDgkhCgGbgeeBJBSHgb+BHwC7AD+vSeu0AH4uDqP5aP/mYX2oYDR3lk++p/5uNpaz3kPqzlIXsJ8cuzGg8wyl2wgL2E+YTUHTU7mv8LqVuaMs6GNMFoGqxfPYc44G4XVrSYn819Jq5u5zvNMdxozELJvu4O5zvOUtLpNTuYfd/NBmp3zudiYBMC86Stpds7H3Wy9v19VVbWkjh7BaM04sTlv4VJSR4+gqiqwlrCZTpWeZP7Me4iPNNq3N946gfkz7+FU6cke2V5QJ227W6AnbQFY33EJ00Pvdl8gE3z9mTUAvLbiEZOTBO/GdwsAOHZXmslJgrfpb/cCsPBLb5mcJDgvvvgdAB5//FcmJwle8TM7AZi4YtZVvrL3O/NSMQBDHp0Y0PrddtJWURRFuT6ogq8oitJHqIKvKIrSR6iCryiK0keogu+PijzIec54tbiC8np++8EJCsrrzY4SlPyGZlaVnya/odnsKEFraDiAy7WahoYDZkcJWkVFBTk5OVRUVJgdJWht5Y00flBBW3mj2VGCpiZeXauKPNh4jzH0W3MaIwItOjGqoLyef83Oxe3RcdptvLw005JTo/Ibmvla0QnadYnDJnhj8hjLToxqaDjAgcJvoutubDYnU1L/YNmJURUVFWzcuBGv14umaSxatMiyE6Payhupyz6E9OgIu43BSycQkhBpdqyAqSP8a+XKMYq99BqvrhyzEwUst+wsbo+OLqHdo5NbdtbsSAHZe/4C7brEC7Trkr3nL5gdKWD19fvQdTego+vt1NfvMztSwFwuF16vFyklXq8Xl8tldqSAtZU1ID06SJAenbayBrMjBUUV/GuVON04shea8Zo43exEActMGoTTbkMT4LDbyEwaZHakgEwb0B+HTaABDptg2gDrPk8lOvpmbDYnoGGzOYiOvtnsSAFLTExE0zSEEGiaRmJiotmRAhaSFIWw20CAsNsISYq6+kq9mGrpXKv4DKON48oxir1F2zlgDP5+eWkmuWVnyUwaZMl2DhhDv9+YPIa95y8wbUB/y7ZzwBj6PSX1D9TX7yM6+mbLtnPAGPq9aNEiXC4XiYmJlm3nAIQkRDJ46QTayhoISYqydDsHVMH3T3yGpQt9V2kJ0ZYt9F2lR4VbutB3FRU1xdKFvqv4+HhLF/quQhIiLV/ofVRLR1EUpY9QBV9RFKWPUAVfURSlj1AFX1EUpY+w/Enb8lVfRR/1ZUZ1WfbpO7/E9ulfSXjiPdNy+Svz2S1XDPnO2rCdwupWcp+ab1KqwGTu3tEx5LvzEralH75PSaub3DtnmxcsAK/vmHnFoO+tOctxNx9kwexdJqXy3+pfPHvFkO+tm7Kpqqol63tPmZQqMAefegsxLATQLi0rXrUdWd3GpGfvNS9YAD753Xv0HxtLaJdSfHL3AS4cr+Gmf/tqt2/P8kf4+qgvk1DwLBfPlgNGsU8oeBZ91JdNTuaf1GGhbD+qU+o2rgbwTY367E7AClJCnWxzD6DNbgwJ8U2NMnYC1uIMn0S4ewvhbdVA5+Soz+4Eervhw2MoLD1FS6MD6Jwa9dmdgBWIYSFEV4bjvGjM8ihetZ3oyvCOnYC19B8bS/vOOtrqmwCj2LfvrKP/2J6ZZ31dDED59J1fMqrg5wBICVW2oZx3DOnueD1utbyTbU1TAWPiuxZqs+zw79Bh9ZwaNtr4ICUxos6yw79vbtrJ7f2M8YxSQps70pIDwM+Vjae0ruPnQkJ/EWrZ4d9xDSFMEuMAkFLSrDfhdbSbnCowNq+N/gzAbWsFLzhmDWbknf5dntunBqCMuvu7l95fEGGWLPYAWWL3pfdCw7LFHqC1uvMa/1AuWrbYA+yL6Jyo5PU6LFnsAQYmHbn03iFsli32AJVRbZfet0u3ZYs9gK7ptHgbCZH9uDi8ze9i7w/L9/Cho40jBcW2cUyURwmbcPdlOwGryNqwHdmoY4u0IRt1pg9qY/XiOWbHCsjSD9/nXbfOKHmCMjGGBFt/sm+7w+xYAdmasxwpoaFpCFERZ2jtl8a86SvNjuW3rZuyQZ5iCP05wwWGjghl3sKlZscKSPGq7cgGSV1bFYNDhnMhNoSJT1jzZ+Xk7gO0/08tf28rZGTVOE7uPtBjRT+oI3whxGtCiKKO/1xCiKKO5YlCiItdfu333RP3Sr6e/dp+S3g29nnK054ioeBZPn3nlz21yR7h69k7RoXinDSEOeNsbD+qk7Vhu9nR/Obr2c/V3+Yxxx+Z6zzPNvcAln74vtnR/Obr2btOpVJVMZ1m53zC3VvYmrPc7Gh+8fXsJ8kR3GabSOroERSWnjJ2Ahbj69l/0riPwxc/oplu2IIAAA1hSURBVD6umejKcIpXWe9nxdezPyz2cWZAFY5Zg2nfWcfJ3T3ziOygCr6U8utSyslSysnAn4E3u/xyqe/XpJSPBZXyf2H79K+Upz3F+9ELAKO9U572FLZP/9pTm+wRhdWtzBlnQxthtAtWL57DnHE2CqtbTU7mv5JWN3Od55nuLAIg+7Y7mOs8T0mr2+Rk/nM3H6TZOZ+LjUkAzJu+kmbnfNzNB01O5p+qqlpSR49gtDYSgHkLl5I6egRVVbUmJ/OfrG6jPq6ZWv0kABOfmEN9XDOyuu0qa/Y+F47X4Jg1mIvhLQCMvHMKjlmDuXC8pke21y0nbYUQAjgJ3CGl/LsQIhHYJqVM8ef7BHrSFuDrL30MwGuP3hLQ+r3Fje8WAHDsrjSTkwRv09+MS+QWfuktk5ME78UXvwPA44//yuQkwSl+ZicAE1fMuspX9n6bH30cgAdfetHkJMF77ZnvAfD1Fb8IaP0v+qTtdOC0lPLvXZaNEkIUCiH+KoT4h88SFkI8IoTIF0Lk19Za72hDURTFKq560lYIsRv4vItCfyCl3Nrx/gHglS6/Vg2MlFKeFUKkAW8JIZKllFfMCJNSrgHWgHGE7+9vQFEURbk2Vy34Uso7/7dfF0LYgfuASz0IKWUb0NbxvkAIUQqMBQLr1yiKoihB646Wzp3A/2vv3mObOs84jn8fOxdYgMTlDg2EAOEO5bI2hNB2KutGYYKppUK0K2yDqUL7A7FNy0CVKlVFbFOrDm10Gl0H7cqqlq20gqHeVLrQFCiXhotI0jg4NQ2XQh1TUhqS+N0fPglOggEbk+ODn49k5fj1Mf5x8vo5J6+Pz1thjDnR2iAifUXEbS3nAyOBmgS8llJKqTglouAvpP1wDsDdwCERKQe2AI8bY75KwGvZan9tgL98UM3+2oDdUW7YvmAD62pPsy/YYHeUGxIMHsDne55g8OacxtaV/H4/paWl+P1+u6PcsMba85z/wE9jbadRXMepqzrGnjdeo67q2LVXTnI3/MUrY8ySK7T9m/BpmreM/bUBHnlhN5eaQ2SkuXhlaaFjZ4zaF2zgoU+raQoZ0l3CljtGOHLWqGDwAAcO/oRQ6BIuVwZTJr/s2Bmj/H4/mzZtoqWlBbfbzeLFix07Y1Rj7XnOvnAY0xxC0lz0WTrBsTNG1VUd4/WnVtPS3Iw7LY0FTzzNoIIxdseK2y1xaYWusLvmHJeaQ4QMNDWH2F1zzu5IcSurv0BTyNACNIUMZfUX7I4Ul0BgD6HQJSBEKNREILDH7khx8/l8tLS0YIyhpaUFn89nd6S4NdYEMc0hMGCaQzTWBO2OFDf/0cO0NDdjQiFampvxHz1sd6QbogX/OhXm9yYjzYVbID3NRWF+b7sjxa0opwfpLsENpLuEohxnXlPF47kLlysDcONypePx3GV3pLjl5eXhdrsREdxuN3l5eXZHiltmfjaS5gIBSXORmZ997SclqdxxE3CnpSEuF+60NHLHTbA70g25Ja6l0xWmDvXwytJCdtecozC/t2OHcyA88feWO0ZQVn+BopwejhzOgfCk31Mmv0wgsAeP5y7HDudAeNLvxYsX4/P5yMvLc+xwDoQn/e6zdAKNNUEy87MdO5wDMKhgDAueeBr/0cPkjpvg6OEc0IIfk6lDPY4u9JGmZWc5ttBHys6e4uhCHyk3N9fRhT5S5tBeji70kQYVjHF8oW+lQzpKKZUitOArpVSK0IKvlFIpQgu+UkqlCMd/aLvkH3uZMaL9KZIbSr18VH2OjT+906ZUsSv+4DDFOe0/RC05WMOu+gZ2fc9Zp4Jt/fhRcjxF7dp2VqynPlDG/On/tClVfDZsWMOw/Lx2be+9v5njNT6WLVtlT6g4VD1TSubw9h+i1m4tp9F7noJfRb2YbVJ67Te/JXdc+0nkP964Gf/Rch7+o7MmPtr75hYGDC9o1/b5kUOc8lZx57yHEv56jj/CnzGiN2u2V1BXH54zdUOplzXbKzrtBJJdcU4WGwNBvnGH75ccrGFjINhpJ+AEOZ4imuqepXtz+HLXOyvW01T3bKedgBMMy89jV2klDd8IEC72u0orO+0Ekl3m8F64dgfJuBR+y9duLce1O9hpJ+AEueMmUbZjMxcu1gPhYl+2Y3OnnYATDBhewLbn1vJtQ/jLj58fOcS259Z22gkkSkImQEmUeCdA2VDq5entFW33cz3dGZTTPZHRusRnWVA32MptDP3T3OT3dN7/A2DUhdf5ftPfATCApN9OTtZge0PFqaqqG5UVQ6x7hl69uuPx9Lc1UzwGnslifGAgAAaD6emmW5+eNqeKzxHvhxysfbvtfi9PX3oNdN7vBODbhgucrfXRs09fmhsbmbuihCHjJ8b0b3T1BCi2WjZzeNtyz0y3I4s9wMiI65hlII4t9gCVPRa0LbdIlmOLPUBBweVpJtPTxZHFHuBkv8sdLOTGscUeYPzwe9qW09O7ObbYA3TL6kHPPn35+uyXTLr/gZiLfSwcP4YP4SN8Ab6b5+ETX4BZY/u12wk4RcnBGj4MBBlsXHwhIUbhZu3kfLtjxWVnxXqaAkJD5kSyGg/x9XeKuXf0crtjxeW99zcDlfTrl8aZM83kDslk1n2L7I4Vs9qt5RiCXMwxdK8XLvY3DJ3vvGEQCA/jAAzoP5xTp70MGTqR6Uuc9zuBy8M4hQ8upPyd/5I7duJNK/qOP8JvHbNfNWc0rz1exKo5o1mzvYINpV67o8Wkdcx+iSeb/fdNYoknm42BICUHnTeNQOuYffqglcyf8R/SB62kqe5ZdlastztazFrH7ItnjmL58iconjmKXaWV1k7AOVrH7EOF2RSU3EOoMBvX7iC1W501GTtcHrMvmr2IR9b9iaLZiyjbsbltJ+AkrcV+7ooSZjz8KHNXlLDtubV8fuTQTXk9xxf8j6rPsWrO6LYj+mUzh7Nqzmg+qnbW1Sx31TewxJPddkS/dnI+SzzZ7Kp33vXq6wNlpA9a2XZEf+/o5aQPWkl9oMzmZLE7XuOjeOaotiP6WfctonjmKI7X+OwNFqNG73lChdltR/RD508iVJhNo9d516v3Hy2naPaitiP66UsWUTR7Ef6jztt5nfJWtRuzHzJ+InNXlHDKW3VTXu+W+NBWKaVSWUp9aKuUUuratOArpVSK0IKvlFIpQgu+UkqlCC34SimVIpLqLB0R+RKotTvHNfQBztod4jpozsRzSlbNmXjJnnWoMabvtVZKqoLvBCKy73pOf7Kb5kw8p2TVnInnpKxXo0M6SimVIrTgK6VUitCCH7u/2R3gOmnOxHNKVs2ZeE7KGpWO4SulVIrQI3yllEoRWvCjEJEFInJUREIiMi2iPU9ELorIp9btrxGPTRWRwyJSLSLrRETszGo99jsrT6WI/CCi/YdWW7WIlHRFzg65nhSRLyK24wPXymwXu7fV1YiIz+pzn4rIPqvtNhF5V0Q+s356bMr2ooicEZEjEW1XzCZh66xtfEhEptic0zH9MybGGL1d4QaMAUYBO4FpEe15wJEoz9kLTAcE2AHMtjnrWKAcyASGAV7Abd28QD6QYa0ztou375PAr6/QfsXMNvYD27fVNfL5gD4d2v4AlFjLJcDvbcp2NzAl8v0SLRvwgPWeEaAQ2GNzTkf0z1hveoQfhTHmmDGm8nrXF5GBQC9jzMcm3DNeAubftIARrpJ1HvCqMabRGHMcqAbutG7VxpgaY8wl4FVr3WQQLbNdknlbRTMP2GQtb6KL+mFHxpj/AV91aI6WbR7wkgnbDeRY7ym7ckaTbP0zJlrw4zNMRA6KyIciMtNqGwyciFjnhNVmp8GAP+J+a6Zo7V3tl9af7y9GDDskS7ZWyZanIwO8IyL7ReQXVlt/Y8xJAOtnP9vSdRYtWzJuZyf0z5jcEnPaxktE3gMGXOGh1caYN6M87SQwxBhzTkSmAltFZBzhP0U7StgpUHFmjZbpSjv6hJ+udbXMwPPAU9brPgU8A/yMm7wd45BseTqaYYypE5F+wLsiUmF3oDgl23Z2Sv+MSUoXfGPMrDie0wg0Wsv7RcQLFBDe098esertQF0iclqvFXNWwplyI+5HZorWnjDXm1lENgDbrLtXy2yHZMvTjjGmzvp5RkTeIDy8cFpEBhpjTlrDImdsDdletGxJtZ2NMadbl5O8f8ZEh3RiJCJ9RcRtLecDI4Ea68/Tr0Wk0Do75zEg2pF3V3kLWCgimSIyjHDWvcAnwEgRGSYiGcBCa90u02F89sdA6xkS0TLbxfZtFY2IZIlIz9Zl4H7C2/EtYLG12mLs74eRomV7C3jMOlunEAi2Dv3YwUH9MzZ2f2qcrDfCv+QThI/mTwNvW+0PAkcJf1J/APhRxHOmEe4YXuDPWF9ssyur9dhqK08lEWcNET4rosp6bLUN2/dl4DBwiPCbaOC1MtvYF2zdVlfJlW/1w3KrT6622nsD7wOfWT9vsynfvwgPgTZZ/fPn0bIRHir5i7WNDxNxtplNOR3TP2O56TdtlVIqReiQjlJKpQgt+EoplSK04CulVIrQgq+UUilCC75SSqUILfhKKZUitOArpVSK0IKvlFIp4v/lL813OxHKhwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "dso = latlon_to_scrip(nx=6, ny=9, lon0=-180.)\n",
    "dims_grid = tuple(dso.grid_dims.values[::-1])\n",
    "ny, nx = dims_grid\n",
    "\n",
    "LON = dso.grid_center_lon.values.reshape(dims_grid)\n",
    "LAT = dso.grid_center_lat.values.reshape(dims_grid)\n",
    "\n",
    "x_corner = dso.grid_corner_lon.values.reshape(dims_grid+(4,))\n",
    "y_corner = dso.grid_corner_lat.values.reshape(dims_grid+(4,))\n",
    "\n",
    "ind = np.array([0, 1, 2, 3, 0])\n",
    "for i, j in product(range(nx), range(ny)):\n",
    "    p = plt.plot(LON[j, i], LAT[j, i], '.')\n",
    "    plt.plot(x_corner[j, i, ind], y_corner[j, i, ind], 'x-', color=p[0].get_color())\n",
    "    \n",
    "dso.info()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate grid file for ETOPO1 dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "writing /glade/work/mclong/esmlab-regrid/latlon_0.0166x0.0166_180W.nc\n",
      "xarray.Dataset {\n",
      "dimensions:\n",
      "\tgrid_corners = 4 ;\n",
      "\tgrid_rank = 2 ;\n",
      "\tgrid_size = 233280000 ;\n",
      "\n",
      "variables:\n",
      "\tint32 grid_dims(grid_rank) ;\n",
      "\tfloat64 grid_center_lat(grid_size) ;\n",
      "\t\tgrid_center_lat:units = degrees ;\n",
      "\tfloat64 grid_center_lon(grid_size) ;\n",
      "\t\tgrid_center_lon:units = degrees ;\n",
      "\tfloat64 grid_corner_lat(grid_size, grid_corners) ;\n",
      "\t\tgrid_corner_lat:units = degrees ;\n",
      "\tfloat64 grid_corner_lon(grid_size, grid_corners) ;\n",
      "\t\tgrid_corner_lon:units = degrees ;\n",
      "\tint32 grid_imask(grid_size) ;\n",
      "\t\tgrid_imask:units = unitless ;\n",
      "\tfloat64 grid_area(grid_size) ;\n",
      "\t\tgrid_area:units = radians^2 ;\n",
      "\t\tgrid_area:long_name = area weights ;\n",
      "\n",
      "// global attributes:\n",
      "\t:title = 0.016666666666666666 x 0.016666666666666666 (lat x lon) grid ;\n",
      "\t:created_by = latlon_to_scrip ;\n",
      "\t:date_created = 2019-04-07 11:48:46.217945 ;\n",
      "\t:conventions = SCRIP ;\n",
      "}CPU times: user 22.2 s, sys: 24.1 s, total: 46.3 s\n",
      "Wall time: 46.7 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "file_out = '/glade/work/mclong/esmlab-regrid/latlon_0.0166x0.0166_180W.nc'\n",
    "dso = latlon_to_scrip(nx=21600, ny=10800, lon0=-180., file_out=file_out)\n",
    "\n",
    "# extract grid coordinates to make check below\n",
    "dims_grid = tuple(dso.grid_dims.values[::-1])\n",
    "ny, nx = dims_grid\n",
    "lat = dso.grid_center_lat.values.reshape(dims_grid)[:, 0]\n",
    "lon = dso.grid_center_lon.values.reshape(dims_grid)[0, :]\n",
    "\n",
    "# check that this conforms to the grid of the original dataset\n",
    "file_src_data = '/glade/work/mclong/etopo1/ETOPO1_Bed_c_gmt4.nc'\n",
    "ds = xr.open_dataset(file_src_data)\n",
    "np.testing.assert_allclose(lat, ds.y.values)\n",
    "np.testing.assert_allclose(lon, ds.x.values)\n",
    "\n",
    "dso.info()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:dev]",
   "language": "python",
   "name": "conda-env-dev-py"
  },
  "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}