elasticsearch-geo-aggregations.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Geospatial Elasticsearch Geo Aggregations with Python"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import json\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import pyproj\n",
    "import shapely.geometry\n",
    "import geopandas as gpd\n",
    "\n",
    "from pprint import pprint\n",
    "from dotenv import load_dotenv\n",
    "from elasticsearch import Elasticsearch\n",
    "\n",
    "load_dotenv(\".env\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#filepath = \"/home/hadoop/sgds/data/naturalearth/\" \\\n",
    "#           \"ne_10m_admin_0_countries.shp\"\n",
    "filepath = \"/home/hadoop/sgds/data/naturalearth/\" \\\n",
    "           \"ne_110m_admin_0_countries.shp\"\n",
    "\n",
    "gdf_ne = gpd.read_file(filepath)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "es = Elasticsearch(hosts=['http://{}:{}@{}:{}'.format(\n",
    "    os.environ[\"ES_USERNAME\"],\n",
    "    os.environ[\"ES_PASSWORD\"],\n",
    "    os.environ[\"ES_HOST\"],\n",
    "    os.environ[\"ES_PORT\"],\n",
    ")])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'count': 24330238,\n",
       " '_shards': {'total': 5, 'successful': 5, 'skipped': 0, 'failed': 0}}"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "es.count(index='gdelt_custom_2020')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Geohash Grid Aggregation\n",
    "\n",
    "- [Geohash Grid Aggregation](https://www.elastic.co/guide/en/elasticsearch/reference/current/search-aggregations-bucket-geohashgrid-aggregation.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "total 10000\n",
      "CPU times: user 0 ns, sys: 3.8 ms, total: 3.8 ms\n",
      "Wall time: 49.5 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "result = es.search(index='gdelt_custom_2020', body={\n",
    "    #\"size\": 10000,\n",
    "    \"size\": 0,\n",
    "    \"query\": {\n",
    "        \"bool\": {\n",
    "            \"must\": [],\n",
    "            \"must_not\": []\n",
    "        }\n",
    "    },\n",
    "    \"aggs\": {\n",
    "        \"scene-grid\": {\n",
    "            \"geohash_grid\": {\n",
    "                \"field\": \"location\",\n",
    "                \"precision\": 2\n",
    "            }\n",
    "        }\n",
    "    }\n",
    "})\n",
    "\n",
    "print('total', result['hits']['total']['value'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'pandas.core.frame.DataFrame'>\n",
      "RangeIndex: 293 entries, 0 to 292\n",
      "Data columns (total 2 columns):\n",
      " #   Column     Non-Null Count  Dtype \n",
      "---  ------     --------------  ----- \n",
      " 0   key        293 non-null    object\n",
      " 1   doc_count  293 non-null    int64 \n",
      "dtypes: int64(1), object(1)\n",
      "memory usage: 19.3 KB\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>key</th>\n",
       "      <th>doc_count</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>dr</td>\n",
       "      <td>6483</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>sv</td>\n",
       "      <td>5964</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>dq</td>\n",
       "      <td>4413</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>gc</td>\n",
       "      <td>4072</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>dp</td>\n",
       "      <td>3603</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  key  doc_count\n",
       "0  dr       6483\n",
       "1  sv       5964\n",
       "2  dq       4413\n",
       "3  gc       4072\n",
       "4  dp       3603"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df = pd.DataFrame(result['aggregations']['scene-grid']['buckets'])\n",
    "df.info(memory_usage='deep')\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "import geohash\n",
    "import shapely.geometry"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 6.88 ms, sys: 3.54 ms, total: 10.4 ms\n",
      "Wall time: 7.64 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "def get_geohash_geometry(gh_code):\n",
    "    coords = geohash.bbox(gh_code)\n",
    "    # (minx, miny, maxx, maxy)\n",
    "    extent = [coords['w'], coords['s'], \n",
    "              coords['e'], coords['n']]\n",
    "    return shapely.geometry.box(*extent)\n",
    "\n",
    "df['geom'] = df['key'].apply(get_geohash_geometry)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtEAAAGBCAYAAACgr+1tAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAe/UlEQVR4nO3df6yld10n8PenLTVdxAXsUPoDKNttWKqx1Z20VVih6mBpdMcudtPGYKOYCQZ2xcjG7pKoWbOJP1CJCjZVKyXRstuFkcYtP0Y0Ftad2ikppaX8KBWXcRo6gIIuTUztZ/+4Z+RyuXfu/fY+c8+Ze1+v5Oae8zzf53s+873fc+57nnnm+VZ3BwAA2LhT5l0AAACcbIRoAAAYJEQDAMAgIRoAAAYJ0QAAMEiIBgCAQSdtiK6qm6vq0aq6fwNtf62q7p19faKq/nYragQAYHuqk/U+0VX1nUn+PsnbuvubB477D0m+tbt/9IQVBwDAtnbSnonu7juTfGH5tqq6oKreU1X3VNUHqupfrXLodUlu3ZIiAQDYlk6bdwETuynJq7v7k1V1WZK3JPmuYzur6nlJnp/kT+ZUHwAA28C2CdFV9fVJviPJbVV1bPPXrWh2bZL/2d3/uJW1AQCwvWybEJ2lS1P+trsvOU6ba5O8ZovqAQBgmzppr4leqbu/lOQvq+qaJKklFx/bX1UvSPKMJP9nTiUCALBNnLQhuqpuzVIgfkFVHa6qVyX5oSSvqqoPJ3kgyd5lh1yX5O19st6OBACAhXHS3uIOAADm5aQ9Ew0AAPNyUv7HwjPPPLPPP//8eZcBAMA2d88993yuu3et3H5Shujzzz8/hw4dmncZAABsc1X1V6ttdzkHAAAMEqIBAGCQEA0AAIOEaAAAGCREAwDAoElCdFXdXFWPVtX9a+z/oaq6b/b15yuW4/50VX2kqu6tKrfcAABg4U11JvqtSa48zv6/TPKS7v6WJD+f5KYV+6/o7ku6e/dE9QAAwAkzyX2iu/vOqjr/OPv/fNnTg0nOm+J1AQBgHuZxTfSrkrx72fNO8r6quqeq9q11UFXtq6pDVXXo6NGjJ7xIAABYy5auWFhVV2QpRL942eYXdfeRqnpWkgNV9bHuvnPlsd19U2aXgezevbu3pGAAAFjFlp2JrqpvSfI7SfZ29+ePbe/uI7PvjybZn+TSraoJAACejC0J0VX13CTvTPLK7v7Esu1PraqnHXuc5GVJVr3DBwAALIpJLueoqluTvDTJmVV1OMnPJnlKknT3jUl+Jsk3JnlLVSXJ47M7cZyVZP9s22lJ/qC73zNFTQAAcKJMdXeO69bZ/2NJfmyV7Q8nufhrj2DEnlOumXcJJ9yBJ26bdwksoCnnvjk2ZlE/d3bCz9G852Sx3eeqFQsBAGCQEA0AAIOEaAAAGCREAwDAICEaAAAGCdEAADBIiAYAgEFCNAAADBKiAQBgkBANAACDhGgAABgkRAMAwCAhGgAABgnRAAAwSIgGAIBBQjQAAAwSogEAYJAQDQAAg4RoAAAYJEQDAMAgIRoAAAYJ0QAAMEiIBgCAQUI0AAAMEqIBAGCQEA0AAIOEaAAAGCREAwDAoNPmXQDb14Enbpt3CWxzU86xPadcM1lfU9oJf0bG+GydnynfQzvh57jd/4yTnImuqpur6tGqun+N/VVVv15VD1XVfVX1bcv2XVlVH5/tu2GKegAA4ESa6nKOtya58jj7X57kwtnXviS/lSRVdWqSN8/2X5Tkuqq6aKKaAADghJgkRHf3nUm+cJwme5O8rZccTPL0qjo7yaVJHuruh7v7H5K8fdYWAAAW1lb9x8Jzk3xm2fPDs21rbf8aVbWvqg5V1aGjR4+esEIBAGA9WxWia5VtfZztX7ux+6bu3t3du3ft2jVpcQAAMGKr7s5xOMlzlj0/L8mRJKevsR0AABbWVp2Jvj3JD8/u0nF5ki929yNJ7k5yYVU9v6pOT3LtrC0AACysSc5EV9WtSV6a5MyqOpzkZ5M8JUm6+8YkdyS5KslDSb6c5Edm+x6vqtcmeW+SU5Pc3N0PTFETAACcKJOE6O6+bp39neQ1a+y7I0shGwAATgqW/QYAgEFCNAAADBKiAQBgkBANAACDhGgAABgkRAMAwCAhGgAABgnRAAAwSIgGAIBBQjQAAAwSogEAYJAQDQAAg6q7513DsN27d/ehQ4fmXcam7Dnlmsn6euzqyybr64z9d03W1yI78MRtk/X14le8cbK+PviO10/W15RzbMrxYsyUP0e2hynfj4v6u2jKz0LYrKq6p7t3r9zuTDQAAAwSogEAYJAQDQAAg4RoAAAYJEQDAMAgIRoAAAYJ0QAAMEiIBgCAQUI0AAAMEqIBAGCQEA0AAIOEaAAAGCREAwDAICEaAAAGCdEAADBIiAYAgEGThOiqurKqPl5VD1XVDavs/09Vde/s6/6q+seqeuZs36er6iOzfYemqAcAAE6k0zbbQVWdmuTNSfYkOZzk7qq6vbs/eqxNd/9ykl+etf/+JD/Z3V9Y1s0V3f25zdYCAABbYdMhOsmlSR7q7oeTpKrenmRvko+u0f66JLdO8LrMfPAdr593Cavac8o18y5hTVPWdsZkPU1b14Enbpusr0Wtayc45+A3TNbXBw5eNFlf5/xZT9bXGfvvmqyvx66+bLK+pqxrSov6fpyyrkX9M8JyU1zOcW6Szyx7fni27WtU1T9LcmWSdyzb3EneV1X3VNW+tV6kqvZV1aGqOnT06NEJygYAgCdnihBdq2xb6xTF9yf53ysu5XhRd39bkpcneU1VfedqB3b3Td29u7t379q1a3MVAwDAJkwRog8nec6y5+clObJG22uz4lKO7j4y+/5okv1ZujwEAAAW1hQh+u4kF1bV86vq9CwF5dtXNqqqf57kJUnetWzbU6vqacceJ3lZkvsnqAkAAE6YTf/Hwu5+vKpem+S9SU5NcnN3P1BVr57tv3HW9Ook7+vu/7fs8LOS7K+qY7X8QXe/Z7M1AQDAiTTF3TnS3XckuWPFthtXPH9rkreu2PZwkounqAEAALaKFQsBAGCQEA0AAIOEaAAAGCREAwDAICEaAAAGCdEAADBIiAYAgEFCNAAADBKiAQBg0CQrFjLuwBO3zbuEk8qn3nT5pP1d8LqDk/X12NWXTdbXGfvvmqyvF7/ijZP1dcZkPSV7Trlmsr52wvvoyOVfmq6zN03X1ZSmfH9P+d6e0pRzdcr30JQW9c+4E8aL+XAmGgAABgnRAAAwSIgGAIBBQjQAAAwSogEAYJAQDQAAg4RoAAAYJEQDAMAgIRoAAAYJ0QAAMEiIBgCAQUI0AAAMEqIBAGDQafMuADbinD/rSft77OrLJuvrgp9+cLK+8tPfMFlXv/fcmybr63v3XzxZXweeuG2yvnaCSefq6w5O1teUdU39/l5Ee065Zt4lrGpR65rSp950+WR9TTlXpxx7n6vz4Uw0AAAMEqIBAGCQEA0AAIOEaAAAGCREAwDAoElCdFVdWVUfr6qHquqGVfa/tKq+WFX3zr5+ZqPHAgDAotn0Le6q6tQkb06yJ8nhJHdX1e3d/dEVTT/Q3d/3JI8FAICFMcWZ6EuTPNTdD3f3PyR5e5K9W3AsAADMxRQh+twkn1n2/PBs20rfXlUfrqp3V9U3DR6bqtpXVYeq6tDRo0cnKBsAAJ6cKUJ0rbJt5ZI+H0ryvO6+OMlvJPnDgWOXNnbf1N27u3v3rl27nnSxAACwWVOE6MNJnrPs+XlJjixv0N1f6u6/nz2+I8lTqurMjRwLAACLZooQfXeSC6vq+VV1epJrk9y+vEFVPbuqavb40tnrfn4jxwIAwKLZ9N05uvvxqnptkvcmOTXJzd39QFW9erb/xiQ/mOTHq+rxJI8luba7O8mqx262JhbDgSdum6yv6//iVZP1lSSf+sUXTtbXBw5eNFlfU/qRCfs68MTvTtgbI87Yf9e8S1jVotZ1zsFvmKyvKd/bF7zu4GR9PXb1ZZP1NeXPccrP/H/x678yWV9Tjv2n3nT5ZH3lJRP2xVxsOkQn/3SJxh0rtt247PFvJvnNjR4LAACLzIqFAAAwSIgGAIBBQjQAAAwSogEAYJAQDQAAg4RoAAAYJEQDAMAgIRoAAAYJ0QAAMEiIBgCAQZMs+w0n2qd+8YXzLmFN/+byj07W1wcOXjRZX7dc+ruT9bXnlGsm6+vAE7dN1hfz89jVl827hFX93nNvmqyvfznh+3HK8brgpx+crK8j+yfratLPiQsm62laF7zu4GR9+Sw8+TkTDQAAg4RoAAAYJEQDAMAgIRoAAAYJ0QAAMEiIBgCAQUI0AAAMEqIBAGCQEA0AAIOEaAAAGGTZb04KH3zH6yft78WveONkfR25/EuT9ZU3TdeVpbo5kc7Yf9e8S1jV9+6/eLK+zrm6J+trSp/6xRdO1teRN9VkfZ3zZ9ON15Tza8ol1xd13jMfzkQDAMAgIRoAAAYJ0QAAMEiIBgCAQUI0AAAMEqIBAGCQEA0AAIMmCdFVdWVVfbyqHqqqG1bZ/0NVdd/s68+r6uJl+z5dVR+pqnur6tAU9QAAwIm06cVWqurUJG9OsifJ4SR3V9Xt3f3RZc3+MslLuvtvqurlSW5Ksvzu51d09+c2WwsAAGyFKVYsvDTJQ939cJJU1duT7E3yTyG6u/98WfuDSc6b4HVPalOuJjelnbIy3aKuOnXB6w5O1teUq3RZ/XDMor6/p7SoP8cpx37KlVKnrOtTb7p8sr6m/MzZCRZ13jMfU1zOcW6Szyx7fni2bS2vSvLuZc87yfuq6p6q2rfWQVW1r6oOVdWho0ePbqpgAADYjCnORNcq23rVhlVXZClEv3jZ5hd195GqelaSA1X1se6+82s67L4pS5eBZPfu3av2DwAAW2GKM9GHkzxn2fPzkhxZ2aiqviXJ7yTZ292fP7a9u4/Mvj+aZH+WLg8BAICFNUWIvjvJhVX1/Ko6Pcm1SW5f3qCqnpvknUle2d2fWLb9qVX1tGOPk7wsyf0T1AQAACfMpi/n6O7Hq+q1Sd6b5NQkN3f3A1X16tn+G5P8TJJvTPKWqkqSx7t7d5KzkuyfbTstyR9093s2WxMAAJxIU1wTne6+I8kdK7bduOzxjyX5sVWOezjJxSu3AwDAIrNiIQAADBKiAQBgkBANAACDJrkmmnFWPZqvKcd/UVenW9TV1pifRV3Fcsq6zpisp8VdqfPFr3jjZH0tqinnxJQW9bNwUX+nbfes40w0AAAMEqIBAGCQEA0AAIOEaAAAGCREAwDAICEaAAAGCdEAADBIiAYAgEFCNAAADLJiIWzSoq4Utagray2qRR2vKVduO/KSmqyvC/ZP1lXO2H/XdJ0xZNLVIif8Oe6EObHdV/NLtv/qh85EAwDAICEaAAAGCdEAADBIiAYAgEFCNAAADBKiAQBgkBANAACDhGgAABgkRAMAwCArFgJkMVfDSqZd8WvKVQaZnylX81vUFVcXdSXFRf2cmNKizolF5Ew0AAAMEqIBAGCQEA0AAIOEaAAAGCREAwDAoElCdFVdWVUfr6qHquqGVfZXVf36bP99VfVtGz0WAAAWzaZDdFWdmuTNSV6e5KIk11XVRSuavTzJhbOvfUl+a+BYAABYKFOcib40yUPd/XB3/0OStyfZu6LN3iRv6yUHkzy9qs7e4LEAALBQplhs5dwkn1n2/HCSlXdJX63NuRs8NklSVfuydBY7z33uczdX8ZM05U3Dd8IN23eKRZ0Xi1rXlHbCn3FR69ruiygsskWdE4v6+bUTLOpn4aLO1alMcSa6VtnWG2yzkWOXNnbf1N27u3v3rl27BksEAIDpTHEm+nCS5yx7fl6SIxtsc/oGjgUAgIUyxZnou5NcWFXPr6rTk1yb5PYVbW5P8sOzu3RcnuSL3f3IBo8FAICFsukz0d39eFW9Nsl7k5ya5ObufqCqXj3bf2OSO5JcleShJF9O8iPHO3azNQEAwIk0xeUc6e47shSUl2+7cdnjTvKajR4LAACLzIqFAAAwSIgGAIBBQjQAAAwSogEAYNAk/7Fwp1jU1Zi2+4pAi24njL/Vw9iJFvW97ffH/PgsHLPd56oz0QAAMEiIBgCAQUI0AAAMEqIBAGCQEA0AAIOEaAAAGCREAwDAICEaAAAGWWxlgJuss1Mt6kJDsBMt4qITU1vUP+Oifn4t6nhtd85EAwDAICEaAAAGCdEAADBIiAYAgEFCNAAADBKiAQBgkBANAACDhGgAABgkRAMAwCArFs6J1YU40XbCKoPeR5xI5v32sKg/xylN+Wc0vzbOmWgAABgkRAMAwCAhGgAABgnRAAAwSIgGAIBBmwrRVfXMqjpQVZ+cfX/GKm2eU1V/WlUPVtUDVfUTy/b9XFX9dVXdO/u6ajP1AADAVtjsmegbkry/uy9M8v7Z85UeT/JT3f3CJJcneU1VXbRs/6919yWzrzs2WQ8AAJxwmw3Re5PcMnt8S5IfWNmgux/p7g/NHv9dkgeTnLvJ1wUAgLnZbIg+q7sfSZbCcpJnHa9xVZ2f5FuT3LVs82ur6r6qunm1y0GWHbuvqg5V1aGjR49usmwAAHjyqruP36Dqj5M8e5Vdb0hyS3c/fVnbv+nuVYNwVX19kj9L8t+6+52zbWcl+VySTvLzSc7u7h9dr+jdu3f3oUOH1mvGnFlBidUs6uph5tiYKX+Oj1192WR9ffAdr5+sL7aHRf3MmXLen7H/rvUbbZDPwq9VVfd09+6V29dd9ru7v+c4nX62qs7u7keq6uwkj67R7ilJ3pHk948F6Fnfn13W5reT/NF69QAAwLxt9nKO25NcP3t8fZJ3rWxQVZXkd5M82N2/umLf2cueXp3k/k3WAwAAJ9xmQ/QvJNlTVZ9Msmf2PFV1TlUdu9PGi5K8Msl3rXIru1+qqo9U1X1Jrkjyk5usBwAATrh1L+c4nu7+fJLvXmX7kSRXzR5/MEmtcfwrN/P6AAAwD1YsBACAQUI0AAAMEqIBAGCQEA0AAIOEaAAAGCREAwDAoE3d4g6Ox9Kh87WoS91OyRwDTlZTLtU9pSl/d2z3z2hnogEAYJAQDQAAg4RoAAAYJEQDAMAgIRoAAAYJ0QAAMEiIBgCAQUI0AAAMEqIBAGCQEA0AAIMs+w3b1JTLre6EJcSZnw++4/XzLgE2ZMrP1Re/4o2T9XXkJTVZXw//x5+arK/tvoS4M9EAADBIiAYAgEFCNAAADBKiAQBgkBANAACDhGgAABgkRAMAwCAhGgAABgnRAAAwSIgGAIBBQjQAAAw6bTMHV9Uzk/z3JOcn+XSSf9/df7NKu08n+bsk/5jk8e7ePXL8otjua8DDWsxXVjInOFks6lw9Y/9dk/V1wf7Jukr+43RdPXb1ZdN1toA2eyb6hiTv7+4Lk7x/9nwtV3T3JccC9JM4HgAAFsJmQ/TeJLfMHt+S5Ae2+HgAANhymw3RZ3X3I0ky+/6sNdp1kvdV1T1Vte9JHJ+q2ldVh6rq0NGjRzdZNgAAPHnrXhNdVX+c5Nmr7HrDwOu8qLuPVNWzkhyoqo91950Dx6e7b0pyU5Ls3r27R44FAIAprRuiu/t71tpXVZ+tqrO7+5GqOjvJo2v0cWT2/dGq2p/k0iR3JtnQ8QAAsEg2eznH7Umunz2+Psm7VjaoqqdW1dOOPU7ysiT3b/R4AABYNJsN0b+QZE9VfTLJntnzVNU5VXXHrM1ZST5YVR9O8hdJ/ld3v+d4xwMAwCLb1H2iu/vzSb57le1Hklw1e/xwkotHjgcAgEVmxUIAABgkRAMAwCAhGgAABgnRAAAwaFP/sXCnOfDEbfMuAQDYBhY1U+w55ZrJ+jpjsp4WkzPRAAAwSIgGAIBBQjQAAAwSogEAYJAQDQAAg4RoAAAYJEQDAMAgIRoAAAYJ0QAAMEiIBgCAQUI0AAAMEqIBAGCQEA0AAINOm3cBwImx55RrJuvrwBO3TdYXAIvL5/3GORMNAACDhGgAABgkRAMAwCAhGgAABgnRAAAwSIgGAIBBQjQAAAwSogEAYJAQDQAAg4RoAAAYtKkQXVXPrKoDVfXJ2fdnrNLmBVV177KvL1XV62b7fq6q/nrZvqs2Uw8AAGyFzZ6JviHJ+7v7wiTvnz3/Kt398e6+pLsvSfKvk3w5yf5lTX7t2P7uvmOT9QAAwAm32RC9N8kts8e3JPmBddp/d5JPdfdfbfJ1AQBgbjYbos/q7keSZPb9Weu0vzbJrSu2vbaq7quqm1e7HOSYqtpXVYeq6tDRo0c3VzUAAGzCuiG6qv64qu5f5WvvyAtV1elJ/m2S25Zt/q0kFyS5JMkjSX5lreO7+6bu3t3du3ft2jXy0gAAMKnT1mvQ3d+z1r6q+mxVnd3dj1TV2UkePU5XL0/yoe7+7LK+/+lxVf12kj/aWNnAeg48cdv6jQCAJ2Wzl3PcnuT62ePrk7zrOG2vy4pLOWbB+5irk9y/yXoAAOCE22yI/oUke6rqk0n2zJ6nqs6pqn+600ZV/bPZ/neuOP6XquojVXVfkiuS/OQm6wEAgBNu3cs5jqe7P5+lO26s3H4kyVXLnn85yTeu0u6Vm3l9AACYBysWAgDAICEaAAAGCdEAADBIiAYAgEFCNAAADBKiAQBgkBANAACDhGgAABgkRAMAwCAhGgAABgnRAAAwqLp73jUMq6qjSf5qTi9/ZpLPzem1T0bGa4zxGmO8xhivMcZrjPEaY7zGzHO8ntfdu1ZuPClD9DxV1aHu3j3vOk4WxmuM8RpjvMYYrzHGa4zxGmO8xizieLmcAwAABgnRAAAwSIged9O8CzjJGK8xxmuM8RpjvMYYrzHGa4zxGrNw4+WaaAAAGORMNAAADBKiAQBgkBC9jqr65ar6WFXdV1X7q+rpa7S7sqo+XlUPVdUNW13noqiqa6rqgap6oqrWvBVNVX26qj5SVfdW1aGtrHGRDIyX+ZWkqp5ZVQeq6pOz789Yo92OnV/rzZVa8uuz/fdV1bfNo85FsYHxemlVfXE2l+6tqp+ZR52LoqpurqpHq+r+NfabX8tsYLzMr5mqek5V/WlVPTj7vfgTq7RZqPklRK/vQJJv7u5vSfKJJP95ZYOqOjXJm5O8PMlFSa6rqou2tMrFcX+Sf5fkzg20vaK7L1m0+z5usXXHy/z6KjckeX93X5jk/bPna9lx82uDc+XlSS6cfe1L8ltbWuQCGXhvfWA2ly7p7v+6pUUunrcmufI4+82vr/bWHH+8EvPrmMeT/FR3vzDJ5Ules+ifX0L0Orr7fd39+OzpwSTnrdLs0iQPdffD3f0PSd6eZO9W1bhIuvvB7v74vOs4WWxwvMyvr9ib5JbZ41uS/MAca1lEG5kre5O8rZccTPL0qjp7qwtdEN5bg7r7ziRfOE4T82uZDYwXM939SHd/aPb475I8mOTcFc0Wan4J0WN+NMm7V9l+bpLPLHt+OF/7g+erdZL3VdU9VbVv3sUsOPPrK87q7keSpQ/cJM9ao91OnV8bmSvm01dsdCy+vao+XFXvrqpv2prSTlrm1zjza4WqOj/Jtya5a8WuhZpfp83rhRdJVf1xkmevsusN3f2uWZs3ZOmfGn5/tS5W2bZt7x24kfHagBd195GqelaSA1X1sdnf2LedCcbL/FryhoFudsz8WmEjc2VHzad1bGQsPpTked3991V1VZI/zNI/JbM682uM+bVCVX19knckeV13f2nl7lUOmdv8EqKTdPf3HG9/VV2f5PuSfHevfmPtw0mes+z5eUmOTFfhYllvvDbYx5HZ90eran+W/ll1W4acCcbL/Jqpqs9W1dnd/cjsn/AeXaOPHTO/VtjIXNlR82kd647F8l/i3X1HVb2lqs7s7s9tUY0nG/NrgPn11arqKVkK0L/f3e9cpclCzS+Xc6yjqq5M8tNJ/m13f3mNZncnubCqnl9Vpye5NsntW1XjyaaqnlpVTzv2OMnLsvQf7Fid+fUVtye5fvb4+iRfcyZ/h8+vjcyV25P88Ox/uV+e5IvHLpHZgdYdr6p6dlXV7PGlWfq9+fktr/TkYX4NML++YjYOv5vkwe7+1TWaLdT8ciZ6fb+Z5Ouy9E/CSXKwu19dVeck+Z3uvqq7H6+q1yZ5b5JTk9zc3Q/Mr+T5qaqrk/xGkl1J/ldV3dvd37t8vJKclWT/bDxPS/IH3f2euRU9RxsZL/Prq/xCkv9RVa9K8n+TXJMk5teSteZKVb16tv/GJHckuSrJQ0m+nORH5lXvvG1wvH4wyY9X1eNJHkty7Rr/IrkjVNWtSV6a5MyqOpzkZ5M8JTG/VrOB8TK/vuJFSV6Z5CNVde9s239J8txkMeeXZb8BAGCQyzkAAGCQEA0AAIOEaAAAGCREAwDAICEaAIBtp6purqpHq2rd25xW1a9V1b2zr09U1d+ue4y7cwAAsN1U1Xcm+fskb+vubx447j8k+dbu/tHjtXMmGgCAbae770zyheXbquqCqnpPVd1TVR+oqn+1yqHXJbl1vf4ttgIAwE5xU5JXd/cnq+qyJG9J8l3HdlbV85I8P8mfrNeREA0AwLZXVV+f5DuS3DZb1TZZWpV6uWuT/M/u/sf1+hOiAQDYCU5J8rfdfclx2lyb5DUb7QwAALa17v5Skr+sqmuSpJZcfGx/Vb0gyTOS/J+N9CdEAwCw7VTVrVkKxC+oqsNV9aokP5TkVVX14SQPJNm77JDrkry9N3jrOre4AwCAQc5EAwDAICEaAAAGCdEAADBIiAYAgEFCNAAADBKiAQBgkBANAACD/j9S7jW/NEhVeAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x864 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "gdf = gpd.GeoDataFrame(df, geometry='geom', crs='epsg:4326')\n",
    "\n",
    "mask = gdf.geom.apply(lambda g: g.centroid.y > -80)\n",
    "\n",
    "gdf[mask].to_crs(epsg=3857).plot(figsize=(12, 12), \n",
    "                                 column='doc_count', \n",
    "                                 cmap='viridis', \n",
    "                                 scheme='naturalbreaks', k=5);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# GeoTile Grid Aggregation\n",
    "\n",
    "> A multi-bucket aggregation that works on `geo_point` fields and groups points into buckets that represent cells in a grid. The resulting grid can be sparse and only contains cells that have matching data. Each cell corresponds to a [map tile](https://en.wikipedia.org/wiki/Tiled_web_map) as used by many online map sites. Each cell is labeled using a \"{zoom}/{x}/{y}\" format, where zoom is equal to the user-specified precision.\n",
    "\n",
    "- [GeoTile Grid Aggregation](https://www.elastic.co/guide/en/elasticsearch/reference/current/search-aggregations-bucket-geotilegrid-aggregation.html)\n",
    "- OpenStreetMap Wiki - [Slippy Map](https://wiki.openstreetmap.org/wiki/Slippy_Map)\n",
    "- OpenStreetMap Wiki - [Slippy map tilenames](https://wiki.openstreetmap.org/wiki/Slippy_map_tilenames)\n",
    "- Github - [mapbox/mercantile](https://github.com/mapbox/mercantile)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "total 10000\n",
      "CPU times: user 8.1 ms, sys: 0 ns, total: 8.1 ms\n",
      "Wall time: 208 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "EXTENT = [-126.6,23.1,-64.4,49.3]\n",
    "\n",
    "result = es.search(index='gdelt_custom_2020', body={\n",
    "    #\"size\": 10000,\n",
    "    \"size\": 0,\n",
    "    \"query\": {\n",
    "        \"bool\": {\n",
    "            \"must\": [],\n",
    "            \"must_not\": [],\n",
    "            #\"filter\" : {\n",
    "            #    \"geo_bounding_box\" : {\n",
    "            #        \"location\" : {\n",
    "            #            \"top_left\" : [EXTENT[0], EXTENT[3]],\n",
    "            #            \"bottom_right\" : [EXTENT[2], EXTENT[1]]\n",
    "            #        }\n",
    "            #    }\n",
    "            #}\n",
    "        }\n",
    "    },\n",
    "    \"aggs\": {\n",
    "        \"scene-grid\": {\n",
    "            \"geotile_grid\": {\n",
    "                \"field\": \"location\",\n",
    "                \"precision\": 6\n",
    "            }\n",
    "        }\n",
    "    }\n",
    "})\n",
    "\n",
    "print('total', result['hits']['total']['value'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'pandas.core.frame.DataFrame'>\n",
      "RangeIndex: 512 entries, 0 to 511\n",
      "Data columns (total 2 columns):\n",
      " #   Column     Non-Null Count  Dtype \n",
      "---  ------     --------------  ----- \n",
      " 0   key        512 non-null    object\n",
      " 1   doc_count  512 non-null    int64 \n",
      "dtypes: int64(1), object(1)\n",
      "memory usage: 36.1 KB\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>key</th>\n",
       "      <th>doc_count</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>6/18/24</td>\n",
       "      <td>5709</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>6/39/25</td>\n",
       "      <td>3796</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>6/41/25</td>\n",
       "      <td>3118</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>6/18/23</td>\n",
       "      <td>2707</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>6/19/23</td>\n",
       "      <td>2365</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       key  doc_count\n",
       "0  6/18/24       5709\n",
       "1  6/39/25       3796\n",
       "2  6/41/25       3118\n",
       "3  6/18/23       2707\n",
       "4  6/19/23       2365"
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df = pd.DataFrame(result['aggregations']['scene-grid']['buckets'])\n",
    "df.info(memory_usage='deep')\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [],
   "source": [
    "df[['z', 'x', 'y']] = df['key'].apply(\n",
    "    lambda k: pd.Series([int(k) for k in k.split('/')], \n",
    "                        index=['z', 'x', 'y']))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [],
   "source": [
    "def xyz_tile_to_bbox(x, y, z):\n",
    "    n = 2**z\n",
    "    lon_deg0 = x / n * 360.0 - 180.0\n",
    "    lat_rad0 = np.arctan(np.sinh(np.pi * (1 - 2 * y / n)))\n",
    "    lat_deg0 = lat_rad0 * 180.0 / np.pi\n",
    "    lon_deg1 = (x + 1) / n * 360.0 - 180.0\n",
    "    lat_rad1 = np.arctan(np.sinh(np.pi * (1 - 2 * (y + 1) / n)))\n",
    "    lat_deg1 = lat_rad1 * 180.0 / np.pi\n",
    "    \n",
    "    # [W, S, E, N]\n",
    "    extent = [lon_deg0, lat_deg1, lon_deg1, lat_deg0]\n",
    "    return shapely.geometry.box(*extent)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['geom'] = df.apply(\n",
    "    lambda row: xyz_tile_to_bbox(row['x'], row['y'], row['z']), axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtEAAAF7CAYAAAAUptcSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAf70lEQVR4nO3df4xlZ3kf8O8TG1cpoQXihdiYxdSySGgUO8nIuAUFnMTUWG03NEGyFRE3Aa2ocFtHoYrbSEnU/pM2SROlARw3sTBSYtQUXKzWARxaxdDUqdfIgB0DNgaX7Vo4/HRTkJDTp3/M3Xg8zOyc13N27p3dz0ca7b3nvOfcd955752vX585T3V3AACA6b5l2R0AAID9RogGAIBBQjQAAAwSogEAYJAQDQAAg4RoAAAYtG9DdFXdVFWPVdV9E9r+WlXdu/j6VFV9ZS/6CADAqan2632iq+oHkvx5knd293cPHPePk3xvd//USescAACntH27Et3ddyb50sZtVXVBVb2vqu6pqg9V1XducejVSW7Zk04CAHBKOnPZHZjZjUne1N0PVtXLkrwtyQ8e31lVL0ry4iT/dUn9AwDgFHDKhOiq+rYkfzvJ71fV8c1/ZVOzq5L8x+7+i73sGwAAp5ZTJkRn/dKUr3T3xSdoc1WSN+9RfwAAOEXt22uiN+vux5N8pqpelyS17qLj+6vqJUmek+R/LKmLAACcIvZtiK6qW7IeiF9SVUer6g1JfjzJG6rqo0nuT3JowyFXJ3lX79fbkQAAsDL27S3uAABgWfbtSjQAACzLvvzDwrPPPrvPP//8ZXcDAIBT3D333POF7j6wefu+DNHnn39+jhw5suxuAABwiquqR7ba7nIOAAAYNEuIrqqbquqxqrpvm/0/XlUfW3z98aZbz322qj5eVfdWleVlAABW3lwr0e9IcsUJ9n8mySu7+3uS/Kusl+fe6LLuvri712bqDwAAnDSzXBPd3XdW1fkn2P/HG57eleS8OV4XAACWYRnXRL8hyR9seN5JPlBV91TV4e0OqqrDVXWkqo782Z/92UnvJAAAbGdP785RVZdlPUS/YsPml3f3sap6XpI7quoT3X3n5mO7+8YsLgNZW1tTIQYAgKXZs5XoqvqeJL+d5FB3f/H49u4+tvj3sSS3Jrlkr/oEAABPx56E6Ko6mOQ9SV7f3Z/asP2ZVfWs44+TvDrJlnf4AACAVTHL5RxVdUuSVyU5u6qOJvmFJM9Iku6+IcnPJ/n2JG+rqiR5YnEnjucnuXWx7cwkv9fd75ujTwAAcLLMdXeOq3fY/8Ykb9xi+8NJLvrmIwAAYHWpWAgAAIOEaAAAGCREAwDAoD29TzTAnC7/ltft2OaO//f7e9ATAE43VqIBAGCQEA0AAIOEaAAAGCREAwDAICEaAAAGCdEAADBIiAYAgEFCNAAADBKiAQBgkIqFp4ApVdsSlds49ZjTyzNntUiVJ8cYL/aLUz2fWIkGAIBBQjQAAAwSogEAYJAQDQAAg4RoAAAYJEQDAMAgIRoAAAYJ0QAAMEixlVPAMm5SPvUG6lPs15usA/PwGTDGeMFqsBINAACDhGgAABgkRAMAwCAhGgAABgnRAAAwSIgGAIBBs4Toqrqpqh6rqvu22V9V9RtV9VBVfayqvm/Dviuq6pOLfdfP0R8AADiZ5lqJfkeSK06w/zVJLlx8HU7y9iSpqjOSvHWx/6VJrq6ql87UJwAAOClmCdHdfWeSL52gyaEk7+x1dyV5dlWdk+SSJA9198Pd/Y0k71q0BQCAlbVXFQtfkORzG54fXWzbavvLtjpBVR3O+ip2Dh48eHJ6eYqbUmVwaiUsFbNge3O+11bVfu8/nEynw2fAFKf697hXf1hYW2zrE2z/5o3dN3b3WnevHThwYNbOAQDAiL1aiT6a5IUbnp+X5FiSs7bZDgAAK2uvVqJvS/ITi7t0XJrkq939aJK7k1xYVS+uqrOSXLVoCwAAK2uWleiquiXJq5KcXVVHk/xCkmckSXffkOT2JFcmeSjJ15L85GLfE1V1bZL3JzkjyU3dff8cfQIAgJNllhDd3VfvsL+TvHmbfbdnPWQDAMC+oGIhAAAMEqIBAGCQEA0AAIP26hZ3rIA5b3q+1zeSn/J6c7/m6cC4jpk6XvvZqn6P5iCwaqxEAwDAICEaAAAGCdEAADBIiAYAgEFCNAAADBKiAQBgkBANAACDhGgAABgkRAMAwCAVC1fcXlcPUxVs3Ct+9Fd2bPPhd79lD3oyzs97zKqO1+lQeXLOz0LVVJdnyudlsrqfmVP5eZ8erEQDAMAgIRoAAAYJ0QAAMEiIBgCAQUI0AAAMEqIBAGCQEA0AAIOEaAAAGFTdvew+DFtbW+sjR44suxsrYxk3+1/VwgerSkEGnq69fq8t472910Wlptrrz8xlvP9XtV+wSqrqnu5e27zdSjQAAAwSogEAYJAQDQAAg4RoAAAYJEQDAMAgIRoAAAbNEqKr6oqq+mRVPVRV12+x/59V1b2Lr/uq6i+q6rmLfZ+tqo8v9rlvHQAAK+/M3Z6gqs5I8tYklyc5muTuqrqtu//0eJvu/uUkv7xo//eS/HR3f2nDaS7r7i/sti8AALAX5liJviTJQ939cHd/I8m7khw6Qfurk9wyw+sCAMBS7HolOskLknxuw/OjSV62VcOq+qtJrkhy7YbNneQDVdVJfqu7b9zm2MNJDifJwYMHZ+j2cs1ZwW4Z1aTmrHy2qtXK5rSMapF7PS9WtV/73ad//dLZzrXX77VlvLe//totf/08xYff/ZZJ5/LZBJzIHCvRtcW27WqJ/70k/33TpRwv7+7vS/KaJG+uqh/Y6sDuvrG717p77cCBA7vrMQAA7MIcIfpokhdueH5ekmPbtL0qmy7l6O5ji38fS3Jr1i8PAQCAlTVHiL47yYVV9eKqOivrQfm2zY2q6q8neWWS927Y9syqetbxx0leneS+GfoEAAAnza6vie7uJ6rq2iTvT3JGkpu6+/6qetNi/w2Lpq9N8oHu/r8bDn9+klur6nhffq+737fbPgEAwMk0xx8WprtvT3L7pm03bHr+jiTv2LTt4SQXzdEHAADYKyoWAgDAICEaAAAGCdEAADBolmuiGXc63OB+6vd4OhQ0mGpVx2LOfk051+nw/pjqguvu2rHN1IIsUwqRTPWtt/7JbOeaYur3OGW85iwMtIwiQ6v6Htrrfs35ueQzh6fDSjQAAAwSogEAYJAQDQAAg4RoAAAYJEQDAMAgIRoAAAYJ0QAAMEiIBgCAQUI0AAAMqu5edh+Gra2t9ZEjR5bdDWYyperU1Epre11FbZVNGbNVHS/Vw8ZMrdw2perfuX803++EVZ1fczodKgMuw6p+j8uoUMnyVdU93b22ebuVaAAAGCREAwDAICEaAAAGCdEAADBIiAYAgEFCNAAADBKiAQBgkBANAACDzlx2B2BOyyjKMvU1pzj2ytqxzdRiGBf87AM7N/rZvzbpXMcufXzHNooLnBxTiztMccF1d812rjnnPevm/FnPea6p7+05C6RMaTdn4ZM5x2vq+ZYxrszLSjQAAAwSogEAYJAQDQAAg4RoAAAYJEQDAMAgIRoAAAbNEqKr6oqq+mRVPVRV12+x/1VV9dWqunfx9fNTjwUAgFWz6/tEV9UZSd6a5PIkR5PcXVW3dfefbmr6oe7+u0/zWAAAWBlzrERfkuSh7n64u7+R5F1JDu3BsQAAsBRzVCx8QZLPbXh+NMlWpaz+VlV9NMmxJG/p7vsHjk1VHU5yOEkOHjw4Q7fZT6ZU8kumV2Sbs9ralMqAn/7X3zXb600dizsPfmjHNn/n3IsmnUs1rNW2jEqdc55rijnfs3P2fe5Kd3PZ7+/ZOav0/Y3f+NUd21ww6UzJp3/90kntpvwuWtW5w3RzrERv9Rt9c13ijyR5UXdflOTfJflPA8eub+y+sbvXunvtwIEDT7uzAACwW3OE6KNJXrjh+XlZX23+S939eHf/+eLx7UmeUVVnTzkWAABWzRwh+u4kF1bVi6vqrCRXJbltY4Oq+o6qqsXjSxav+8UpxwIAwKrZ9TXR3f1EVV2b5P1JzkhyU3ffX1VvWuy/IcmPJflHVfVEkq8nuaq7O8mWx+62TwAAcDLN8YeFxy/RuH3Tths2PP7NJL859VgAAFhlKhYCAMAgIRoAAAYJ0QAAMKjW/75vf1lbW+sjR44suxvsoSk3y0+Sc/9ovvk8tSDDnEUgphZSmWLKzf6nFiuYs/ABY+YsyDDlZzT19abM+70uyJIk597112Y717FLH5/t9W6+5Hd2bDPn2M9ZoGoZ7+1JBVImFtdaVVMLtzz8T37mJPeEnVTVPd29tnm7lWgAABgkRAMAwCAhGgAABgnRAAAwSIgGAIBBQjQAAAwSogEAYJAQDQAAg4RoAAAYdOayOwBTTK3YdPl181V3m2pKZbCplRSnVOCaWuVqzipjc1a6U9lwfsuoPLnXVSxf8aO/MrHlAzu2mFKJMJlYGfDSaVUZL8/efzZNMeV7nLNq5lQX7Pkr7j2VCPc/K9EAADBIiAYAgEFCNAAADBKiAQBgkBANAACDhGgAABgkRAMAwCAhGgAABlX3tCIQq2Rtba2PHDmy7G7sib0uaHC6mDKuU4oQJMm33rpzsYU5C1gso18sz37+DJizAM+cBT+mjtf0Ai87m/J+XIapnydzWcY4zPk9ztn/VX3f8s2q6p7uXtu83Uo0AAAMEqIBAGCQEA0AAIOEaAAAGCREAwDAICEaAAAGzRKiq+qKqvpkVT1UVddvsf/Hq+pji68/rqqLNuz7bFV9vKrurarT4751AADsa2fu9gRVdUaStya5PMnRJHdX1W3d/acbmn0mySu7+8tV9ZokNybZeOPGy7r7C7vtCwAA7IU5VqIvSfJQdz/c3d9I8q4khzY26O4/7u4vL57eleS8GV4XAACWYtcVC6vqx5Jc0d1vXDx/fZKXdfe127R/S5Lv3ND+M0m+nKST/FZ337jNcYeTHE6SgwcPfv8jjzyyq37vF8uo0sW6OceeJ5mH65ZRzW+vz7UMq/qZOWe/plTgO/bKmnSuC667a7fdWXmnw7zn5NquYuGuL+dIstU7dctkXlWXJXlDklds2Pzy7j5WVc9LckdVfaK77/ymE66H6xuT9bLfu+82AAA8PXNcznE0yQs3PD8vybHNjarqe5L8dpJD3f3F49u7+9ji38eS3Jr1y0MAAGBlzRGi705yYVW9uKrOSnJVkts2Nqiqg0nek+T13f2pDdufWVXPOv44yauT3DdDnwAA4KTZ9eUc3f1EVV2b5P1JzkhyU3ffX1VvWuy/IcnPJ/n2JG+rqiR5YnFtyfOT3LrYdmaS3+vu9+22TwAAcDLNcU10uvv2JLdv2nbDhsdvTPLGLY57OMlFm7cDAMAqU7EQAAAGCdEAADBIiAYAgEG7LrayDGtra33kyJFldwMmW9UCEHNSrOBJ+7lQz5RCHknyrbf+yUnuyf6wqvP5FT/6K5PaTfk5Tv0ezfsn7XWhnlUtDLSq749R2xVbsRINAACDhGgAABgkRAMAwCAhGgAABgnRAAAwSIgGAIBBQjQAAAwSogEAYJAQDQAAg1QshBWynyt+TbWfK1jNWZVxzp/11Gprx15Zs73mBdfdtWObVa2iNqdlzOcpYzF1Tkzx4Xe/ZbZzLaManrlzcpxOlQ1VLAQAgJkI0QAAMEiIBgCAQUI0AAAMEqIBAGCQEA0AAIOEaAAAGCREAwDAoDOX3QHgSXtdpGOqVe3XXpuzIMAyCpFccOvObeYs0jFncRrGzFkgZRmmzMNV/Sw8Xfi9YCUaAACGCdEAADBIiAYAgEFCNAAADBKiAQBgkBANAACDZgnRVXVFVX2yqh6qquu32F9V9RuL/R+rqu+beiwAAKyaXYfoqjojyVuTvCbJS5NcXVUv3dTsNUkuXHwdTvL2gWMBAGClzLESfUmSh7r74e7+RpJ3JTm0qc2hJO/sdXcleXZVnTPxWAAAWClzVCx8QZLPbXh+NMnmUkNbtXnBxGOTJFV1OOur2Dl48ODuenySTanQo+rR6WXOym1T585eV4paxpz2Xlu338d+VauazVmRbVWrXe61OefEfp9fc1rV6qD7ea5OMcdKdG2xrSe2mXLs+sbuG7t7rbvXDhw4MNhFAACYzxwr0UeTvHDD8/OSHJvY5qwJxwIAwEqZYyX67iQXVtWLq+qsJFcluW1Tm9uS/MTiLh2XJvlqdz868VgAAFgpu16J7u4nquraJO9PckaSm7r7/qp602L/DUluT3JlkoeSfC3JT57o2N32CQAATqY5LudId9+e9aC8cdsNGx53kjdPPRYAAFaZioUAADBIiAYAgEFCNAAADJrlmmieas6biysmcWrY7z+jvS58sN/Hi9PHfi90s58to4jKMoqaKBazbhXntJVoAAAYJEQDAMAgIRoAAAYJ0QAAMEiIBgCAQUI0AAAMEqIBAGCQEA0AAIOEaAAAGKRi4UmwXyvvwHbmnK9TzjVnVbBlVBhj3apWWtvvP+v93v8p9vv3uNc5YL+P135lJRoAAAYJ0QAAMEiIBgCAQUI0AAAMEqIBAGCQEA0AAIOEaAAAGCREAwDAIMVWBiyjcIAbqPN0zTl39nMBoVXt1+lg6tivalGWve6Xubo8qzoHp1JUajmsRAMAwCAhGgAABgnRAAAwSIgGAIBBQjQAAAwSogEAYNCuQnRVPbeq7qiqBxf/PmeLNi+sqv9WVQ9U1f1V9U837PvFqvrfVXXv4uvK3fQHAAD2wm5Xoq9P8sHuvjDJBxfPN3siyc9093cluTTJm6vqpRv2/1p3X7z4un2X/QEAgJNutyH6UJKbF49vTvIjmxt096Pd/ZHF4/+T5IEkL9jl6wIAwNJUdz/9g6u+0t3P3vD8y939TZd0bNh/fpI7k3x3dz9eVb+Y5B8meTzJkayvWH95m2MPJzmcJAcPHvz+Rx555Gn3m72hghJbUQXu9KHKK6tgVasRfv21L5vU7ltv/ZOT3JOn8h76ZlV1T3evbd6+40p0Vf1hVd23xdehwQ58W5J3J7muux9fbH57kguSXJzk0SS/ut3x3X1jd69199qBAwdGXhoAAGZ15k4NuvuHt9tXVZ+vqnO6+9GqOifJY9u0e0bWA/Tvdvd7Npz78xva/Psk/3mk8wAAsAy7vSb6tiTXLB5fk+S9mxtUVSX5nSQPdPe/3bTvnA1PX5vkvl32BwAATrrdhuhfSnJ5VT2Y5PLF81TVuVV1/E4bL0/y+iQ/uMWt7P5NVX28qj6W5LIkP73L/gAAwEm34+UcJ9LdX0zyQ1tsP5bkysXjDyepbY5//W5eHwAAlkHFQgAAGCREAwDAICEaAAAG7eqaaDgRN2xffatahMDcAU51U4uoTPk8nPOzXKG06axEAwDAICEaAAAGCdEAADBIiAYAgEFCNAAADBKiAQBgkBANAACDhGgAABgkRAMAwCAVCwEAJpizSt/UyoBT2n361y+ddK4LrrtrxzbL+B6nWMUKiVaiAQBgkBANAACDhGgAABgkRAMAwCAhGgAABgnRAAAwSIgGAIBBQjQAAAwSogEAYJCKhXAaW0ZlqlWsOsXJ4WfNKtjv8/Drr33Zjm3O/aOedK79Pharxko0AAAMEqIBAGCQEA0AAIOEaAAAGCREAwDAoF2F6Kp6blXdUVUPLv59zjbtPltVH6+qe6vqyOjxAACwSna7En19kg9294VJPrh4vp3Luvvi7l57mscDAMBK2G2IPpTk5sXjm5P8yB4fDwAAe263xVae392PJkl3P1pVz9umXSf5QFV1kt/q7hsHj09VHU5yOEkOHjy4y24/PYpJwPbMe4Dppn5mTs0e7L0dQ3RV/WGS79hi188NvM7Lu/vYIiTfUVWf6O47B47PInjfmCRra2vTSvMAAMBJsGOI7u4f3m5fVX2+qs5ZrCKfk+Sxbc5xbPHvY1V1a5JLktyZZNLxAACwSnZ7TfRtSa5ZPL4myXs3N6iqZ1bVs44/TvLqJPdNPR4AAFbNbkP0LyW5vKoeTHL54nmq6tyqun3R5vlJPlxVH03yP5P8l+5+34mOBwCAVbarPyzs7i8m+aEtth9LcuXi8cNJLho5HgAAVpmKhQAAMEiIBgCAQUI0AAAMEqIBAGDQbisWnlZUZAMA9tKqZo85Kymu6ve4EyvRAAAwSIgGAIBBQjQAAAwSogEAYJAQDQAAg4RoAAAYJEQDAMAgIRoAAAYptgIAwJD9WiBlTlaiAQBgkBANAACDhGgAABgkRAMAwCAhGgAABgnRAAAwSIgGAIBBQjQAAAwSogEAYJCKhXAau/xbXjepncpUAPBUVqIBAGCQEA0AAIOEaAAAGCREAwDAICEaAAAG7SpEV9Vzq+qOqnpw8e9ztmjzkqq6d8PX41V13WLfL1bV/96w78rd9AcAAPbCbleir0/ywe6+MMkHF8+fors/2d0Xd/fFSb4/ydeS3Lqhya8d39/dt++yPwAAcNLtNkQfSnLz4vHNSX5kh/Y/lOTT3f3ILl8XAACWZrfFVp7f3Y8mSXc/WlXP26H9VUlu2bTt2qr6iSRHkvxMd395qwOr6nCSw0ly8ODB3fUaSKKICgA8XTuuRFfVH1bVfVt8HRp5oao6K8nfT7Lxt/bbk1yQ5OIkjyb51e2O7+4bu3utu9cOHDgw8tIAADCrHVeiu/uHt9tXVZ+vqnMWq9DnJHnsBKd6TZKPdPfnN5z7Lx9X1b9P8p+ndRsAAJZnt9dE35bkmsXja5K89wRtr86mSzkWwfu41ya5b5f9AQCAk263IfqXklxeVQ8muXzxPFV1blX95Z02quqvLva/Z9Px/6aqPl5VH0tyWZKf3mV/AADgpNvVHxZ29xezfseNzduPJblyw/OvJfn2Ldq9fjevDwAAy6BiIQAADBKiAQBgkBANAACDhGgAABgkRAMAwCAhGgAABgnRAAAwSIgGAIBBQjQAAAyq7l52H4ZV1Z8leWRJL392ki8s6bX3I+M1xniNMV5jjNcY4zXGeI0xXmOWOV4v6u4DmzfuyxC9TFV1pLvXlt2P/cJ4jTFeY4zXGOM1xniNMV5jjNeYVRwvl3MAAMAgIRoAAAYJ0eNuXHYH9hnjNcZ4jTFeY4zXGOM1xniNMV5jVm68XBMNAACDrEQDAMAgIRoAAAYJ0Tuoql+uqk9U1ceq6taqevY27a6oqk9W1UNVdf1e93NVVNXrqur+qvp/VbXtrWiq6rNV9fGqureqjuxlH1fJwHiZX0mq6rlVdUdVPbj49znbtDtt59dOc6XW/cZi/8eq6vuW0c9VMWG8XlVVX13MpXur6ueX0c9VUVU3VdVjVXXfNvvNrw0mjJf5tVBVL6yq/1ZVDyx+L/7TLdqs1PwSond2R5Lv7u7vSfKpJP98c4OqOiPJW5O8JslLk1xdVS/d016ujvuS/IMkd05oe1l3X7xq933cYzuOl/n1FNcn+WB3X5jkg4vn2znt5tfEufKaJBcuvg4nefuednKFDLy3PrSYSxd397/c006unnckueIE+82vp3pHTjxeifl13BNJfqa7vyvJpUnevOqfX0L0Drr7A939xOLpXUnO26LZJUke6u6Hu/sbSd6V5NBe9XGVdPcD3f3JZfdjv5g4XubXkw4luXnx+OYkP7LEvqyiKXPlUJJ39rq7kjy7qs7Z646uCO+tQd19Z5IvnaCJ+bXBhPFiobsf7e6PLB7/nyQPJHnBpmYrNb+E6DE/leQPttj+giSf2/D8aL75B89TdZIPVNU9VXV42Z1ZcebXk57f3Y8m6x+4SZ63TbvTdX5NmSvm05OmjsXfqqqPVtUfVNXf3Juu7Vvm1zjza5OqOj/J9yb5k027Vmp+nbmsF14lVfWHSb5ji10/193vXbT5uaz/r4bf3eoUW2w7Ze8dOGW8Jnh5dx+rqucluaOqPrH4L/ZTzgzjZX6t+7mB05w282uTKXPltJpPO5gyFh9J8qLu/vOqujLJf8r6/0pma+bXGPNrk6r6tiTvTnJddz++efcWhyxtfgnRSbr7h0+0v6quSfJ3k/xQb31j7aNJXrjh+XlJjs3Xw9Wy03hNPMexxb+PVdWtWf/fqqdkyJlhvMyvhar6fFWd092PLv4X3mPbnOO0mV+bTJkrp9V82sGOY7Hxl3h3315Vb6uqs7v7C3vUx/3G/Bpgfj1VVT0j6wH6d7v7PVs0Wan55XKOHVTVFUl+Nsnf7+6vbdPs7iQXVtWLq+qsJFcluW2v+rjfVNUzq+pZxx8neXXW/8COrZlfT7otyTWLx9ck+aaV/NN8fk2ZK7cl+YnFX7lfmuSrxy+ROQ3tOF5V9R1VVYvHl2T99+YX97yn+4f5NcD8etJiHH4nyQPd/W+3abZS88tK9M5+M8lfyfr/Ek6Su7r7TVV1bpLf7u4ru/uJqro2yfuTnJHkpu6+f3ldXp6qem2Sf5fkQJL/UlX3dvff2TheSZ6f5NbFeJ6Z5Pe6+31L6/QSTRkv8+spfinJf6iqNyT5X0lelyTm17rt5kpVvWmx/4Yktye5MslDSb6W5CeX1d9lmzheP5bkH1XVE0m+nuSqbf6P5Gmhqm5J8qokZ1fV0SS/kOQZifm1lQnjZX496eVJXp/k41V172Lbv0hyMFnN+aXsNwAADHI5BwAADBKiAQBgkBANAACDhGgAABgkRAMAcMqpqpuq6rGq2vE2p1X1a1V17+LrU1X1lR2PcXcOAABONVX1A0n+PMk7u/u7B477x0m+t7t/6kTtrEQDAHDK6e47k3xp47aquqCq3ldV91TVh6rqO7c49Ookt+x0fsVWAAA4XdyY5E3d/WBVvSzJ25L84PGdVfWiJC9O8l93OpEQDQDAKa+qvi3J307y+4uqtsl6VeqNrkryH7v7L3Y6nxANAMDp4FuSfKW7Lz5Bm6uSvHnqyQAA4JTW3Y8n+UxVvS5Jat1Fx/dX1UuSPCfJ/5hyPiEaAIBTTlXdkvVA/JKqOlpVb0jy40neUFUfTXJ/kkMbDrk6ybt64q3r3OIOAAAGWYkGAIBBQjQAAAwSogEAYJAQDQAAg4RoAAAYJEQDAMAgIRoAAAb9f40djpuFby+BAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x864 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "gdf = gpd.GeoDataFrame(df, geometry='geom', crs='epsg:4326')\n",
    "\n",
    "mask = gdf.geom.apply(lambda g: g.centroid.y > -80)\n",
    "\n",
    "gdf[mask].to_crs(epsg=3857).plot(figsize=(12, 12), \n",
    "                                 column='doc_count', \n",
    "                                 cmap='viridis', \n",
    "                                 scheme='naturalbreaks', k=5);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Geo Bounds Aggregation\n",
    "\n",
    "> A metric aggregation that computes the bounding box containing all geo values for a field.\n",
    "\n",
    "- [Geo Bounds Aggregation](https://www.elastic.co/guide/en/elasticsearch/reference/current/search-aggregations-metrics-geobounds-aggregation.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "total 10000\n",
      "CPU times: user 3.82 ms, sys: 396 µs, total: 4.21 ms\n",
      "Wall time: 359 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "result = es.search(index='gdelt_custom_2020', body={\n",
    "    #\"size\": 10000,\n",
    "    \"size\": 0,\n",
    "    \"aggs\": {\n",
    "        \"viewport\": {\n",
    "            \"geo_bounds\": {\n",
    "                \"field\": \"location\",    \n",
    "                \"wrap_longitude\": True  \n",
    "            }\n",
    "        }\n",
    "    }\n",
    "})\n",
    "\n",
    "print('total', result['hits']['total']['value'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'viewport': {'bounds': {'bottom_right': {'lat': -85.62210003845394,\n",
      "                                          'lon': 179.33299999684095},\n",
      "                         'top_left': {'lat': 71.9999999832362,\n",
      "                                      'lon': -176.5330000128597}}}}\n"
     ]
    }
   ],
   "source": [
    "pprint(result['aggregations'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Geo Centroid Aggregation\n",
    "\n",
    "> A metric aggregation that computes the weighted [centroid](https://en.wikipedia.org/wiki/Centroid) from all coordinate values for a [Geo-point datatype](https://www.elastic.co/guide/en/elasticsearch/reference/6.8/geo-point.html) field.\n",
    "\n",
    "- [Geo Centroid Aggregation](https://www.elastic.co/guide/en/elasticsearch/reference/6.8/search-aggregations-metrics-geocentroid-aggregation.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "total 10000\n",
      "CPU times: user 4.18 ms, sys: 0 ns, total: 4.18 ms\n",
      "Wall time: 216 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "result = es.search(index='gdelt_custom_2020', body={\n",
    "    #\"size\": 10000,\n",
    "    \"size\": 0,\n",
    "    \"aggs\": {\n",
    "        \"centroid\": {\n",
    "            \"geo_centroid\": {\n",
    "                \"field\": \"location\",    \n",
    "            }\n",
    "        }\n",
    "    }\n",
    "})\n",
    "\n",
    "print('total', result['hits']['total']['value'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'centroid': {'count': 97148,\n",
      "              'location': {'lat': 30.46608442130391,\n",
      "                           'lon': 0.35541025303929935}}}\n"
     ]
    }
   ],
   "source": [
    "pprint(result['aggregations'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Geo Distance Aggregation\n",
    "\n",
    "> A multi-bucket aggregation that works on `geo_point` fields and conceptually works very similar to the [range](https://www.elastic.co/guide/en/elasticsearch/reference/current/search-aggregations-bucket-range-aggregation.html) aggregation. The user can define a point of origin and a set of distance range buckets. The aggregation evaluate the distance of each document value from the origin point and determines the buckets it belongs to based on the ranges \n",
    "\n",
    "- [Geo Distance Aggregation](https://www.elastic.co/guide/en/elasticsearch/reference/current/search-aggregations-bucket-geodistance-aggregation.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "total 10000\n",
      "CPU times: user 6.03 ms, sys: 204 µs, total: 6.24 ms\n",
      "Wall time: 27.5 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "lat, lon = 48.864716, 2.349014\n",
    "\n",
    "result = es.search(index='gdelt_custom_2020', body={\n",
    "    #\"size\": 10000,\n",
    "    \"size\": 0,\n",
    "    \"aggs\": {\n",
    "        \"rings\": {\n",
    "            \"geo_distance\": {\n",
    "                \"field\": \"location\",    \n",
    "                \"origin\": {\"lat\": lat, \"lon\": lon},\n",
    "                \"unit\": \"km\",\n",
    "                \"ranges\": [\n",
    "                    {\"to\":   1000},\n",
    "                    {\"from\": 1000, \"to\": 2000},\n",
    "                    {\"from\": 2000}\n",
    "                ]\n",
    "            }\n",
    "        }\n",
    "    }\n",
    "})\n",
    "\n",
    "print('total', result['hits']['total']['value'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'rings': {'buckets': [{'doc_count': 7492,\n",
      "                        'from': 0.0,\n",
      "                        'key': '*-1000.0',\n",
      "                        'to': 1000.0},\n",
      "                       {'doc_count': 3015,\n",
      "                        'from': 1000.0,\n",
      "                        'key': '1000.0-2000.0',\n",
      "                        'to': 2000.0},\n",
      "                       {'doc_count': 86641,\n",
      "                        'from': 2000.0,\n",
      "                        'key': '2000.0-*'}]}}\n"
     ]
    }
   ],
   "source": [
    "pprint(result['aggregations'])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "gis",
   "language": "python",
   "name": "gis"
  },
  "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}