Geodesic Distances in Python.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import json, vincenty"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## We live on a sphere, and that makes this hard\n",
    "\n",
    "Given a bunch of points, compute some distances with shapely and pyproj\n",
    "\n",
    "For simplicity, we keep all of our points in WGS84, we can use vincenty distance for super accurate results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Let's load a user and get their points in WGS84"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "points = []\n",
    "with open('/data/chime/geo2/NYC/ZoneA/alex_tonk.geojsonl') as inFile:\n",
    "    for idx, line in enumerate(inFile):\n",
    "        t = json.loads(line)\n",
    "        points.append(t['geometry']['coordinates'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "points is now an array of points in lon, lat (x,y) format"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[-73.96226984, 40.78033596], [-73.96387997, 40.77846793]]"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "points[:2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use Vincenty Distance.\n",
    "\n",
    "NOTE: To use Vincenty distance, the format is: (lat1, lon1, lat2, lon2), so reverse the GeoJSON representation of (x,y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Distance between [-73.96226984, 40.78033596] and [-73.96226984, 40.78033596]: 0.0 kilometers\n",
      "Distance between [-73.96226984, 40.78033596] and [-73.96387997, 40.77846793]: 0.248006 kilometers\n",
      "Distance between [-73.96226984, 40.78033596] and [-73.81858043, 40.75890912]: 12.362479 kilometers\n",
      "Distance between [-73.96387997, 40.77846793] and [-73.96226984, 40.78033596]: 0.248006 kilometers\n",
      "Distance between [-73.96387997, 40.77846793] and [-73.96387997, 40.77846793]: 0.0 kilometers\n",
      "Distance between [-73.96387997, 40.77846793] and [-73.81858043, 40.75890912]: 12.458236 kilometers\n",
      "Distance between [-73.81858043, 40.75890912] and [-73.96226984, 40.78033596]: 12.362479 kilometers\n",
      "Distance between [-73.81858043, 40.75890912] and [-73.96387997, 40.77846793]: 12.458236 kilometers\n",
      "Distance between [-73.81858043, 40.75890912] and [-73.81858043, 40.75890912]: 0.0 kilometers\n"
     ]
    }
   ],
   "source": [
    "for p1 in points[:3]:\n",
    "    for p2 in points[:3]:\n",
    "        distance = vincenty.vincenty(list(reversed(p1)), list(reversed(p2)))\n",
    "        print(\"Distance between {0} and {1}: {2} kilometers\".format(p1,p2,distance))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following link confirms that Google agrees with this distance: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "http://townsendjennings.com/geo/?geojson={\"type\":\"FeatureCollection\",\"features\":[{\"geometry\":{\"coordinates\":[-73.96226984,40.78033596],\"type\":\"Point\"},\"properties\":{},\"type\":\"Feature\"},{\"geometry\":{\"coordinates\":[-73.96387997,40.77846793],\"type\":\"Point\"},\"properties\":{},\"type\":\"Feature\"},{\"geometry\":{\"coordinates\":[-73.81858043,40.75890912],\"type\":\"Point\"},\"properties\":{},\"type\":\"Feature\"}]}\n"
     ]
    }
   ],
   "source": [
    "feats = []\n",
    "for p in points[:3]: \n",
    "    feats.append({\"type\":\"Feature\",\"properties\":{},\"geometry\":{\"type\":\"Point\", \"coordinates\":[p[0],p[1]]}})\n",
    "print(\"http://townsendjennings.com/geo/?geojson=\"+json.dumps({\"type\":\"FeatureCollection\",\"features\":feats}).replace(\" \",\"\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Geometric Centorid?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "27"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(points)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-73.95001657629629, 40.73189305185184)"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum(p[0] for p in points)/len(points), sum(p[1] for p in points)/len(points)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from shapely import geometry"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "pts = geometry.MultiPoint([geometry.Point(x) for x in points])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "POINT (-73.95001657629629 40.73189305185184)\n"
     ]
    }
   ],
   "source": [
    "print(pts.centroid)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now add some units?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Eh... that's probably ok'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "IPython (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.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}