cythonized_geopandas_performance.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import sys\n",
    "sys.path.insert(0, '../../')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import geopandas as gpd\n",
    "import numpy as np\n",
    "import shapely\n",
    "import xarray as xr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from emiprep.regrid.cdo import _metgrid_to_cdo_grid_info_extraction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<xarray.Dataset>\n",
       "Dimensions:          (grid_corners: 4, grid_size: 44100, grid_xsize: 210, grid_ysize: 210)\n",
       "Coordinates:\n",
       "  * grid_size        (grid_size) int64 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 ...\n",
       "Dimensions without coordinates: grid_corners, grid_xsize, grid_ysize\n",
       "Data variables:\n",
       "    grid_center_lat  (grid_ysize, grid_xsize) float32 36.3869 36.4172 ...\n",
       "    grid_center_lon  (grid_ysize, grid_xsize) float32 -12.1618 -12.0042 ...\n",
       "    grid_corner_lat  (grid_ysize, grid_xsize, grid_corners) float32 36.3083 ...\n",
       "    grid_corner_lon  (grid_ysize, grid_xsize, grid_corners) float32 -12.2217 ...\n",
       "    grid_imask       (grid_size) int64 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ...\n",
       "    dummydata        (grid_ysize, grid_xsize) float32 0.0 0.0 0.0 0.0 0.0 ..."
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ds = _metgrid_to_cdo_grid_info_extraction('../../tests/testdata/metgrid_coords.nc')\n",
    "ds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def mk_geoms(ds):\n",
    "    lat_crnr = ds['grid_corner_lat'].values.reshape(-1, 4)\n",
    "    lon_crnr = ds['grid_corner_lon'].values.reshape(-1, 4)\n",
    "    geoms = []\n",
    "    for lons, lats in zip(lon_crnr, lat_crnr):\n",
    "        geoms.append(shapely.geometry.Polygon(((x, y) for x, y in zip(lons, lats))).convex_hull)\n",
    "    return geoms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 36.30832291,  36.33865356,  36.46554184,  36.4351387 ], dtype=float32)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ds['grid_corner_lat'].values.reshape(-1, 4)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 36.30832291,  36.33865356,  36.46554184,  36.4351387 ], dtype=float32)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ds['grid_corner_lat'].values[0,0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "geoms = mk_geoms(ds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.4 s ± 95.8 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%timeit geoms = mk_geoms(ds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0        POLYGON ((-12.2216796875 36.30832290649414, -1...\n",
       "1        POLYGON ((-12.06427001953125 36.33865356445312...\n",
       "2        POLYGON ((-11.90667724609375 36.36869812011719...\n",
       "3        POLYGON ((-11.74893188476562 36.39848709106445...\n",
       "4        POLYGON ((-11.59103393554688 36.4279899597168,...\n",
       "5        POLYGON ((-11.4329833984375 36.45722961425781,...\n",
       "6        POLYGON ((-11.2747802734375 36.48619842529297,...\n",
       "7        POLYGON ((-11.11639404296875 36.51488494873047...\n",
       "8        POLYGON ((-10.95785522460938 36.54330444335938...\n",
       "9        POLYGON ((-10.7991943359375 36.57144546508789,...\n",
       "10       POLYGON ((-10.64035034179688 36.59931182861328...\n",
       "11       POLYGON ((-10.48138427734375 36.62689971923828...\n",
       "12       POLYGON ((-10.32223510742188 36.65422058105469...\n",
       "13       POLYGON ((-10.1629638671875 36.68126678466797,...\n",
       "14       POLYGON ((-10.0035400390625 36.7080192565918, ...\n",
       "15       POLYGON ((-9.843963623046875 36.7344970703125,...\n",
       "16       POLYGON ((-9.684234619140625 36.76070785522461...\n",
       "17       POLYGON ((-9.524383544921875 36.78664398193359...\n",
       "18       POLYGON ((-9.3643798828125 36.81229019165039, ...\n",
       "19       POLYGON ((-9.2042236328125 36.83766174316406, ...\n",
       "20       POLYGON ((-9.0439453125 36.86275482177734, -9....\n",
       "21       POLYGON ((-8.883514404296875 36.8875617980957,...\n",
       "22       POLYGON ((-8.72296142578125 36.91208648681641,...\n",
       "23       POLYGON ((-8.562286376953125 36.93633651733398...\n",
       "24       POLYGON ((-8.401458740234375 36.96030044555664...\n",
       "25       POLYGON ((-8.240509033203125 36.98399353027344...\n",
       "26       POLYGON ((-8.07940673828125 37.00739288330078,...\n",
       "27       POLYGON ((-7.918212890625 37.03050994873047, -...\n",
       "28       POLYGON ((-7.756866455078125 37.0533447265625,...\n",
       "29       POLYGON ((-7.59539794921875 37.07590103149414,...\n",
       "                               ...                        \n",
       "44070    POLYGON ((28.32315063476562 64.23828125, 28.03...\n",
       "44071    POLYGON ((28.60980224609375 64.19645690917969,...\n",
       "44072    POLYGON ((28.89572143554688 64.15412902832031,...\n",
       "44073    POLYGON ((29.18084716796875 64.11131286621094,...\n",
       "44074    POLYGON ((29.4652099609375 64.06800079345703, ...\n",
       "44075    POLYGON ((29.74880981445312 64.02419281005859,...\n",
       "44076    POLYGON ((30.0316162109375 63.97989273071289, ...\n",
       "44077    POLYGON ((30.31362915039062 63.93510055541992,...\n",
       "44078    POLYGON ((30.5948486328125 63.88982772827148, ...\n",
       "44079    POLYGON ((30.87527465820312 63.84406280517578,...\n",
       "44080    POLYGON ((31.1549072265625 63.79781723022461, ...\n",
       "44081    POLYGON ((31.4337158203125 63.75108337402344, ...\n",
       "44082    POLYGON ((31.71170043945312 63.70387649536133,...\n",
       "44083    POLYGON ((31.9888916015625 63.65618133544922, ...\n",
       "44084    POLYGON ((32.2652587890625 63.60801696777344, ...\n",
       "44085    POLYGON ((32.540771484375 63.55937957763672, 3...\n",
       "44086    POLYGON ((32.81549072265625 63.51026153564453,...\n",
       "44087    POLYGON ((33.08935546875 63.46067047119141, 32...\n",
       "44088    POLYGON ((33.36239624023438 63.41061401367188,...\n",
       "44089    POLYGON ((33.63458251953125 63.36008453369141,...\n",
       "44090    POLYGON ((33.90591430664062 63.30909729003906,...\n",
       "44091    POLYGON ((34.17642211914062 63.25764465332031,...\n",
       "44092    POLYGON ((34.446044921875 63.20572662353516, 3...\n",
       "44093    POLYGON ((34.71484375 63.15335083007812, 34.44...\n",
       "44094    POLYGON ((34.98275756835938 63.10051345825195,...\n",
       "44095    POLYGON ((35.24981689453125 63.04721832275391,...\n",
       "44096    POLYGON ((35.51602172851562 62.99347686767578,...\n",
       "44097    POLYGON ((35.78134155273438 62.93927764892578,...\n",
       "44098    POLYGON ((36.04580688476562 62.88462829589844,...\n",
       "44099    POLYGON ((36.30938720703125 62.82952880859375,...\n",
       "Name: geometry, Length: 44100, dtype: object"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "grid = gpd.GeoDataFrame(geometry=geoms)\n",
    "grid['geometry']\n",
    "cell0 = grid['geometry'].loc[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "RangeIndex(start=0, stop=44100, step=1)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "grid.index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "POINT (-12.16183306218986 36.38693706684736)\n"
     ]
    }
   ],
   "source": [
    "print(cell0.centroid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def mk_randomcoords(ds, size=20000):\n",
    "    np.random.seed(123)\n",
    "    random_lons = np.random.uniform(ds['grid_corner_lon'].values.min(), ds['grid_corner_lon'].values.max(), size)\n",
    "    np.random.seed(321)\n",
    "    random_lats = np.random.uniform(ds['grid_corner_lat'].values.min(), ds['grid_corner_lat'].values.max(), size)\n",
    "    return random_lons, random_lats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "rlons, rlats = mk_randomcoords(ds, 20000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "points = [shapely.geometry.Point(x, y) for x, y in zip(rlons, rlats)]\n",
    "rgeoms = gpd.GeoDataFrame(geometry=points)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "315 ms ± 13.9 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%timeit rgeoms = [shapely.geometry.Point(x, y) for x, y in zip(rlons, rlats)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "joined = gpd.sjoin(rgeoms, grid, op='within')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'POLYGON ((19.13558959960938 57.14743804931641, 18.8924560546875 57.17330169677734, 18.93991088867188 57.30516815185547, 19.183837890625 57.27922439575195, 19.13558959960938 57.14743804931641))'"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "grid['geometry'].loc[31451].wkt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'POINT (18.90326593723728 60.21152971829886)'"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rgeoms['geometry'].loc[1998].wkt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "grid['geometry'].loc[31451].contains(rgeoms['geometry'].loc[1998])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "41 s ± 5.38 s per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%timeit _ = gpd.sjoin(rgeoms, grid, op='within')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'0.2+124.g4db21e7'"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gpd.__version__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'/home2/hilboll/.local/easybuild/software/Anaconda3/4.2.0/envs/emiprep_dev/lib/python3.6/site-packages/geopandas-0.2+124.g4db21e7-py3.6-linux-x86_64.egg/geopandas/__init__.py'"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gpd.__file__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((20000, 1), (44100, 1))"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rgeoms.shape, grid.shape"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}