jgomezdans/xarray-temporal-fit.ipynb

## xarray-temporal-fit.ipynb
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import dask.array as da\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import xarray as xr\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(x,y):\n",
    "    ''' simple least squares fit'''\n",
    "    return np.polyfit(x, y, 1)[0]\n",
    "\n",
    "def g(x,y):\n",
    "    ''' simple least squares fit'''\n",
    "    return np.polyfit(x, y, 1)[1]\n",
    "\n",
    "def get_slope(x, y, dim='time'):\n",
    "    # x = Pixel value, y = a vector containing the date, dim == dimension\n",
    "    DS1 = xr.apply_ufunc(\n",
    "        f, x , y,\n",
    "        input_core_dims=[[dim], [dim]],\n",
    "        vectorize=True, \n",
    "        dask='parallelized',\n",
    "        output_dtypes=[float],\n",
    "        )\n",
    "    return DS1.rename({'field':'p0'})\n",
    "\n",
    "def get_intercept(x, y, dim='time'):\n",
    "    DS2 = xr.apply_ufunc(\n",
    "        g, x , y,\n",
    "        input_core_dims=[[dim], [dim]],\n",
    "        vectorize=True, \n",
    "        dask='parallelized',\n",
    "        output_dtypes=[float],\n",
    "        )\n",
    "    return DS2.rename({'field':'p1'})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<xarray.Dataset>\n",
      "Dimensions:  (lat: 3, lon: 4, time: 5)\n",
      "Coordinates:\n",
      "  * time     (time) int64 0 1 2 3 4\n",
      "Dimensions without coordinates: lat, lon\n",
      "Data variables:\n",
      "    field    (time, lat, lon) float64 dask.array<shape=(5, 3, 4), chunksize=(1, 3, 4)>\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([[[0.14839575, 0.4925233 , 0.44542527, 0.19103748],\n",
       "        [0.61474645, 0.61937579, 0.47499865, 0.8210922 ],\n",
       "        [0.33607644, 0.34127119, 0.74631819, 0.99013419]],\n",
       "\n",
       "       [[1.66166221, 1.18530321, 1.63593271, 1.82127854],\n",
       "        [1.32384755, 1.14598867, 1.50385216, 1.34005422],\n",
       "        [1.18531918, 1.51841097, 1.84682563, 1.1976601 ]],\n",
       "\n",
       "       [[2.14621515, 2.5799377 , 2.8170234 , 2.81498151],\n",
       "        [2.31027184, 2.70145675, 2.19297758, 2.87434952],\n",
       "        [2.59357375, 2.03970408, 2.24268609, 2.24780081]],\n",
       "\n",
       "       [[3.24278498, 3.0836049 , 3.51158625, 3.23522646],\n",
       "        [3.71424349, 3.93648939, 3.24764248, 3.23555226],\n",
       "        [3.58591781, 3.77918273, 3.23180579, 3.75257605]],\n",
       "\n",
       "       [[4.2034127 , 4.4668951 , 4.1021313 , 4.0664199 ],\n",
       "        [4.82141666, 4.89286553, 4.77839903, 4.29896122],\n",
       "        [4.69204694, 4.11801724, 4.80787637, 4.83397485]]])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Give each pixel a different set of noise\n",
    "nt, ny, nx = 5, 4, 3\n",
    "images = np.ones([ny,nx,nt])\n",
    "values = np.arange(nt)\n",
    "noise = np.random.random([ny,nx,nt])\n",
    "# NOTE: more slope variability:\n",
    "#cube = (images*values*noise).T\n",
    "# alternatively:\n",
    "cube = ((images*values) + noise).T\n",
    "\n",
    "#time = pd.date_range(start='1950-01-01', periods=nt, freq='10D')\n",
    "time = np.arange(nt)\n",
    "DA = xr.DataArray(da.from_array(cube, chunks=(1, nx, ny)),\n",
    "                   dims=('time', 'lat', 'lon'),\n",
    "                   coords={'time': time})\n",
    "DS = DA.to_dataset(name='field')\n",
    "print(DS)\n",
    "DS['field'].values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds = DS.chunk(chunks={'time':-1, 'lat':3, 'lon':2})\n",
    "ds1 = get_slope(ds.time, ds, 'time').compute()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get intercept as well\n",
    "ds2 = get_intercept(ds.time, ds, 'time').compute()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Can merge these into a single datasets\n",
    "R = xr.merge([ds1,ds2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the results to make sure things look okay\n",
    "# So to plot points and line:\n",
    "def plot_pixel(x=0,y=0):\n",
    "    data = ds['field'].isel(dict(lon=x,lat=y)).compute()\n",
    "    #print(data)\n",
    "    coef = R.isel(dict(lon=x,lat=y)).to_array().values\n",
    "    #print(coef)\n",
    "    plt.title(f'Lon={x} Lat={y}')\n",
    "    plt.scatter(data.time.values, data.values, label='data')\n",
    "    plt.plot(data.time.values, np.polyval(coef, data.time.values), label='fit')\n",
    "    plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAEICAYAAAB/Dx7IAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAAIABJREFUeJzt3XlcVXX+x/HXF0UBBVFwA1RcEZVcwiXLshUrc6tpmvaacppmn8kZMdssy8aZpmaqaWyfaZ1JJJfUNsv2cklBEfcFcEEUAVlk+f7+AP0ZoV6Se89d3s/Hg0cX7vfe+/bEfXM498O5xlqLiIj4jiCnA4iISOOouEVEfIyKW0TEx6i4RUR8jIpbRMTHqLhFRHyMiltExMeouKXJGWO2G2MucvDxrzXG7DDGHDbGpBtj2p1krTXG9PoBj/GD/43GmDnGmGxjTI0x5uYfch8S2FTc4leMMf2BfwE3AB2BUuBpR0N93xrgTmCV00HEN6m4xaOMMbcbYzYbYw4YY+YbY2KOu84aY+4wxmwyxhw0xjxljDGNfIjrgAXW2uXW2hLgHmCSMSa8kTl7GmM+NMYUGGP2G2NeNcZE1l33H6ArsMAYU2KM+WNj7tta+5S19gOgvDG3EzlKxS0eY4y5AHgEuBroDOwA3qi3bCwwFBhYty6l7rbnGGMKT/JxTt3t+1O7RwuAtXYLcATo09i4dVljgESgC3B/3X3eAOwErrDWtrbW/rku48nyTW3k44ucUHOnA0hAuQ54wVq7CsAYkwocNMbEW2u3162ZZa0tBAqNMcuAQcASa+2nQKQLj9EaOFTva4eARu1xW2s3A5vrPs03xjwG3HeK27iST+S0aY9bPCmG2r1sAOoOZRQAscet2XPc5VJqi7gxSoCIel+LAIobcyfGmA7GmDeMMbnGmCLgFSC6kVlE3ELFLZ6UB3Q7+okxphUQBeSe6obGmFF1x5NP9DGqbuk6ag+zHL1dD6AlsLGRWR8BLHCGtTYCuJ7awydHfe+0mqfIN62Rjy9yQjpUIu4SbIwJOe7zKuA14A1jzGtAFvAw8NVxh0lOyFr7Ca7tfb8KfFFX5KuAGUCatfZke9wt6mWtpPbQyiFqD9nEAlPq3WYv0KNeRpd+OzDGtKB2p8nw/9vpiLW2xpXbi2iPW9zlHaDsuI/76yYp7gHmAruBnsA1Tfmg1tp1wB3UFvg+agv4zlPcbF29rLcADwBDqC3vRUBavds8Akyve+HxrkbGfLfucUYCc+oun9vI+5AAZvRGCiIivkV73CIiPkbFLSLiY1TcIiI+RsUtIuJjXBoHNMZsp/YPGKqBKmtt8snWR0dH2/j4+NMOJyISKFauXLnfWtvelbWNmeM+31q735WF8fHxrFixohF3LSIS2IwxO069qpYOlYiI+BhXi9sC7xpjVhpjJje0wBgz2RizwhizIj8/v+kSiojId7ha3Gdba4cAlwK/MMZ876+8rLVzrLXJ1trk9u1dOkwjIiI/gEvHuK21eXX/3WeMmQcMA5Y35oEqKyvJycmhvNx/zx0fEhJCXFwcwcHBTkcRET92yuKuO4NbkLW2uO7yJdSeuKdRcnJyCA8PJz4+nsa/qYn3s9ZSUFBATk4O3bt3dzqOiPgxVw6VdAQ+NcasAb4GFllrlzT2gcrLy4mKivLL0gYwxhAVFeXXv1GIiHc45R63tXYrx53f+HT4a2kf5e//PhH5vvTVucxemk1eYRkxkaFMSUlgwuDYU9/wNOh83CIiP1D66lxS0zIoq6wGILewjNS0DAC3lndAzXH//e9/JzExkbZt2zJr1iwA0tPTWb9+vcPJRMQXzV6afay0jyqrrGb20my3Pm5A7XE//fTTLF68+DsvHqanpzN27Fj69evnYDIR8UV5hWWN+npTCZg97jvuuIOtW7cybtw4/va3v/HLX/6Szz//nPnz5zNlyhQGDRrEli1bnI4pIj4kJjK0UV9vKo7scT+wYB3r84qa9D77xURw3xX9T3j9M888w5IlS1i2bBkLFy4EYOTIkYwbN46xY8dy1VVXNWkeEfF/U1ISvnOMGyA0uBlTUhLc+rgBdahERKQpHX0B8s9LNpB3qJyYNiH8cUxf/5wqOdmesYiIL+nbOZwOESHkHSpn5qQkzk/o4PbHDPg97vDwcIqLi52OISI+pryymn98uIl/fbyVNqHBPHHNIEb38cx5mgLmxckTueaaa5g9ezaDBw/Wi5Mi4pLPt+zn0ic+4allW5gwOJb3f38e4wfFeuyP8Iy1tsnvNDk52dZ/I4WsrCwSExOb/LG8TaD8O0UCUWHpEWYuyuJ/K3PoFhXGwxOTOLtXdJPctzFm5aneXeyogD9UIiJyKtZaFqzdzYwF6zhYWsnPR/fkNxf2JiS4mSN5VNwiIieRc7CU6emZfJSdz8C4Nvz71uH0i4lwNJOKW0SkAdU1lhc/28Zf392IMXDv2H7cNDKeZkHOn0xOxS0iUs+6vEOkpmWwNucQF/TtwIMTBhDr5r+GbAwVt4hInbIj1Tz+wUae+2QbbcNa8OS1g7k8qbPXnbJZxS0iAny6aT/T5mWw80Ap1wztQuqlibQJ8863IQzY4r7//vtp3bo1d911V4PXp6en06dPH501UMTPHTh8hIcWrSdtVS49olvx+u0jOKtnlNOxTipgi/tUdLpXEf9mrSX921weXJhFUVklv7qgF784v5djI36N4bV/OZm+OpezZ31I96mLOHvWh6Svzj3t+5w5cyYJCQlcdNFFZGfXnuj82WefZejQoQwcOJArr7yS0tLSBk/32tA6EfFNOwtKufGFr/ndm2voFhXGol+P4g+XJPhEaYOXFvfRtwPKLSzD8v9vB3Q65b1y5UreeOMNVq9eTVpaGt988w0AkyZN4ptvvmHNmjUkJiby/PPPHzvd6+zZs/n222/p2bNng+tExLdUVdcwZ/kWLnn8Y1bvLGTG+P68dcdIEjqFOx2tUbzyUMnJ3g7oh54u8ZNPPmHixImEhYUBMG7cOAAyMzOZPn06hYWFlJSUkJKS0uDtXV0nIt4pI+cQU9PWsi6viIsSO/LghP50buM9I36N4ZXF7a63A2popOfmm28mPT2dgQMH8tJLL/HRRx81eFtX14mIdyk9UsVj727khc+2EdW6Jf+8bghjBnTyuhG/xvDKQyXueDugc889l3nz5lFWVkZxcTELFiwAoLi4mM6dO1NZWcmrr756bH39072eaJ2IeK+PN+Zzyd+W89yn27hmWFfe//15XOqFc9mN5ZXFPSUlgdB6LxKc7tsBDRkyhB//+McMGjSIK6+8klGjRgHw4IMPMnz4cC6++GL69u17bH39072eaJ2IeJ/9JRX85o3V3PTC17RsHsR/f3YWD09Mok2od85lN5bXntY1fXUus5dmk1dYRkxkKFNSEtz+dkBNQad1FXGOtZa3VuYw850sDldUcefoXtx5fk9aNvf+aRG/OK3rhMGxPlHUIuIdtu8/zLR5GXy+pYDkbm15ZFISvTv61rSIq7y2uEVEXFFZXcOzn2zlifc30aJZEA9NGMC1w7oS5AVn8XMXjxa3tdbnXxQ4GXccdhKRE/t2VyFT565lw55ixvTvxAPj+9MxIsTpWG7nseIOCQmhoKCAqKgovyxvay0FBQWEhPj/N42I00oqqvjru9m89Pl2OoaH8K8bziSlfyenY3mMx4o7Li6OnJwc8vPzPfWQHhcSEkJcXJzTMUT82ocb9jJ9Xia7i8q5YUQ3pqQkEB7iH9MirvJYcQcHB9O9e3dPPZyI+Jl9xeU8sGA9i9bupk/H1rx17Vmc2a2d07EcoRcnRcSrWWv574pdzFyURXllDX+4uA8/O68nLZp75Z+heISKW0S81pb8EqalZfDVtgMM696ORyYl0bN9a6djOU7FLSJe50hVDf/6eAv/WLaZkOZBzJqUxNXJXfx6xK8xVNwi4lVW7jhIatpaNu4t4fIzOnPfFf3oEK5preO5XNzGmGbACiDXWjvWfZFEJBAVl1cye2k2//lyB50iQnjuxmQu6tfR6VheqTF73L8BsoAIN2URkQD17ro93Pv2OvYWl3PTWfHclZJA65Y6IHAiLm0ZY0wccDkwE/i9WxOJSMDYW1TO/fPXsThzD307hfPP64cwuGtbp2N5PVd/pD0O/BE44RlbjDGTgckAXbt2Pf1kIuK3amosr3+zk1mLN1BRVcOUlAQmn9uD4GaBO+LXGKcsbmPMWGCftXalMWb0idZZa+cAc6D2tK5NllBE/MrmfcWkpmXwzfaDjOwZxcyJSXSPbuV0LJ/iyh732cA4Y8xlQAgQYYx5xVp7vXujiYg/qaiq5ullW3j6o82EtWjO7KvO4Koz4/zy3EXudsrittamAqkAdXvcd6m0RaQxvtl+gKlz17Il/zDjBsZw7xX9iG7d0ulYPksv24qI2xwqq+TRJRt47audxEaG8uItQzk/oYPTsXxeo4rbWvsR8JFbkoiI37DWsiRzD/fNX8f+kgpuO6c7v7u4D6004tcktBVFpEntPlTGvW+v4731e+nXOYLnbkrmjLhIp2P5FRW3iDSJmhrLK1/t4M9LsqmqqSH10r7cek53jfi5gYpbRE5b9p5iUtPWsmpnIef0imbmxAF0i9KIn7uouEXkByuvrOapZZt55uMttG7ZnMeuHsjEwbEa8XMzFbeIfE/66lxmL80mr7CMmMhQpqQkMGFw7HfWfLm1gGlpGWzdf5hJg2O5+/JEojTi5xEqbhH5jvTVuaSmZVBWWQ1AbmEZqWkZAEwYHMuh0koefieLN1fsoku7UP596zDO7dPeycgBR8UtIt8xe2n2sdI+qqyymj8v2UCzIMMDC9ZzsPQIPzu3B7+5qDdhLVQjnqYtLiLfkVdY1vDXD5Xzq9dXkxTbhpduGcqA2DYeTiZHqbhF5DtiIkPJbaC8DXD35YncPDKe5hrxc5S2voh8x5SUBEKDm33na0EGpo9N5LZRPVTaXkD/B0TkO8YM6MTIXlHHPm8bFsxffzSQn57Tw8FUcjwdKhGRYz7bvJ9p8zLYUVDKVWfGcfdlibRt1cLpWFKPiltEOHj4CA8tymLuqhy6RYXx2m3DGdkr2ulYcgIqbpEAZq1l/po8ZixYz6GySu4c3ZNfX9ibkHrHuMW7qLhFAtSuA6VMT8/k4435DIxrwyu3DSexc4TTscQFKm6RAFNVXcOLn23nsfc2Ygzcd0U/bjwrnmZBOr+Ir1BxiwSQzNxDpKZlkJF7iAv6duDBCQOIjQx1OpY0kopbJACUHanmb+9v5PlPt9E2rAVPXjuYy5M66yx+PkrFLeLnlm/M5+70DHYdKOOaoV1IvTSRNmHBTseS06DiFvFTBSUVPLQoi3mrc+kR3YrXbx/BWT2jTn1D8XoqbhE/Y60lbVUuDy1aT3F5Fb+6oBe/OL+XRvz8iIpbxI/sLCjl7vQMPtm0nyFdI3lk0hkkdAp3OpY0MRW3iB+oqq7huU+38fj7G2keFMSM8f25fng3gjTi55dU3CI+bm1OIVPnZrB+dxEX9+vIjPH96dxGI37+TMUt4qMOV1Tx2HsbefGzbUS3bskz1w8hpX8njfgFABW3iA9alr2P6fMyyS0s49rhXfnTmL60CdWIX6BQcYv4kPziCmYsXM+CNXn0bN+K/91xFkPj2zkdSzxMxS3iA6y1/G9lDjMXZVF6pIrfXtSbn4/uScvmGvELRCpuES+3bf9hpqVl8MXWApK7tWXWlUn06qARv0Cm4hbxUpXVNcxZvpUnPthEy2ZBzJw4gJ8M7aoRP1Fxi3ij1TsPkpqWwYY9xYzp34kHxvenY0SI07HES6i4RbxISUUVf1mazctfbKdjeAj/uuFMUvp3cjqWeBkVt4iX+CBrL/ekZ7K7qJwbRnRjSkoC4SEa8ZPvU3GLOGxfcTkPzF/Poozd9OnYmreuHcmZ3do6HUu8mIpbxCE1NZY3V+zikXeyKK+s4Q8X9+Fn5/WkRfMgp6OJl1Nxizhg874Sps3L4OttBxjevR0PT0qiZ/vWTscSH3HK4jbGhADLgZZ169+y1t7n7mAi/uhIVQ3PfLyFJz/cTEhwELMmJXF1cheN+EmjuLLHXQFcYK0tMcYEA58aYxZba790czYRv7JyxwGmzs1g074SLj+jM/dd0Y8O4Rrxk8Y7ZXFbay1QUvdpcN2HdWcoEX9SVF7J7CXZvPLVDjpHhPD8TclcmNjR6Vjiw1w6xm2MaQasBHoBT1lrv2pgzWRgMkDXrl2bMqOIz1q6bg/3vp3JvuIKbjornrtSEmjdUi8tyelx6TvIWlsNDDLGRALzjDEDrLWZ9dbMAeYAJCcna49cAtreonLufTuTpev20rdTOP+6IZlBXSKdjiV+olE/+q21hcaYj4AxQOYplosEnJoay2tf7+TRxRuoqK5hSkoCk8/tQXAzjfhJ03FlqqQ9UFlX2qHARcCjbk8m4mM27S0mNS2DFTsOMrJnFDMnJtE9upXTscQPubLH3Rl4ue44dxDwX2vtQvfGEvEdFVXVPL1sC09/tJmwFs2ZfdUZXHVmnN5CTNzGlamStcBgD2QR8TlfbztAatpatuQfZvygGO4Z24/o1i2djiV+Ti9vi/wAh8oqmbV4A69/vZPYyFBevGUo5yd0cDqWBAgVt0gjWGtZkrmH++avY39JBbed053fX9KHsBZ6Konn6LtNxEW7D5VxT/o63s/aS/+YCJ6/aShJcW2cjiUBSMUtcgrVNZZXvtzB7KXZVNXUkHppX356Tneaa8RPHKLiFjmJDXuKSE3LYPXOQkb1jmbmhCS6RoU5HUsCnIpbpAHlldU8+eFmnvl4C+EhzXns6oFMHByrET/xCipukXq+2FLAtHkZbNt/mEmDY5k+th/tWrVwOpbIMSpukTqFpUd45J0NvLliF13ahfKfnw5jVO/2TscS+R4VtwQ8ay0L1+7mgQXrOFhayc/O7cFvL+pDaItmTkcTaZCKWwJabmEZ96Rn8uGGfSTFtuGlW4YxIFYjfuLdVNwSkKprLC9/vp2/vJuNtTD98kRuHhmvET/xCSpuCTjr84pITVvLmpxDnNenPQ9NGECXdhrxE9+h4paAUV5ZzRMfbGLO8q1EhgbzxDWDGDcwRiN+4nNU3BIQPtu8n2nzMthRUMqPzozj7ssTiQzTiJ/4JhW3+LWDh4/w0KIs5q7KIT4qjNduG87IXtFOxxI5LSpu8UvWWt7+No8ZC9dTVFbJnaN78usLexMSrBE/8X0qbvE7uw6Ucnd6Jss35jOwSySzJiWR2DnC6VgiTUbFLX6jqrqGFz/bzmPvbcQYuO+Kftx4VjzNgvTio/gXFbf4hczcQ0xNW0tmbhEX9u3AjAkDiI0MdTqWiFuouMWnlR6p4vH3N/H8p9toG9aCJ68dzOVJnTXiJ35NxS0+a/nGfO5Oz2DXgTKuGdqF1EsTaRMW7HQsEbdTcYvPKSip4KFFWcxbnUuP6Fa8MXkEI3pEOR1LxGNU3OIzrLWkrcrloUXrKamo4tcX9OLO83u5NOKXvjqX2UuzySssIyYylCkpCUwYHOuB1CJNT8UtPmFHwWHunpfJp5v3M6RrJLOuPIM+HcNdum366lxS0zIoq6wGas8ImJqWAaDyFp+k4havVlldw/OfbuPx9zfSPCiIB8f357rh3QhqxIjf7KXZx0r7qLLKamYvzVZxi09ScYvXWptTyJ/mZpC1u4iL+3Vkxvj+dG7T+BG/vMKyRn1dxNupuMXrHK6o4rH3NvLiZ9uIbt2SZ64fwpgBnX/w/cVEhpLbQEnHaM5bfJSKW7zKsux9TJ+XSW5hGdcN78qfLu1LRMjpjfhNSUn4zjFugNDgZkxJSTjduCKOUHGLV8gvrmDGwvUsWJNHrw6t+d8dZzE0vl2T3PfR49iaKhF/oeIWR1lr+d/KHGYuyqLsSDW/vag3Px/dk5bNm/YsfhMGx6qoxW+ouMUx2/YfZlpaBl9sLWBofFsemZRErw6ujfiJBDIVt3jckaoanv1kK098sImWzYKYOXEAPxnatVEjfiKBTMUtHrV650Gmzs0ge28xlyV14r4r+tMxIsTpWCI+RcUtHlFSUcVflmbz8hfb6RgewrM3JnNxv45OxxLxSSpucbv31+/lnrcz2VNUzo0junFXSgLhpzniJxLITlncxpguwL+BTkANMMda+4S7g4nv21dUzgML1rMoYzd9OrbmyWtHcma3tk7HEvF5ruxxVwF/sNauMsaEAyuNMe9Za9e7OZv4qJoay5srdvHwO1lUVNVw1yV9mHxuT1o0D3I6mohfOGVxW2t3A7vrLhcbY7KAWEDFLd+zeV8J0+Zl8PW2Awzv3o5HJiXRo31rp2OJ+JVGHeM2xsQDg4GvGrhuMjAZoGvXrk0QTXzJkaoanvl4C09+uJmQ4CAevTKJq5O76C3ERNzA5eI2xrQG5gK/tdYW1b/eWjsHmAOQnJxsmyyheL2VOw4wdW4Gm/aVMPaMztx7RT86hGvET8RdXCpuY0wwtaX9qrU2zb2RxFcUlVcye0k2r3y1g84RIbxwczIX9NWIn4i7uTJVYoDngSxr7WPujyS+YOm6Pdz7dib5xRXcMrI7f7ikD61aarpUxBNceaadDdwAZBhjvq372jRr7TvuiyXeam9ROfe+ncnSdXvp2ymcOTckM7BLpNOxRAKKK1MlnwJ6hSnA1dRYXvt6J48u3sCR6hr+NKYvt43qTnAzjfiJeJp+t5VT2rS3mNS0DFbsOMjZvaKYOSGJ+OhWTscSCVgqbjmhiqpqnlq2hX9+tJlWLZvzlx8N5MohsRrxE3GYilsa9PW2A0xNW8vW/MNMGBTD9LH9iG7d0ulYIoKKW+o5VFbJrMUbeP3rncS1DeWlW4YyOqGD07FE5DgqbgFq30JsceYe7pu/joKSCm4f1Z3fXdyHsBb6FhHxNnpWCnmFZdz79jrez9pL/5gIXrhpKElxbZyOJSInoOIOYNU1lle+3MGfl2yg2lqmXdaXW8/uTnON+Il4NRV3gNqwp4jUtAxW7yxkVO9oHp6YRJd2YU7HEhEXqLgDTHllNU9+uJlnPt5CRGgwj/94EOMHxWjET8SHqLgDyBdbCpg2L4Nt+w8zaUgs0y/vR7tWLZyOJSKNpOIOAIWlR3jknQ28uWIXXduF8cpPh3NO72inY4nID6Ti9mPWWhau3c0DC9ZxsLSSO87ryW8u7E1oi2ZORxOR06Di9lO5hWXck57Jhxv2cUZcG16+dRj9YzTiJ+IPVNx+prrG8vLn2/nLu9kA3DO2HzePjKdZkF58FPEXKm4/sj6viNS0tazJOcTohPY8NGEAcW014ifib1TcfqC8sprH39/Es59spW1YMH//yWCuOKOzRvxE/JSK28d9umk/d6dnsKOglKuT45h2WSKRYRrxE/FnKm4fdfDwER5alMXcVTnER4Xx2u3DGdlTI34igUDF7WOstbz9bR4zFq6nqKySX5zfk19d0JuQYI34iQQKFbcP2XWglLvTM1m+MZ9BXSKZdWUSfTtFOB1LRDxMxe0DqqprePGz7Tz23kaCDDwwrj/Xj+imET+RAKXi9nKZuYeYmraWzNwiLkrswIzxA4iJDHU6log4SMXtpUqPVPH4+5t4/tNttGvVgqevG8KlAzodG/FLX53L7KXZ5BWWERMZypSUBCYMjnU4tYh4gorbCy3fmM/d6RnsOlDGT4Z1YeqYRNqEBR+7Pn11LqlpGZRVVgO1f96empYBoPIWCQAqbi9SUFLBQ4uymLc6lx7tW/Hm5BEM7xH1vXWzl2YfK+2jyiqrmb00W8UtEgBU3F7AWsvcVbk8tGg9hyuq+PWFvblzdM8TjvjlFZY16usi4l9U3A7bUXCYafMy+GxzAWd2a8sjk5Lo0zH8pLeJiQwlt4GS1ouWIoFBxe2QyuoanvtkG4+/v5HgZkE8OGEA1w3rSpALI35TUhK+c4wbIDS4GVNSEtwZWUS8hIrbAWt2FTI1LYOs3UVc0q8jM8YPoFObEJdvf/Q4tqZKRAKTituDDldU8dd3N/LS59toH96SZ64/kzEDOv2g+5owOFZFLRKgVNwesmzDPqanZ5JbWMb1I7ryxzF9iQgJPvUNRUTqUXG7WX5xBTMWrmfBmjx6d2jNW3ecRXJ8O6djiYgPU3G7ibWW/63IYeY7WZQdqeZ3F/XhjtE9aNlcZ/ETkdOj4naDrfklTJuXwZdbDzAsvh0PT0qiV4fWTscSET+h4m5CR6pqmLN8C3//cDMtmwfxyKQkfpzcxaURPxERV6m4m8iqnQdJnZtB9t5iLkvqxP1X9KdDhOsjfiIirjplcRtjXgDGAvustQPcH8m3lFRUMXvJBv795Q46RYTw7I3JXNyvo9OxRMSPubLH/RLwJPBv90bxPe+t38u9b2eyp6icG0d0466UBMI14icibnbK4rbWLjfGxLs/iu/YV1TO/QvW8U7GHhI6hvPUdUMY0rWt07FEJEA02TFuY8xkYDJA165dm+puvUpNjeWNb3bxyOIsKqpqmJKSwO2jetCieZDT0UQkgDRZcVtr5wBzAJKTk21T3a+32LyvhGlpGXy9/QAjerTj4YlJ9GivET8R8TxNlZzCkaoa/vnRFp5atpnQFs3485Vn8KPkuGNvISYi4mkq7pNYsf0AqWkZbNpXwhUDY7h3bD/ah7d0OpaIBDhXxgFfB0YD0caYHOA+a+3z7g7mpKLySh5dvIFXv9pJbGQoL948lPP7dnA6logI4NpUyU88EcRbLMncw33zM8kvruDWs7vzh0v60KqlfjEREe+hRqqz51A5976dybvr95LYOYI5NyQzsEuk07FERL4n4Iu7psby6lc7eHRJNpXVNfxpTF9uG9Wd4GYa8RMR7xTQxb1xbzGpaRms3HGQs3tF8fDEJLpFtXI6lojISQVkcZdXVvP0ss388+MttGrZnL/+aCCThsRqxE9EfELAFfdXWwtInZfB1vzDTBgUwz1j+xHVWiN+IuI7Aqa4D5VVMmtxFq9/vYu4tqG8fOswzuvT3ulYIiKN5vfFba3lnYw93L9gHQUlFdw+qju/u7gPYS38/p8uIn7Kr9srr7CMe9Iz+WDDPvrHRPDizUMZENvG6VgiIqfFa4o7fXUus5dmk1dr0ulzAAAGxklEQVRYRkxkKFNSEpgwOPYH3Vd1jeU/X2xn9tJsqq3l7ssSueXseJprxE9E/IBXFHf66lxS0zIoq6wGILewjNS0DIBGl/eGPUVMnZvBt7sKGdU7mocnJtGlXViTZxYRcYpXFPfspdnHSvuosspqZi/Ndrm4yyur+fsHm5izfCsRocE8/uNBjB8UoxE/EfE7XlHceYVljfp6fZ9v2c+0tAy2F5Ry5ZA4pl+eSNtWLZoyooiI1/CK4o6JDCW3gZKOiQw96e0KS48wc1EW/1uZQ9d2Ybzy0+Gc0zvaXTFFRLyCVxT3lJSE7xzjBggNbsaUlIQG11trWbB2NzMWrONgaSV3nNeT31zYm9AWzTwVWUTEMV5R3EePY7syVZJzsJTp6Zl8lJ3PwLg2/PvW4fSLifB0ZBERx3hFcUNteZ/shcjqGsuLn23jr+9uxBi4Z2w/bh4ZT7MgvfgoIoHFa4r7ZNblHSI1LYO1OYc4P6E9D04YQFxbjfiJSGDy6uIuO1LN4x9s5LlPttE2LJh//GQwY8/orBE/EQloXlvcn2zK5+55mew8UMrVyXFMuyyRyDCN+ImIeF1xHzh8hIcWridtdS7do1vx2u3DGdlTI34iIkd5VXF/saWAX7y2iqKySn5xfk9+dUFvQoI14icicjyvKu7u0a3oHxPB3Zcn0reTRvxERBriVcXdqU0I//npcKdjiIh4NZ3nVETEx6i4RUR8jIpbRMTHqLhFRHyMiltExMeouEVEfIyKW0TEx6i4RUR8jLHWNv2dGpMP7PiBN48G9jdhnKaiXI2jXI2jXI3jjblON1M3a217Vxa6pbhPhzFmhbU22ekc9SlX4yhX4yhX43hjLk9m0qESEREfo+IWEfEx3ljcc5wOcALK1TjK1TjK1TjemMtjmbzuGLeIiJycN+5xi4jISai4RUR8jGPFbYwZY4zJNsZsNsZMbeD6lsaYN+uu/8oYE+8luW42xuQbY76t+7jNA5leMMbsM8ZknuB6Y4z5e13mtcaYIe7O5GKu0caYQ8dtq3s9lKuLMWaZMSbLGLPOGPObBtZ4fJu5mMvj28wYE2KM+doYs6Yu1wMNrPHo89HFTB5/Lh732M2MMauNMQsbuM7928pa6/EPoBmwBegBtADWAP3qrbkTeKbu8jXAm16S62bgSQ9vr3OBIUDmCa6/DFgMGGAE8JWX5BoNLHTg+6szMKTucjiwsYH/jx7fZi7m8vg2q9sGresuBwNfASPqrfHo89HFTB5/Lh732L8HXmvo/5UntpVTe9zDgM3W2q3W2iPAG8D4emvGAy/XXX4LuNAYY7wgl8dZa5cDB06yZDzwb1vrSyDSGNPZC3I5wlq721q7qu5yMZAFxNZb5vFt5mIuj6vbBiV1nwbXfdSfWvDo89HFTI4wxsQBlwPPnWCJ27eVU8UdC+w67vMcvv8NfGyNtbYKOAREeUEugCvrfr1+yxjTxc2ZXOFqbiecVffr7mJjTH9PP3jdr6mDqd1jO56j2+wkucCBbVb3q/+3wD7gPWvtCbeXp56PLmQCZ56LjwN/BGpOcL3bt5VTxd3QT5/6P01dWdPUXHnMBUC8tfYM4H3+/yerk5zYVq5YRe35FwYC/wDSPfngxpjWwFzgt9baovpXN3ATj2yzU+RyZJtZa6uttYOAOGCYMWZAvSUe314uZPL4c9EYMxbYZ61debJlDXytSbeVU8WdAxz/0zEOyDvRGmNMc6AN7v+1/JS5rLUF1tqKuk+fBc50cyZXuLI9Pc5aW3T0111r7TtAsDEm2hOPbYwJprYcX7XWpjWwxJFtdqpcTm6zuscsBD4CxtS7yonn40kzOfRcPBsYZ4zZTu2h1AuMMa/UW+P2beVUcX8D9DbGdDfGtKD2AP78emvmAzfVXb4K+NDWHe13Mle946DjqD1O6bT5wI11kxIjgEPW2t1OhzLGdDp6bM8YM4za77cCDzyuAZ4Hsqy1j51gmce3mSu5nNhmxpj2xpjIusuhwEXAhnrLPPp8dCWTE89Fa22qtTbOWhtPbT98aK29vt4yt2+r5k15Z66y1lYZY34JLKV2kuMFa+06Y8wMYIW1dj613+D/McZspvan1TVekuvXxphxQFVdrpvdncsY8zq10wbRxpgc4D5qX6zBWvsM8A61UxKbgVLgFndncjHXVcDPjTFVQBlwjQd++ELtXtENQEbdMVKAaUDX47I5sc1cyeXENusMvGyMaUbtD4r/WmsXOvx8dCWTx5+LJ+LpbaU/eRcR8TH6y0kRER+j4hYR8TEqbhERH6PiFhHxMSpuEREfo+IWEfExKm4RER/zfxrUPw6q0hYiAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_pixel(0,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Slope map')"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEWCAYAAABxMXBSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAAHNFJREFUeJzt3Xu4ZFV95vHvS9OACnJrhZZuQBy8tAQBO1zUydMz3oBR2ogGMFEw5unRyKNmJHlQM5BBE3HGmIyiYqv9cJGgBo12YhMGudhRh0uLNLcWaBmUphGElgaUAH3OO3/sfaQo6pyzd9dl16nzfnj2c2rXXlX7V9VU/Wqttddask1ERERV2zQdQEREzCxJHBERUUsSR0RE1JLEERERtSRxRERELUkcERFRSxJHDAVJJ0n6ftNxRMT0kjhiYCS9StIPJW2WtEnSDyT9btNxRUQ92zYdQMwOkp4N/AvwHuDrwHbAfwQeazKuiKgvNY4YlBcC2L7Q9pjtR23/H9s3dCos6RWSri1rJ9dKekXLsSslfVzSNeXxb0vareX44WXN5kFJayUtmSwoSXdK+nNJN0j6taQvS9pD0sWSHpb0XUm7tpT/R0m/KM+7WtJLW46dI+lsSZeWj/2epH26fN8ihk4SRwzKbcCYpHMlHdX6ZdyuTALfAT4N7A58CviOpN1bir0D+GPgecCWsiyS9iof+zFgN+AU4BuSnjNFbMcCr6VIbm8ELgY+DMyj+Iy8r6XsxcD+wHOB64AL2p7rD4GPlo+9vsPxiBkviSMGwvZDwKsAA18EfilppaQ9OhT/L8Dtts+3vcX2hcBPKL7UJ5xv+ybbvwb+O/AHkuYAfwSssr3K9rjtS4E1wNFThPcZ2/favhv4N+Bq2z+2/RjwT8DBLa9jhe2Hy2N/BbxM0s4tz/Ud26vL4x8BjpC0sPo7FTH8kjhiYGyvs32S7QXAARS1hb/vUPR5wM/a7vsZsFfL/l1tx+ZS/MrfB3hr2Uz1oKQHKRLW/ClCu7fl9qMd9ncEkDRH0pmSfirpIeDOssy8TnHZfgTYVL6eiJGRxBGNsP0T4ByKBNJuI0UCaLU3cHfL/sK2Y08A91N8cZ9ve5eW7Vm2z+xB2G8DlgKvAXYG9i3vV6e4JO1I0Vy2sQfnjhgaSRwxEJJeLOmDkhaU+wuBE4CrOhRfBbxQ0tskbSvpOGARxVVZE/5I0iJJzwTOAC6yPQZ8BXijpNeXNYQdJC2ZOG+XdqK4CuwB4JnA33Qoc3R52fF2FH0dV9u+q0O5iBkriSMG5WHgMOBqSb+mSBg3AR9sL2j7AeAN5bEHgL8A3mD7/pZi51PUWH4B7EDZgV1+SS+l6Nz+JUUN5M/pzf/r51E0i90N3ELnpPcPwOkUTVQvp+gsjxgpykJOMdNIuhL4iu0vNR1LK0nnABts/2XTsUT0U2ocERFRS2OJQ9JCSVdIWifpZknv71BmSTnQ6vpyO62JWCMi4kmNNVVJmg/Mt32dpJ2AHwFvsn1LS5klwCm239BIkBER8TSN1Ths32P7uvL2w8A6nnqdfkREDKGhmORQ0r4Uo3Ov7nD4CElrKa6FP8X2zR0evwxYBrDDM/XyfV4wt3/BziA/u++5TYcwNOY83nQEw+OJZ+eCmAmP/+zu+21PNR3NtF7/n57lBzaNVSr7oxseu8T2kd2cbxg0njjKQVLfAD5QTkvR6jpgH9uPSDoa+BbFPEFPYXs5sBzgxQdu7y+u7MUl+zPfuz93ctMhDI2d76z2wZ4N7j4q78WEn//xqe0zFNT2wKYxrrlk70pl58y/fd70pYZfo1dVSZpLkTQusP3N9uO2HyqnbcD2KmCupJF44yNiNBgYr/jfqGisxiFJwJeBdbY/NUmZPYF7bVvSoRSJ7oEBhhkRMSVjnvDsqsU12VT1SuDtwI2Sri/v+zDFvEPYPht4C/AeSVsoJps73hmxGBFDZpRqE1U0ljhsf5+nTg7XqcxZwFmDiSgioj5jxmbZ79nGO8cjIma6cZI4IiKiIgNjSRwREVFHahwREVGZgSfSxxEREVUZp6kqIiJqMIzNrryRxBER0Y1i5PjsksQREdEVMTb1kLSRk8QREdGFonM8iSMiIioqxnEkcURERA3jqXFERERVqXFEREQtRow1u7TRwCVxRER0KU1VERFRmRGPe07TYQxUEkdERBeKAYBpqoqIiBrSOR4REZXZYsypcURERA3jqXFERERVRef47PoqnV2vNiKix2Zj5/jserUREX0wZlXapiNphaT7JN00yfEXS/q/kh6TdErbsSMl3SppvaRTe/TSOkriiIjowsTI8SpbBecAR05xfBPwPuCTrXdKmgN8FjgKWAScIGnRVr2gCpI4IiK6NO5tKm3Tsb2aIjlMdvw+29cCT7QdOhRYb/sO248DXwWWdvGSppQ+joiILhSTHFb+DT5P0pqW/eW2l/cgjL2Au1r2NwCH9eB5O0riiIjoghFPVJ9y5H7bi/sQRqcOlL6thJ7EERHRBZthGAC4AVjYsr8A2NivkzX2aiUtlHSFpHWSbpb0/g5lJOnT5VUCN0g6pIlYIyImJ8Yrbn10LbC/pOdL2g44HljZr5M1WePYAnzQ9nWSdgJ+JOlS27e0lDkK2L/cDgM+Tx/b7SIi6jK9q3FIuhBYQtEXsgE4HZgLYPtsSXsCa4BnA+OSPgAssv2QpJOBS4A5wArbN/ckqA4aSxy27wHuKW8/LGkdRQdPa+JYCpxn28BVknaRNL98bETEUOjVQk62T5jm+C8omqE6HVsFrOpJINMYij4OSfsCBwNXtx3qdKXAXpQJJyKiaUZZyGnQJO0IfAP4gO2H2g93eMjTrhSQtAxYBrDH8xp/SRExixh4YpbNVdXopQCS5lIkjQtsf7NDkUpXCthebnux7cW77N741Q0RMauIsYrbqGjyqioBXwbW2f7UJMVWAu8or646HNic/o2IGCamdyPHZ4om61evBN4O3Cjp+vK+DwN7Q3EFAUVHz9HAeuA3wDsbiDMiYkqjVJuoosmrqr5P5z6M1jIG3juYiCIi6rM1UrWJKmZXj05ERI8VneOVpxwZCUkcERFdyZrjERFRQ9E5nj6OiIiooVcjx2eKJI6IiC5k5HhERNQ2nhpHRERUZcMT40kcERFRUdFUlcQRERE1ZOR4RERUlstxIyKipjRVRURETX1eT3zoJHFERHShuKoqc1VFRERFGQAYERG1pakqIiIqy1VVERFRW66qioiIymyxJYkjIiLqSFNVRERUlj6OiIioLYkjIiIqyziOiIioLeM4IiKiMhu2ZCGniIioI01VERFRWfo4IiKiNidxREREHbOtc7zRHh1JKyTdJ+mmSY4vkbRZ0vXldtqgY4yImIpd9HFU2aZT4TtRkj4tab2kGyQd0nJsrOW7cmUPX+LTNF3jOAc4CzhvijL/ZvsNgwknIqIuMda7q6rOYervxKOA/cvtMODz5V+AR20f1KtAptJojcP2amBTkzFERHTLVqVt+ueZ9jtxKXCeC1cBu0ia36OXUVnTNY4qjpC0FtgInGL75vYCkpYBywB2mv9MvvXgywcc4nC68b99rukQhsYLv3di0yEMjReesLbpEIbGz3vwHAOeq2ov4K6W/Q3lffcAO0haA2wBzrT9rX4FMeyJ4zpgH9uPSDoa+BZFFe0pbC8HlgPsuWg3DzbEiJjVXPRzVDSv/HKfsLz8/qqqU4aaOPvetjdK2g+4XNKNtn9a47krG+rEYfuhlturJH1O0jzb9zcZV0REqxpXVd1ve3EXp9oALGzZX0DRGoPtib93SLoSOBjoS+IY6nHykvaUpPL2oRTxPtBsVBERT3LZOV5l64GVwDvKq6sOBzbbvkfSrpK2B5A0D3glcEsvTthJozUOSRcCSyiqbxuA04G5ALbPBt4CvEfSFuBR4Hi7RqUwImIAevWtVOE7cRVwNLAe+A3wzvKhLwG+IGmc4gf2mbZHM3HYPmGa42dRXJoWETG0ejVyvMJ3ooH3drj/h8Dv9CSICoa6jyMiYtjZmXIkIiJqyiSHERFRy2zreU3iiIjoghHjWcgpIiLqmGUVjiSOiIiupHM8IiJqm2VVjiSOiIgupcYRERGVGRgfT+KIiIiqDKTGERERdWQcR0RE1JPEERER1VVbFnaUJHFERHQrNY6IiKjM4FxVFRER9SRxREREHWmqioiIWpI4IiKisgwAjIiIujIAMCIi6slVVRERUYdS44iIiMpMOscjIqIOpXM8IiJqSo0jIiJqGW86gMHapkohSZ+ocl9ExKwzMY6jyjYiKiUO4LUd7juq25NLWiHpPkk3TXJckj4tab2kGyQd0u05IyJ6Ta62DQNJ20r6r5L+tfxeXSvpYknvljS3ynNM2VQl6T3AnwL7Sbqh5dBOwA+2PvTfOgc4CzhvkuNHAfuX22HA58u/ERHDY0iSQkXnAw8CfwVsKO9bAJwIfAU4bronmK6P4x+Ai4GPA6e23P+w7U01g30a26sl7TtFkaXAebYNXCVpF0nzbd/T7bkjImapQ2y/qO2+DRTfsbdVeYIpm6psb7Z9p+0TbP8MeJQit+4oae+tCrmevYC7WvY3lPc9haRlktZIWvObBx8bQFgREU+aSU1VwK8kvVXSb7//JW0j6TjgV1WeoGrn+Bsl3Q78P+B7wJ0UNZF+69Sb9LS33/Zy24ttL37mLtsPIKyIiJIpphypsg2H44G3APdKuq2sZfwCeHN5bFpVO8c/BhwO3Gb7+cCr6U0fx3Q2AAtb9hcAGwdw3oiI6lxxGwJlK9JxwD7AF4DbgO8DVwGVugGqJo4nbD8AbCNpG9tXAAdtRcx1rQTeUV5ddTiwOf0bETFsZlhT1YRzgJcAfwt8prx9fpUHVh0A+KCkHYHVwAWS7gO21I/zqSRdCCwB5knaAJwOzAWwfTawCjgaWA/8Bnhnt+eMiOi54UsKVbzI9sta9q+QtLbKA6smjqXAvwN/BvwhsDNwRq0QO7B9wjTHDby32/NERPRVjxKHpBXAG4D7bB/Q4biA/03xg/o3wEm2ryuPnQj8ZVn0Y7bPneZ0P5Z0uO2ryscfRsUuiEqJw/avW3anCyYiYtbocTPUOWzF2DZJu1G02CymSGM/krTS9lRXSR1G0RXw83J/b2CdpBspfrcfONkDpxsA+DCdc6nKJ372VI+PiJgVenTF1NaObaNo8r90YnydpEuBI4ELp3iuI7c2zikTh+2dtvaJIyJmixo1jnmS1rTsL7e9vMapJhvbVmnMW6tybN5Wyey4ERHdqp447re9uIszTTa2rdKYt16pejluRER0UvFS3B71g0w2tm2gY96SOCIiujW4AYCTjW27BHidpF0l7Qq8rryvL9JUFRHRJfVoIaetHdtme5OkjwLXlk91Ri8mop1MEkdExJDoZmyb7RXAin7E1S6JIyKiWzNz5PhWS+KIiOjGcM5D1VdJHBER3UriiIiIWpI4IiKiKtG7q6pmiiSOiIhupI8jIiJqS+KIiIhakjgiIqKONFVFREQ9SRwREVGZc1VVRETUlRpHRETUkT6OiIioJ4kjIiIq690iTTNGEkdERBdEmqoiIqKmJI6IiKgniSMiImqZZYljmyZPLulISbdKWi/p1A7HT5L0S0nXl9ufNBFnRMSkytlxq2yjorEah6Q5wGeB1wIbgGslrbR9S1vRr9k+eeABRkRUNUJJoYomaxyHAutt32H7ceCrwNIG44mI2Coar7aNiib7OPYC7mrZ3wAc1qHcsZJ+D7gN+DPbd7UXkLQMWAYwd8ddufwLh/ch3JnngF2OaDqEofGCizc1HcLQ0Etf1HQIw+Om3jzNKDVDVdFkjUMd7mt/+/8Z2Nf2gcB3gXM7PZHt5bYX21687TOe1eMwIyKm4BrbiGgycWwAFrbsLwA2thaw/YDtx8rdLwIvH1BsERHVJXEMzLXA/pKeL2k74HhgZWsBSfNbdo8B1g0wvoiIaU2MHM9VVQNge4ukk4FLgDnACts3SzoDWGN7JfA+SccAW4BNwElNxRsRMRmNj1BWqKDRAYC2VwGr2u47reX2h4APDTquiIjKRqwZqoqMHI+I6NIoNUNVkcQREdGtJI6IiKgjNY6IiKgniSMiIirzaE0nUkUSR0REF7ICYERE1OfZlTkaXY8jImIU9HLkeIV1ivaRdJmkGyRdKWlBy7GxlvWLVrY/tldS44iI6EYPBwBWXKfok8B5ts+V9J+BjwNvL489avug3kQzudQ4IiK61MP1OKqsU7QIuKy8fUWH432XxBER0aUaiWOepDUt27K2p+q0TtFebWXWAseWt38f2EnS7uX+DuXzXiXpTb19lU9KU1VERDdMnc7x+20vnuJ4lXWKTgHOknQSsBq4m2IiWIC9bW+UtB9wuaQbbf+0anBVJXFERHSph5fjVlmnaCPwZgBJOwLH2t7ccgzbd0i6EjgY6HniSFNVRES3ereQU5V1iuZJmvju/hCworx/V0nbT5QBXgm0dqr3TBJHREQXermQk+0twMQ6ReuAr0+sU1SuTQSwBLhV0m3AHsBfl/e/BFgjaS1Fp/mZbVdj9UyaqiIiumH3dCGnCusUXQRc1OFxPwR+p2eBTCGJIyKiW7Nr4HgSR0REtzJXVUREVGcga45HREQtsytvJHFERHQrTVUREVFLL6+qmgmSOCIiutHD2XFniiSOiIguFAMAZ1fmSOKIiOhW1hyPiIg6UuOIiIjqZmEfR6OTHFZYW3d7SV8rj18tad/BRxkRMZVirqoq26hoLHG0rK17FMVSiCdIWtRW7F3Ar2z/B+DvgE8MNsqIiArsatuIaLLGUWVt3aXAueXti4BXS+q0QlZERDPc0zXHZ4QmE0eVtXV/W6acp34zsDsREcNkltU4muwcr7K2bpUylAu+LwOYu+Ou3UcWEVHH6OSESpqscUy7tm5rGUnbAjsDm9qfyPZy24ttL972Gc/qU7gREZ1pfLzSNiqaTBzTrq1b7p9Y3n4LcLk9QvW9iJj5TDEAsMo2IhprqrK9RdLE2rpzgBUTa+sCa2yvBL4MnC9pPUVN4/im4o2I6EQ4AwAHqcLauv8OvHXQcUVE1JLEERERtSRxREREZRN9HLNIEkdERJdG6YqpKpI4IiK6MlqD+6pI4oiI6IZJ4oiIiJpmV0tVEkdERLcyjiMiIupJ4oiIiMpsGJtdbVVJHBER3UqNIyIiakniiIiIygyM0HriVSRxRER0xeD0cURERFUmneMREVHTLOvjaHIFwIiI0WBX2yqQdKSkWyWtl3Rqh+P7SLpM0g2SrpS0oOXYiZJuL7cT2x/bK0kcERFdqZg0KiQOSXOAzwJHAYuAEyQtaiv2SeA82wcCZwAfLx+7G3A6cBhwKHC6pF179jJbJHFERHTDwPh4tW16hwLrbd9h+3Hgq8DStjKLgMvK21e0HH89cKntTbZ/BVwKHNnty+skiSMiolvVaxzzJK1p2Za1PdNewF0t+xvK+1qtBY4tb/8+sJOk3Ss+tifSOR4R0ZVaU47cb3vxFMfV+QRPcQpwlqSTgNXA3cCWio/tiSSOiIhuGNy7cRwbgIUt+wuAjU85nb0ReDOApB2BY21vlrQBWNL22Ct7FVirNFVFRHRr3NW26V0L7C/p+ZK2A44HVrYWkDRP0sR394eAFeXtS4DXSdq17BR/XXlfzyVxRER0q0dXVdneApxM8YW/Dvi67ZslnSHpmLLYEuBWSbcBewB/XT52E/BRiuRzLXBGeV/PpakqIqIbdtUrpio+nVcBq9ruO63l9kXARZM8dgVP1kD6JokjIqJbs2zkeBJHRERXjMfGmg5ioJI4IiK6kWnVIyKitkyr3n/lnCpfA/YF7gT+oBwi315uDLix3P257WPay0RENMmAZ1mNo6nLcU8FLrO9P8WcK0+bAbL0qO2Dyi1JIyKGj8uFnKpsI6KpxLEUOLe8fS7wpobiiIjomsfGKm2jQm7gMjJJD9repWX/V7afNv2vpC3A9RTzsJxp+1uTPN8yYGKysAOAm3ofdW3zgPsTAzAccQxDDDAccQxDDDAccbzI9k7dPIGkf6V4LVXcb7svM9YOUt8Sh6TvAnt2OPQR4NyKieN5tjdK2g+4HHi17Z9Oc94100wiNhDDEMcwxDAscQxDDMMSxzDEMCxxDEMMM1HfOsdtv2ayY5LulTTf9j2S5gP3TfIcG8u/d0i6EjgYmDJxREREfzXVx7ESmFjW8ETg2+0Fyom6ti9vzwNeCdwysAgjIqKjphLHmcBrJd0OvLbcR9JiSV8qy7wEWCNpLcUqV2farpI4lvcj4K0wDHEMQwwwHHEMQwwwHHEMQwwwHHEMQwwzTiOd4xERMXNlWvWIiKgliSMiImqZ8YlD0m6SLpV0e/n3aZf1luXGJF1fbis7ldnK8x8p6VZJ6yU9bQS8pO0lfa08frWkfXt17hoxnCTply2v/0/6EMMKSfdJ6jiGRoVPlzHeIOmQBmJYImlzy/twWqdyPYhjoaQrJK2TdLOk93co09f3o2IMfX8/JO0g6RpJa8s4/keHMn39jFSMoe+fkZFie0ZvwP8ETi1vnwp8YpJyj/Th3HMoLg/eD9gOWAssaivzp8DZ5e3jga81EMNJwFl9/nf4PeAQ4KZJjh8NXAwIOBy4uoEYlgD/0s/3oTzPfOCQ8vZOwG0d/k36+n5UjKHv70f5+nYsb88FrgYObyvT789IlRj6/hkZpW3G1zhodvqSQ4H1tu+w/Tjw1TKeyeK7CHi1JA04hr6zvRqYapnKpcB5LlwF7FKO4RlkDANh+x7b15W3H6ZYAnSvtmJ9fT8qxtB35et7pNydW27tV+T09TNSMYaoYRQSxx6274HiwwI8d5JyO0haI+kqSb1KLnsBd7Xsb+DpH87flnGxnvBmYPcenb9qDADHlk0iF0la2MPzV1U1zn47omyyuFjSS/t9srLZ5WCKX7mtBvZ+TBEDDOD9kDRH0vUUA30vtT3pe9Gnz0iVGKD5z8iMMSMSh6TvSrqpw1bnl/XeLqYWeBvw95Je0IvQOtzX/kumSpl+x/DPwL62DwS+y5O/7gap3+9DFdcB+9h+GfAZoOPcZ70iaUfgG8AHbD/UfrjDQ3r+fkwTw0DeD9tjtg8CFgCHSjqgPcxODxtwDMPwGZkxZkTisP0a2wd02L4N3DtRxVfF6UuAKyl+gXVrA9D6y2QBsHGyMpK2BXamt80p08Zg+wHbj5W7XwRe3sPzV1Xlveor2w9NNFnYXgXMVTErQc9JmkvxhX2B7W92KNL392O6GAb5fpTneJDis9c+yV+/PyPTxjAkn5EZY0Ykjmk0OX3JtcD+kp4vaTuKjr32K7Za43sLcLntXv6amjaGtrbzYyjauwdtJfCO8mqiw4HNE02MgyJpz4m2c0mHUvz//0AfziPgy8A625+apFhf348qMQzi/ZD0HEm7lLefAbwG+Elbsb5+RqrEMCSfkZmj6d75bjeKttDLgNvLv7uV9y8GvlTefgXFSoJry7/v6uH5j6a4YuWnwEfK+84Ajilv7wD8I7AeuAbYrw/vwXQxfBy4uXz9VwAv7kMMFwL3AE9Q/IJ8F/Bu4N3lcQGfLWO8EVjcQAwnt7wPVwGv6NP/k6+iaGq5gWJZgOvLf6OBvR8VY+j7+wEcCPy4jOMm4LRBf0YqxtD3z8gobZlyJCIiahmFpqqIiBigJI6IiKgliSMiImpJ4oiIiFqSOCIiopYkjhgpkh6ZvlREdCOJIyIiakniiJFUjsj+X+WcZjdKOq68f4mkK8uJ7H4i6YIez1YcMfK2bTqAiD55M3AQ8DJgHnCtpNXlsYOBl1LMDfUDiilovt9EkBEzUWocMapeBVzoYlbUe4HvAb9bHrvG9gbb4xRTcezbUIwRM1ISR4yqqZqfHmu5PUZq3hG1JHHEqFoNHFcu4PMcimVlr2k4poiRkF9aMar+CTiCYrZTA39h+xeSXtxsWBEzX2bHjYiIWtJUFRERtSRxRERELUkcERFRSxJHRETUksQRERG1JHFEREQtSRwREVHL/wfSYcqbHbJuAQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "ds1.p0.plot.imshow()\n",
    "plt.title('Slope map')"
   ]
  },
  {
   "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.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}