Reference cycleway query for a given polygon area

example_notebook.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "verbal-burton",
   "metadata": {},
   "outputs": [],
   "source": [
    "import requests\n",
    "from shapely.geometry import shape\n",
    "from shapely.ops import unary_union\n",
    "\n",
    "# query oakland boundary\n",
    "resp = requests.get(\"https://raw.githubusercontent.com/maxogden/oakland-boundaries-geojson/master/districts.geojson\")\n",
    "oak_bounds = unary_union([shape(f[\"geometry\"]) for f in resp.json()[\"features\"]])\n",
    "poly_str = ' '.join([' '.join([str(x) for x in c[::-1]]) for c in oak_bounds.__geo_interface__[\"coordinates\"][0]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "afraid-afternoon",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Query:\n",
      "[out:json][timeout:65];(way[\"highway\"=\"cycleway\"](poly:'37.757221 -122.12024 37.754154 -122.11842 37.749592 -122.1161 37.740585 -122.11472 37.73999 -122.11498 37.739231 -122.11556 37.73671 -122.11536 37.734005 -122.11677 37.732559 -122.11798 37.730892 -122.12099 37.729176 -122.12367 37.729099 -122.12469 37.729149 -122.12481 37.729275 -122.12513 37.729351 -122.12543 37.729477 -122.126419 37.729588 -122.12687 37.730106 -122.13058 37.73142 -122.13178 37.731491 -122.13147 37.731869 -122.13103 37.732269 -122.13074 37.73277 -122.13029 37.73312 -122.13042 37.733624 -122.130363 37.733788 -122.130638 37.733688 -122.131676 37.733788 -122.13203 37.734177 -122.13407 37.734001 -122.13496 37.732494 -122.13608 37.731602 -122.13715 37.731091 -122.13798 37.73094 -122.13837 37.7315 -122.13849 37.731625 -122.13822 37.73198 -122.13845 37.731972 -122.13847 37.73333 -122.13888 37.734688 -122.1393 37.73445 -122.140205 37.73534 -122.140465 37.737324 -122.14152 37.737843 -122.14233 37.73805 -122.142487 37.738186 -122.14259 37.7384 -122.14282 37.739807 -122.1441 37.740276 -122.14452 37.74112 -122.14525 37.741768 -122.14584 37.7425 -122.14652 37.741104 -122.14777 37.74086 -122.14873 37.740623 -122.14954 37.740177 -122.151154 37.739082 -122.155121 37.738106 -122.15854 37.737667 -122.16016 37.735939 -122.16416 37.73723 -122.165131 37.73769 -122.16547 37.737839 -122.16568 37.73506 -122.17156 37.732536 -122.16997 37.732162 -122.169937 37.731735 -122.17093 37.731213 -122.174324 37.7306 -122.1736 37.73061 -122.17337 37.72827 -122.170685 37.72778 -122.17012 37.727856 -122.17083 37.727898 -122.17126 37.727074 -122.17246 37.726673 -122.17339 37.726646 -122.17385 37.72633 -122.17463 37.72544 -122.17529 37.725155 -122.1762 37.725185 -122.1767 37.725231 -122.17728 37.72509 -122.17896 37.725159 -122.17929 37.725384 -122.17961 37.72562 -122.18026 37.725319 -122.18157 37.725391 -122.18294 37.725666 -122.18389 37.725815 -122.1841 37.726322 -122.18516 37.726551 -122.185562 37.726791 -122.18578 37.727856 -122.18731 37.728497 -122.18911 37.72882 -122.19022 37.727974 -122.19093 37.728065 -122.19112 37.727097 -122.19179 37.727795 -122.19311 37.72604 -122.19632 37.726952 -122.19685 37.726479 -122.19797 37.724991 -122.19743 37.725212 -122.19688 37.723606 -122.19611 37.72262 -122.19488 37.722164 -122.194283 37.720531 -122.19228 37.71928 -122.19071 37.71791 -122.189026 37.71781 -122.18914 37.716793 -122.19273 37.718628 -122.19355 37.717823 -122.19562 37.715618 -122.19501 37.715611 -122.19706 37.71309 -122.20181 37.710827 -122.20927 37.702206 -122.20515 37.691132 -122.19973 37.681744 -122.19402 37.683102 -122.1923 37.669247 -122.18118 37.669308 -122.17827 37.66292 -122.18514 37.656223 -122.19217 37.632225 -122.21914 37.633579 -122.22025 37.64321 -122.22817 37.649578 -122.233398 37.672432 -122.25198 37.688889 -122.265648 37.708229 -122.28178 37.707623 -122.2554 37.710789 -122.25151 37.717014 -122.26151 37.72781 -122.24882 37.71862 -122.2361 37.726738 -122.22818 37.72687 -122.22806 37.726555 -122.22631 37.73037 -122.22664 37.73146 -122.22578 37.73307 -122.22552 37.745255 -122.22512 37.745346 -122.22578 37.74559 -122.22697 37.74826 -122.22643 37.754147 -122.22394 37.754429 -122.22378 37.760242 -122.22415 37.762402 -122.224663 37.764397 -122.22515 37.766674 -122.22694 37.768921 -122.23059 37.771732 -122.23656 37.777 -122.24438 37.777355 -122.244881 37.77912 -122.24482 37.780926 -122.24498 37.78342 -122.24618 37.786125 -122.25058 37.78673 -122.25478 37.785427 -122.25668 37.78423 -122.25788 37.788029 -122.26838 37.791649 -122.27463 37.792427 -122.27628 37.79355 -122.28268 37.793827 -122.284279 37.792526 -122.29478 37.794727 -122.314682 37.800629 -122.34028 37.811428 -122.34688 37.815968 -122.34855 37.835728 -122.35588 37.84111 -122.33014 37.837158 -122.3152 37.832355 -122.29433 37.831341 -122.29301 37.830193 -122.29102 37.829445 -122.28766 37.82877 -122.28556 37.827976 -122.28278 37.827633 -122.28126 37.827618 -122.28077 37.82725 -122.27959 37.827187 -122.27872 37.827034 -122.27812 37.82704 -122.27805 37.828735 -122.27857 37.831287 -122.277626 37.83646 -122.27563 37.83649 -122.27615 37.836727 -122.27601 37.836662 -122.27618 37.836723 -122.276413 37.837147 -122.27959 37.836857 -122.28186 37.837151 -122.28213 37.83672 -122.28453 37.8367 -122.28462 37.839779 -122.28558 37.842003 -122.28629 37.84402 -122.28694 37.845772 -122.2875 37.84691 -122.28786 37.849148 -122.288506 37.849972 -122.28892 37.850426 -122.28622 37.850132 -122.28599 37.850628 -122.28411 37.848976 -122.28355 37.849182 -122.28265 37.846973 -122.281799 37.847511 -122.27957 37.847698 -122.27675 37.845959 -122.27624 37.846581 -122.27182 37.846931 -122.26916 37.84808 -122.26961 37.848606 -122.26969 37.848717 -122.26841 37.849987 -122.2686 37.851608 -122.2689 37.851799 -122.26705 37.85299 -122.258331 37.85321 -122.25615 37.85358 -122.25336 37.850883 -122.25305 37.85117 -122.25126 37.851448 -122.24975 37.85181 -122.24614 37.852016 -122.24193 37.852409 -122.23631 37.852547 -122.23451 37.852543 -122.23418 37.855347 -122.23447 37.85749 -122.23494 37.857597 -122.23712 37.857311 -122.23788 37.857182 -122.239357 37.857361 -122.24181 37.85759 -122.24225 37.858738 -122.24335 37.858597 -122.24392 37.859001 -122.24402 37.861496 -122.24416 37.863701 -122.24464 37.864784 -122.244766 37.86646 -122.24486 37.868385 -122.24507 37.870483 -122.24528 37.87443 -122.2455 37.87646 -122.24577 37.879192 -122.24618 37.8802 -122.2463 37.88098 -122.24635 37.881268 -122.2464 37.884769 -122.24656 37.885254 -122.24667 37.885071 -122.24643 37.881924 -122.24198 37.88292 -122.23988 37.8831 -122.23935 37.883297 -122.23857 37.882214 -122.23434 37.88132 -122.23098 37.879227 -122.22578 37.878906 -122.2251 37.878326 -122.22388 37.875725 -122.22171 37.871723 -122.21738 37.86882 -122.216278 37.867867 -122.2179 37.867226 -122.21918 37.865025 -122.221489 37.86443 -122.22039 37.859509 -122.21474 37.857574 -122.21313 37.8568 -122.21185 37.854946 -122.21081 37.854527 -122.21098 37.85169 -122.20849 37.851555 -122.20737 37.85139 -122.20409 37.848682 -122.20174 37.84774 -122.200691 37.847588 -122.20068 37.844826 -122.19678 37.844677 -122.19659 37.84436 -122.19621 37.84359 -122.19527 37.842789 -122.19531 37.842007 -122.1961 37.84159 -122.19521 37.84147 -122.1946 37.841286 -122.19419 37.84116 -122.19412 37.840759 -122.19392 37.839832 -122.19179 37.839626 -122.19158 37.838829 -122.19063 37.838558 -122.19031 37.837627 -122.18887 37.836922 -122.18523 37.835892 -122.18471 37.83392 -122.18418 37.83254 -122.18514 37.832298 -122.18527 37.83211 -122.18537 37.831814 -122.18559 37.831356 -122.18589 37.831013 -122.18616 37.829063 -122.18713 37.827145 -122.18697 37.827065 -122.18681 37.826878 -122.18672 37.825527 -122.18587 37.82483 -122.18614 37.823524 -122.186676 37.820724 -122.18597 37.819527 -122.18148 37.816227 -122.176979 37.814327 -122.17508 37.810429 -122.17448 37.80819 -122.17542 37.80724 -122.17599 37.804722 -122.17773 37.804012 -122.1782 37.803455 -122.17753 37.803314 -122.17736 37.803696 -122.17658 37.80316 -122.17528 37.802235 -122.17569 37.7997 -122.17081 37.799026 -122.168709 37.80065 -122.1685 37.801163 -122.1686 37.801071 -122.169395 37.802315 -122.17019 37.804733 -122.17059 37.80559 -122.16668 37.80511 -122.16478 37.80442 -122.16447 37.804394 -122.16473 37.80323 -122.16427 37.801357 -122.16192 37.800896 -122.16196 37.80074 -122.16223 37.799934 -122.16187 37.799736 -122.16008 37.79814 -122.158577 37.797451 -122.15928 37.796955 -122.158577 37.79736 -122.15787 37.796867 -122.15688 37.796471 -122.15688 37.796394 -122.15147 37.79242 -122.15117 37.79234 -122.146271 37.790379 -122.14408 37.790943 -122.143387 37.790131 -122.1422 37.78767 -122.14406 37.78574 -122.14226 37.778027 -122.13508 37.778328 -122.13438 37.777527 -122.13367 37.77713 -122.13428 37.775063 -122.13236 37.7718 -122.13007 37.769699 -122.128166 37.767605 -122.1262 37.757221 -122.12024');>;);out;\n"
     ]
    }
   ],
   "source": [
    "# assemble the query string/url\n",
    "base_endpoint = \"https://overpass-api.de/api\"\n",
    "url = base_endpoint.rstrip(\"/\") + \"/interpreter\"\n",
    "overpass_settings = \"[out:json][timeout:65]\"\n",
    "osm_filter = 'way[\"highway\"=\"cycleway\"]'\n",
    "polygon_coord_str = \"37.786453 -122.327957 37.786453 -122.239380 37.840699 -122.239380 37.840699 -122.327957 37.786453 -122.327957\"\n",
    "query_str = f\"{overpass_settings};({osm_filter}(poly:'{poly_str}');>;);out;\"\n",
    "print(f\"\\nQuery:\\n{query_str}\")\n",
    "\n",
    "# make the query against overpass\n",
    "resp = requests.get(url, params={\"data\": query_str})\n",
    "assert resp.status_code == 200"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "boxed-navigator",
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "\n",
    "\n",
    "EARTH_CIRCUMFERENCE = 6378137  # earth circumference in meters\n",
    "\n",
    "\n",
    "def great_circle_distance(latlon_a, latlon_b):\n",
    "    \"\"\"Reference calculation via https://gist.github.com/gabesmed/1826175\"\"\"\n",
    "    lat1, lon1 = latlon_a\n",
    "    lat2, lon2 = latlon_b\n",
    "\n",
    "    dLat = math.radians(lat2 - lat1)\n",
    "    dLon = math.radians(lon2 - lon1)\n",
    "    a = (math.sin(dLat / 2) * math.sin(dLat / 2) +\n",
    "            math.cos(math.radians(lat1)) * math.cos(math.radians(lat2)) * \n",
    "            math.sin(dLon / 2) * math.sin(dLon / 2))\n",
    "    c = 2 * math.atan2(math.sqrt(a), math.sqrt(1 - a))\n",
    "    d = EARTH_CIRCUMFERENCE * c\n",
    "    \n",
    "    return d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "german-iraqi",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total distance of cycleway in target area: 35005.33 meters\n"
     ]
    }
   ],
   "source": [
    "# create node id lookup\n",
    "nodes = {n[\"id\"]: [n[\"lat\"], n[\"lon\"]] for n in [e for e in resp.json()[\"elements\"] if e[\"type\"] == \"node\"]}\n",
    "\n",
    "# for each way, get the total length of that cycleway\n",
    "ways = [e for e in resp.json()[\"elements\"] if e[\"type\"] == \"way\"]\n",
    "for way in ways:\n",
    "    way[\"dist_meters\"] = round(sum(\n",
    "        great_circle_distance(nodes[n1], nodes[n2])\n",
    "        for n1, n2 in zip(way[\"nodes\"][:-1], way[\"nodes\"][1:])), 2)\n",
    "    \n",
    "print(f\"Total distance of cycleway in target area: {round(sum(w['dist_meters'] for w in ways), 2)} meters\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "conscious-terrorism",
   "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.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}