working_with_geospatial_data.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Working With Geospatial Data in Python\n",
    "\n",
    "### Chad Hawkins"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# What Are Geospatial Data?\n",
    "\n",
    "* Simply put, data that have a geographic component\n",
    "\n",
    "# What is GIS?\n",
    "\n",
    "* Geographic Information Systems\n",
    "* Systems for storing, displaying and interacting with geospatial data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Types of Geospatial Data\n",
    "\n",
    "* Vector\n",
    "* Raster\n",
    "* Tabular"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Vector Data\n",
    "\n",
    "* Common formats include GeoJSON, Shapefile, and KML\n",
    "* Point, Line, Polygon\n",
    "* Could be derived from GPS or digitzed from imagery\n",
    "* Uses include observations, parcels, roads, soils, administrative boundaries"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "# GeoJSON Example\n",
    "\n",
    "<br>\n",
    "```\n",
    "{\n",
    "  \"type\": \"FeatureCollection\",\n",
    "  \"features\": [\n",
    "    {\n",
    "      \"type\": \"Feature\",\n",
    "      \"geometry\": {\n",
    "        \"type\": \"Point\",\n",
    "        \"coordinates\": [-77.410133, 39.413948]\n",
    "      },\n",
    "      \"properties\": {\n",
    "        \"name\": \"IronNet Office\"\n",
    "      }\n",
    "    }\n",
    "  ]\n",
    "}\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "# KML Example\n",
    "\n",
    "<br>\n",
    "```\n",
    "<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n",
    "<kml xmlns=\"http://www.opengis.net/kml/2.2\">\n",
    "\n",
    "  <Placemark>\n",
    "\n",
    "    <name>A simple placemark on the ground</name>\n",
    "\n",
    "    <Point>\n",
    "\t\t\t<coordinates>8.542952335953721,47.36685263064198,0</coordinates>\n",
    "    </Point>\n",
    "\n",
    "  </Placemark>\n",
    "\n",
    "</kml>\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "# Soils\n",
    "\n",
    "![alt text](images/soils.png \"Soils\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Raster Data\n",
    "\n",
    "* Common formats include GeoTIFF, Erdas Imagine, and GRIB\n",
    "* Used primarily for remotely sensed data (satellite imagery, weather, elevation)\n",
    "* Continuous or categorical data\n",
    "* Efficient for both rendering and storing but ultimately everything is a rectangle"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "# Landcover Classification (Categorical)\n",
    "\n",
    "![alt text](images/landcover_classification.png \"Landcover Classification Raster\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "# Elevation (Continuous)\n",
    "\n",
    "![alt text](images/elevation.png \"Elevation Raster\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Tabular Data\n",
    "\n",
    "* Data such as addresses and zip codes\n",
    "* Not geospatial in the strictest terms but correspond to well understood geographic regions\n",
    "* Addresses are translated to real world coordinates using a process called geocoding"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Map Projections\n",
    "\n",
    "* Transforming locations on a spheroid to coordinates on a plane\n",
    "* All map projections distort some of shape, area, distance, and direction\n",
    "* We must choose which properties are critical to our analysis\n",
    "* Spatial operations must occur between data in the same projection\n",
    "* In this presentation with are working with unprojected data (latitude and longitude)\n",
    "* As such, we cannot compute meaningful areas and distances must be approximated (Haversine formula)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Spatial Operations\n",
    "\n",
    "* Intersect\n",
    "* Intersection/Clip\n",
    "* Within/Contains\n",
    "* Buffer\n",
    "* Merge\n",
    "* Zonal Statistics (raster statistics for vector-defined region)\n",
    "* Vectorize\n",
    "* Rasterize"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Python Packages\n",
    "\n",
    "### Data Processing\n",
    "\n",
    "* Shapely\n",
    "* Fiona\n",
    "* Rasterio\n",
    "* Pyproj\n",
    "* GeoPandas\n",
    "\n",
    "### Web Application\n",
    "* GeoDjango\n",
    "\n",
    "### Mapping\n",
    "* Basemap\n",
    "* Folium"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Related Tools\n",
    "\n",
    "### Integrated Tools\n",
    "* QGIS\n",
    "* ArcGIS\n",
    "\n",
    "### Display\n",
    "* Leaflet\n",
    "* OpenLayers\n",
    "* Google Earth\n",
    "* Mapnik\n",
    "\n",
    "### Processing and Storage\n",
    "* GDAL/OGR\n",
    "* PostgreSQL/PostGIS\n",
    "* R"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Data Sources\n",
    "\n",
    "* [Maryland GIS Data Catalog](https://data.imap.maryland.gov/)\n",
    "* [Frederick County](https://frederickcountymd.gov/5969/Download-GIS-Data)\n",
    "* [ArcGIS Open Data](https://hub.arcgis.com/pages/open-data)\n",
    "* [OpenStreetMap](https://www.openstreetmap.org)\n",
    "* [Census Bureau TIGER/Line](https://www.census.gov/geographies/mapping-files/time-series/geo/tiger-line-file.html)\n",
    "* [USGS](https://earthexplorer.usgs.gov/)\n",
    "* [NWS](https://www.weather.gov/gis/)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Let's Do Some Spatial Analysis\n",
    "\n",
    "## ... of Bigfoot sightings!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# First We Need to Acquire and Prepare Some Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "from collections import defaultdict\n",
    "from copy import deepcopy\n",
    "import datetime\n",
    "import fiona\n",
    "import folium\n",
    "from folium.plugins import MarkerCluster\n",
    "import haversine\n",
    "import pandas as pd\n",
    "from pprint import pprint\n",
    "import requests\n",
    "from shapely.geometry import mapping, Point, shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "url = 'https://opendata.arcgis.com/datasets/d0afc5b29e4346cc9a4cf8e43bcaaed0_0.geojson'\n",
    "\n",
    "response = requests.get(url)\n",
    "data = response.json()\n",
    "\n",
    "print(data['type'])\n",
    "\n",
    "pprint(data['features'][0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# Let's convert TimeWhen to a datetime\n",
    "timewhen_format = '%Y-%m-%d %H:%M:%S'\n",
    "\n",
    "valid_sightings = []\n",
    "sighting_times = []\n",
    "for feature in deepcopy(data['features']):\n",
    "    # Ignore sightings with bad or missing time observations\n",
    "    try:\n",
    "        sighting_time = datetime.datetime.strptime(\n",
    "            feature['properties']['TimeWhen'], \n",
    "            timewhen_format\n",
    "        )\n",
    "    except (TypeError, ValueError):\n",
    "        continue\n",
    "    feature['properties']['TimeWhen'] = sighting_time\n",
    "    valid_sightings.append(feature)\n",
    "    sighting_times.append(sighting_time)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "print('Total sightings: {}'.format(len(data['features'])))\n",
    "print('Valid sightings: {}'.format(len(valid_sightings)))\n",
    "print('Oldest sighting: {}'.format(min(sighting_times)))\n",
    "print('Most recent sighting: {}'.format(max(sighting_times)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# Let's restrict our map to more recent sightings\n",
    "recent_sightings = [\n",
    "    sighting for sighting in valid_sightings if \n",
    "    sighting['properties']['TimeWhen'] >= datetime.datetime(2010, 1, 1)\n",
    "]\n",
    "\n",
    "print('Recent sightings: {}'.format(len(recent_sightings)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Let's Look at These Data on a Map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# Let's start with an empty map\n",
    "# Create a function so we can reuse later\n",
    "def get_map_object():\n",
    "    return folium.Map(\n",
    "        location=[44, -103],\n",
    "        tiles='Stamen Terrain',\n",
    "        zoom_start=4\n",
    "    )\n",
    "\n",
    "m = get_map_object()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "m = get_map_object()\n",
    "\n",
    "for feature in recent_sightings:\n",
    "    lon, lat = feature['geometry']['coordinates']\n",
    "    folium.Marker(\n",
    "        [lat, lon],\n",
    "    ).add_to(m)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# That's a Bit Overwhelming..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# Recreate the map using clusters\n",
    "\n",
    "m = get_map_object()\n",
    "\n",
    "locations = []\n",
    "for feature in recent_sightings:\n",
    "    lon, lat = feature['geometry']['coordinates']\n",
    "    locations.append((lat, lon))\n",
    "cluster = MarkerCluster(locations=locations)\n",
    "\n",
    "m.add_child(cluster)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Much Better... Let's Customize a Bit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "m = get_map_object()\n",
    "\n",
    "locations = []\n",
    "popups = []\n",
    "icons = []\n",
    "for feature in recent_sightings:\n",
    "    lon, lat = feature['geometry']['coordinates']\n",
    "    locations.append((lat, lon))\n",
    "    # Show the Description field when clicked\n",
    "    popup = folium.Popup(feature['properties']['Descr'], max_width=200)\n",
    "    popups.append(popup)\n",
    "    # Use a custom icon\n",
    "    icon = folium.map.Icon(color='gray', icon='paw', prefix='fa')\n",
    "    icons.append(icon)\n",
    "\n",
    "cluster = MarkerCluster(\n",
    "    locations=locations, \n",
    "    popups=popups, \n",
    "    icons=icons, \n",
    "    name='Bigfoot Sightings'\n",
    ")\n",
    "\n",
    "m.add_child(cluster)\n",
    "folium.LayerControl().add_to(m)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Let's Count the Sightings by State\n",
    "\n",
    "* Points on a map can be difficult to properly interpret\n",
    "* It's often useful to aggregate in some fashion\n",
    "* In this case, by state seems reasonable\n",
    "* First we need to associate the point data with a state"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "states = []\n",
    "\n",
    "# Open the shapefile with fiona\n",
    "path = 'tl_2018_us_state/tl_2018_us_state.shp'\n",
    "with fiona.open(path) as shp:\n",
    "    pprint(shp.meta)\n",
    "    for record in shp:\n",
    "        # Store the shapely object for the geometry so we can do our comparison\n",
    "        record['geometry'] = shape(record['geometry'])\n",
    "        states.append(record)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "sightings_by_state = defaultdict(int)\n",
    "\n",
    "for sighting in recent_sightings:\n",
    "    location = shape(sighting['geometry'])\n",
    "    # Loop over the states until we find one that contains our sighting\n",
    "    # This is VERY inefficient, especially with larger datasets\n",
    "    # The proper way to do this is to use an rtree index \n",
    "    # (beyond the scope of this presentation)\n",
    "    for state in states:\n",
    "        if state['geometry'].contains(location):\n",
    "            state_abbr = state['properties']['STUSPS']\n",
    "            sightings_by_state[state_abbr] += 1\n",
    "            break\n",
    "\n",
    "# Create a dataframe\n",
    "df = pd.DataFrame(\n",
    "    [(state_abbr, count) for state_abbr, count in sightings_by_state.items()],\n",
    "    columns=['State', 'Count']\n",
    ")\n",
    "\n",
    "df.sort_values(by='Count', ascending=False).head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# Create GeoJSON\n",
    "for state in states:\n",
    "    # We simplify the geometry here because our boundaries \n",
    "    # are far more detailed than required\n",
    "    geometry = state['geometry'].simplify(0.003, preserve_topology=True)\n",
    "    state['geometry'] = mapping(geometry)\n",
    "    \n",
    "states_geojson = {\n",
    "    'type': 'FeatureCollection', \n",
    "    'features': states\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Let's Create a Choropleth Map\n",
    "\n",
    "* \"A thematic map in which areas are shaded or patterned in proportion to the measurement of the statistical variable being displayed on the map\" [1]\n",
    "* In this case the variable is the number of sightings for each state\n",
    "\n",
    "[1]: https://en.wikipedia.org/wiki/Choropleth_map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "m = get_map_object()\n",
    "\n",
    "folium.Choropleth(\n",
    "    geo_data=states_geojson,\n",
    "    name='Sightings by State',\n",
    "    data=df,\n",
    "    columns=['State', 'Count'],\n",
    "    key_on='feature.properties.STUSPS',\n",
    "    fill_color='YlOrRd',\n",
    "    legend_name='Number of Sightings',\n",
    "    nan_fill_opacity=0.9,\n",
    "    nan_fill_color='white',\n",
    "    fill_opacity=0.9\n",
    ").add_to(m)\n",
    "\n",
    "folium.LayerControl().add_to(m)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Let's Find the Closest Sighting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "ironnet_location = (39.413948, -77.410133)\n",
    "\n",
    "closest_distance = None\n",
    "closest_sighting = None\n",
    "\n",
    "for sighting in valid_sightings:\n",
    "    lon, lat = sighting['geometry']['coordinates']\n",
    "    sighting_location = (lat, lon)\n",
    "    # Get the distance in miles\n",
    "    distance = haversine.haversine(\n",
    "        ironnet_location, \n",
    "        sighting_location, \n",
    "        unit='mi'\n",
    "    )\n",
    "    if closest_distance is None:\n",
    "        closest_distance = distance\n",
    "        closest_sighting = sighting\n",
    "    elif distance < closest_distance:\n",
    "        closest_distance = distance\n",
    "        closest_sighting = sighting\n",
    "\n",
    "desc = closest_sighting['properties']['Descr']\n",
    "print(\n",
    "    'Closest sighting was {} miles away on {}\\n{}'.format(\n",
    "        round(closest_distance, 2), \n",
    "        closest_sighting['properties']['TimeWhen'].date().isoformat(),\n",
    "        desc[desc.index('Report'):desc.index('</b>')]\n",
    "    )\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Questions?\n",
    "\n",
    "* Email: cwh@chadwhawkins.com\n",
    "* GitHub: @chadwhawkins\n",
    "* Frederick Tech Slack"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "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.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}