hex_area.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import h3\n",
    "\n",
    "from math import sin, cos, asin, sqrt, tan, atan, pi\n",
    "\n",
    "def haversine(th1, ph1, th2, ph2):\n",
    "    \"\"\" Haversine between geo points.\n",
    "    \n",
    "    Input/output: both in radians.\n",
    "    \"\"\"\n",
    "    ph1 -= ph2\n",
    "\n",
    "    dz = sin(th1) - sin(th2)\n",
    "    dx = cos(ph1) * cos(th1) - cos(th2)\n",
    "    dy = sin(ph1) * cos(th1)\n",
    "\n",
    "    return asin(sqrt(dx*dx + dy*dy + dz*dz) / 2)*2\n",
    "\n",
    "\n",
    "def haversine_geo(g1, g2):\n",
    "    \"\"\" Haversine distance between two lat/lng points.\n",
    "    \n",
    "    input: degrees\n",
    "    output: radians\n",
    "    \"\"\"\n",
    "    th1, ph1 = g1\n",
    "    th2, ph2 = g2\n",
    "    \n",
    "    c = pi/180\n",
    "    th1 *= c\n",
    "    th2 *= c\n",
    "    ph1 *= c\n",
    "    ph2 *= c\n",
    "    \n",
    "    return haversine(th1, ph1, th2, ph2)\n",
    "    \n",
    "    \n",
    "def area_triangle(g1, g2, g3):\n",
    "    \"\"\" Surface area of a spherical triangle given by three lat/lng points\n",
    "    \n",
    "    input: degrees\n",
    "    ouput: km^2\n",
    "    \"\"\"\n",
    "    \n",
    "    a = haversine_geo(g1, g2)\n",
    "    b = haversine_geo(g2, g3)\n",
    "    c = haversine_geo(g3, g1)\n",
    "    \n",
    "    s = (a + b + c)/2\n",
    "    \n",
    "    T = sqrt(\n",
    "          tan(s/2)\n",
    "        * tan((s-a)/2)\n",
    "        * tan((s-b)/2)\n",
    "        * tan((s-c)/2)\n",
    "    )\n",
    "    \n",
    "    E = 4*atan(T)\n",
    "    \n",
    "    R = 6371.0088\n",
    "    A = (E*R)*R\n",
    "    \n",
    "    return A\n",
    "\n",
    "\n",
    "def cycle_pairs(it):\n",
    "    it = list(it)\n",
    "    it = it + [it[0]]\n",
    "    \n",
    "    return list(zip(it[:-1], it[1:]))\n",
    "\n",
    "def cell_to_triangles(h):\n",
    "    c = h3.h3_to_geo(h)\n",
    "    boundary = h3.h3_to_geo_boundary(h)\n",
    "    \n",
    "    pairs = cycle_pairs(boundary)\n",
    "    triples = [p + (c,)for p in pairs]\n",
    "    \n",
    "    return triples\n",
    "\n",
    "def cell_area(h):\n",
    "    triples = cell_to_triangles(h)\n",
    "\n",
    "    return sum(area_triangle(*t) for t in triples)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4357451.630819805"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Average area of res 0 *hexagons*\n",
    "import numpy as np\n",
    "\n",
    "cells = h3.get_res0_indexes()\n",
    "hexes = {\n",
    "    c\n",
    "    for c in cells\n",
    "    if not h3.h3_is_pentagon(c)\n",
    "}\n",
    "\n",
    "np.mean([cell_area(h) for h in hexes])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4250546.848"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Reported average area of res 0 *hexagons* (why the mismatch?)\n",
    "\n",
    "h3.hex_area(0)"
   ]
  },
  {
   "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}