500Ft.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Primitive KML/Z Loader \n",
    "\n",
    "## General Info \n",
    "\n",
    "This is based off of `fastkml` \n",
    "\n",
    "## Motivation \n",
    "\n",
    "Existing libraries are pythonic -- they load for use in an object-oriented sense. Instead, I needed something that could load these into simple Numpy/Pandas objects and interact with the data in a rational form. Thus, I'm creating this, to hopefully be able to read through KML and eventually KMZ files and create a comprehensive representation of the data stored. It will, however, by necessity, lose some description statistics. \n",
    "\n",
    "## Goals \n",
    "\n",
    "### KML Parsing \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from fastkml import kml\n",
    "import shapely \n",
    "import pandas as pd \n",
    "import numpy as np \n",
    "import lxml \n",
    "import json \n",
    "import matplotlib.pyplot as plt \n",
    "import scipy as sp\n",
    "from scipy.spatial.distance import pdist, squareform\n",
    "from shapely.ops import cascaded_union\n",
    "\n",
    "#garbage managmeent \n",
    "import gc "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "global counter\n",
    "counter = 0 \n",
    "def get_id_from_description(description_string, mark_config = \"OBJECTID\"): \n",
    "    global counter\n",
    "    if description_string is not None: \n",
    "        string = description_string\n",
    "        string = ' '.join(string.replace(\"\\n\", \"\").split())\n",
    "        string = (string[string.find(\"<tr>\"):])\n",
    "        string = \"<root>\" + string + \"</root>\"\n",
    "        tree = lxml.etree.XML(string)\n",
    "        for row in tree.getchildren(): \n",
    "            f = row.getchildren()[0]\n",
    "            if f.text == mark_config: \n",
    "                return int(row.getchildren()[1].text)\n",
    "    counter += 1\n",
    "    return counter \n",
    "\n",
    "def get_coords(kml_object, resolution=1, mark_config=\"Beat\"): \n",
    "    while (isinstance(list(kml_object.features())[0], kml.Document) \\\n",
    "           or isinstance(list(kml_object.features())[0], kml.Folder)): \n",
    "        kml_object = list(kml_object.features())[0]\n",
    "    frames = {} \n",
    "    for k in kml_object.features(): \n",
    "        geom = shapely.geometry.mapping(k.geometry)\n",
    "        #now convert \n",
    "        arr = geom['coordinates']\n",
    "        if (len(arr) == 1): \n",
    "            arr = [arr] \n",
    "        paths = []\n",
    "        for path in arr:\n",
    "            path = path[0]\n",
    "            nparr = np.array(path, dtype=float) \n",
    "            subpath = []\n",
    "            for i in range(0, nparr.shape[0], resolution): \n",
    "                row = nparr[i, :].tolist() \n",
    "                subpath.append({\"lng\": row[0], \"lat\": row[1]})\n",
    "            paths.append(subpath) \n",
    "        geom['coordinates'] = paths \n",
    "        frames[get_id_from_description(k.description, mark_config)] = geom  \n",
    "    return frames "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(\"data/gtboundary_layer.kml\", 'rt', encoding='utf-8') as myfile: \n",
    "    doc = myfile.read() \n",
    "\n",
    "f = kml.KML() \n",
    "f.from_string(doc) \n",
    "\n",
    "frames = get_coords(f, resolution=1, mark_config=\"FID\") "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "#write to file \n",
    "with open(\"shapes.json\", \"w\") as f:\n",
    "    s = json.dumps(frames)\n",
    "    s = \"var shapes = \" + s; \n",
    "    f.write(s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "now we have to conver this to 500 feet boundary \n",
    "\n",
    "first let's calculate what 1 foot is in terms of longitude and latitude \n",
    "\n",
    "This doesn't have to be too exact, so some error is fine \n",
    "\n",
    "Thus, we'll define our vectors as such: \n",
    "\n",
    "POSTIVE X is MORE POSITIVE LONGITUDE \n",
    "\n",
    "POSITIVE Y is MORE POSITIVE LATITUDE \n",
    "\n",
    "We'll be using a reverse haversine implementation \n",
    "\n",
    "The Haversine Formula Says: \n",
    "$$ hav \\left(\\frac d r\\right) = hav(\\phi_2 - \\phi_1) + \\cos(\\phi_1)\\cos(\\phi_2) hav(\\lambda_2 - \\lambda_1)$$\n",
    "where\n",
    "$$ \n",
    "hav(\\theta) = \\sin^2 \\left( \\frac \\theta 2 \\right) = \\frac {1 - cos(\\theta) } {2} $$\n",
    "and $\\phi$ is latitude, $\\lambda$ is longitude, both in radians "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def hav(theta): \n",
    "    return (1 - np.cos(theta)) / 2\n",
    "\n",
    "def haversine(lon1, lat1, lon2, lat2): \n",
    "    #we expect radians as in put \n",
    "    return hav(lat2 - lat1) + np.cos(lat1) * np.cos(lat2) * hav(lon2 - lon1) \n",
    "\n",
    "def calc_distance(**k): \n",
    "    havs = haversine(k[\"lon1\"], k[\"lat1\"], k[\"lon2\"], k[\"lat2\"])\n",
    "    c =  2 * np.arctan2(np.sqrt(havs), np.sqrt(1 - havs)) \n",
    "    invs = 6371e3 * c  \n",
    "    return invs\n",
    "def meters_to_feet(x): \n",
    "    return 3.28084 * x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we'll find the average longitude and latitude of our dataset. Then, we'll compute the distance between longitude, holding lat constant, and vice versa. The inverse will give us a lon-lat vector for unit distance. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "#find vals \n",
    "lats = [] \n",
    "lngs = [] \n",
    "for k in frames:\n",
    "    for b in frames[k]['coordinates']: \n",
    "        for a in b: \n",
    "            lats.append(a['lat'])\n",
    "            lngs.append(a['lng'])      \n",
    "\n",
    "#compute means \n",
    "mean_lat = np.deg2rad(np.mean(lats))\n",
    "mean_lon = np.deg2rad(np.mean(lngs))\n",
    "\n",
    "#and max/min \n",
    "min_lat = np.deg2rad(min(lats)) \n",
    "max_lat = np.deg2rad(max(lats))\n",
    "\n",
    "min_lon = np.deg2rad(min(lngs)) \n",
    "max_lon = np.deg2rad(max(lngs)) \n",
    "\n",
    "#now find the distances \n",
    "movement_along_lat = meters_to_feet(calc_distance(lon1 = min_lon, \n",
    "                                                 lat1 = mean_lat, \n",
    "                                                 lon2 = max_lon, \n",
    "                                                 lat2 = mean_lat))\n",
    "movement_along_lng = meters_to_feet(calc_distance(lon1 = mean_lon, \n",
    "                                                 lat1 = min_lat, \n",
    "                                                 lon2 = mean_lon, \n",
    "                                                 lat2 = max_lat))\n",
    "LonPer500Ft = np.rad2deg((max_lon - min_lon) / (movement_along_lat / 500))\n",
    "LatPer500Ft = np.rad2deg((max_lat - min_lat) / (movement_along_lng / 500))\n",
    "\n",
    "del lats \n",
    "del lngs \n",
    "gc.collect() \n",
    "coordVect =[ LonPer500Ft, LatPer500Ft]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we need to interpolate to fill in any gaps. We'll use a resolution of 250 as half the required space. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 / 34\n",
      "2 / 34\n",
      "3 / 34\n",
      "4 / 34\n",
      "5 / 34\n",
      "6 / 34\n",
      "7 / 34\n",
      "8 / 34\n",
      "9 / 34\n",
      "10 / 34\n",
      "11 / 34\n",
      "12 / 34\n",
      "13 / 34\n",
      "14 / 34\n",
      "15 / 34\n",
      "16 / 34\n",
      "17 / 34\n",
      "18 / 34\n",
      "19 / 34\n",
      "20 / 34\n",
      "21 / 34\n",
      "22 / 34\n",
      "23 / 34\n",
      "24 / 34\n",
      "25 / 34\n",
      "26 / 34\n",
      "27 / 34\n",
      "28 / 34\n",
      "29 / 34\n",
      "30 / 34\n",
      "31 / 34\n",
      "32 / 34\n",
      "33 / 34\n",
      "34 / 34\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Now, linear interpolation to increase resolution: \n",
    "def get_distance_in_feet(p1, p2): \n",
    "    return meters_to_feet(calc_distance(lon1 =  np.deg2rad(p1['lng']), \n",
    "                                       lon2 = np.deg2rad(p2['lng']), \n",
    "                                       lat1 = np.deg2rad(p1['lat']), \n",
    "                                       lat2 = np.deg2rad(p2['lat'])))\n",
    "def interpolate(p1, p2): \n",
    "    return {'lat': (p1['lat'] + p2['lat'])/2, \n",
    "            'lng': (p1['lng'] + p2['lng'])/2}\n",
    "count = 0 \n",
    "RESOLUTION = 250\n",
    "for f in frames: \n",
    "    for vi, v in enumerate(frames[f]['coordinates']): \n",
    "        k = v \n",
    "        finishedInterpolating = False\n",
    "        inter_count = 0 \n",
    "        while not finishedInterpolating:\n",
    "            finishedInterpolating = True \n",
    "            interpolated = [] \n",
    "            for i in range(0, len(k) - 1): \n",
    "                interp = None\n",
    "                if get_distance_in_feet(k[i], k[i + 1]) > RESOLUTION: \n",
    "                    #average coord \n",
    "                    interp = interpolate(k[i], k[i + 1])\n",
    "                interpolated.append(k[i])\n",
    "                if interp is not None: \n",
    "                    finishedInterpolating = False\n",
    "                    interpolated.append(interp) \n",
    "                interpolated.append(k[i + 1]) \n",
    "            k = interpolated \n",
    "            inter_count += 1\n",
    "        frames[f]['coordinates'][vi] = k \n",
    "        count += 1\n",
    "        print(count, \"/\", 34) \n",
    "del interpolated \n",
    "del k \n",
    "del v \n",
    "gc.collect() "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Some functions to use later \n",
    "def zip_coordinates(coords): \n",
    "    xs = list(map(lambda x: x['lng'], coords))\n",
    "    ys = list(map(lambda y: y['lat'], coords))\n",
    "    return (xs, ys) \n",
    "\n",
    "def plot_dict_coors_on_ax(coords, ax, shape=False): \n",
    "    xs, ys = zip_coordinates(coords)\n",
    "    if not shape: \n",
    "        ax.scatter(xs, ys) \n",
    "    else: \n",
    "        ax.plot(xs, ys)\n",
    "    return ax \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plt.close('all')\n",
    "import gc\n",
    "gc.collect()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1\n",
      "2\n",
      "3\n",
      "4\n",
      "5\n",
      "6\n",
      "7\n",
      "8\n",
      "9\n",
      "10\n",
      "11\n",
      "12\n",
      "13\n",
      "14\n",
      "15\n",
      "16\n",
      "17\n",
      "18\n",
      "19\n",
      "20\n",
      "21\n",
      "22\n",
      "23\n",
      "24\n",
      "25\n",
      "26\n",
      "27\n",
      "28\n",
      "29\n",
      "30\n",
      "31\n",
      "32\n",
      "33\n",
      "34\n"
     ]
    }
   ],
   "source": [
    "#Let's go through and regenerate polygons for everything again \n",
    "#Ellipse Template \n",
    "circle = np.array([i for i in np.linspace(0, 2 * np.pi, 50)]).T \n",
    "circle = circle.reshape(circle.shape[0], 1)\n",
    "circle = np.hstack((np.cos(circle), np.sin(circle))) \n",
    "circle = np.multiply(coordVect, circle) \n",
    "final_shapes = []\n",
    "count = 0 \n",
    "for k in frames: \n",
    "    for j in frames[k]['coordinates']: \n",
    "        coordinates = zip_coordinates(j) \n",
    "        #zip em up \n",
    "        zip_coor = list(zip(*coordinates)) \n",
    "        #And now create our shape \n",
    "        ogShape = shapely.geometry.polygon.Polygon(zip_coor)\n",
    "        #And create our circles \n",
    "        circles = [ogShape] \n",
    "        for point in zip_coor: \n",
    "            #generate circle \n",
    "            circles.append(shapely.geometry.polygon.Polygon(np.add(point, circle))) \n",
    "        #Create a Union \n",
    "        union_shape = cascaded_union(circles)\n",
    "        \n",
    "        #Again, extract the lat-lon coordinates from this. \n",
    "        final_shapes.append(union_shape) \n",
    "        count += 1\n",
    "        print(count) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "final_shape = cascaded_union(final_shapes)\n",
    "#convert to our JSON version \n",
    "points_array = [] \n",
    "for k in zip(*final_shape.exterior.xy):\n",
    "    points_array.append({'lng': k[0], 'lat': k[1]})\n",
    "\n",
    "final_frame = {0: {'coordinates': [points_array]}}\n",
    "with open(\"f_shapes.json\", \"w\") as f:\n",
    "    s = json.dumps(final_frame)\n",
    "    s = \"var o_shapes = \" + s; \n",
    "    f.write(s)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:data]",
   "language": "python",
   "name": "conda-env-data-py"
  },
  "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}