notebook_pyconby.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "![Title](img/pycon.png)\n",
    "\n",
    "\n",
    "## Geospatial Data Processing, Analysis and Visualization using Python\n",
    "\n",
    "Prof. Martin Christen<br/>\n",
    "mailto:martin.christen@fhnw.ch<br/>\n",
    "Twitter: @MartinChristen<br/>\n",
    "LinkedIn: https://www.linkedin.com/in/martinchristen/\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### This notebook can be found at: https://github.com/martinchristen/PyConBy2019\n",
    "\n",
    "\n",
    "\n",
    "In this notebook the following will be covered:\n",
    "\n",
    "**Vector File Formats** (Shapefiles, GeoJSON, KML, GeoPackage, ...)<br/>\n",
    "**Raster File Formats** (GeoTiff, JPG, ...)\n",
    "\n",
    "* Shapely: Manipulation and analysis of geometric objects ( https://github.com/Toblerity/Shapely )\n",
    "* Fiona - The pythonic way to handle vector data ( https://github.com/Toblerity/Fiona )\n",
    "* rasterio - The pythonic way to handle raster data ( https://github.com/mapbox/rasterio )\n",
    "* pyproj - transforming spatial reference systems ( https://github.com/jswhit/pyproj ) \n",
    "* Geospatial analysis using GeoPandas ( https://github.com/geopandas/geopandas ) \n",
    "* Creating maps using Folium ( https://github.com/python-visualization/folium ) and/or ipyleaflet ( https://github.com/jupyter-widgets/ipyleaflet )\n",
    "  \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "\n",
    "### Data used in this notebook\n",
    "\n",
    "* Natural Earth Dataset, https://www.naturalearthdata.com/ (public domain)\n",
    "* Blue Marble: Next Generation was produced by Reto Stöckli, NASA Earth Observatory (NASA Goddard Space Flight Center\n",
    "* geonames: http://www.geonames.org/ CC License\n",
    "                                                                                     \n",
    "To keep it simple, I put it on the git repository"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "![Title](img/overview_tools.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## 1. Installation\n",
    "\n",
    "(This tutorial requires anaconda, if you don't have it yet, download it here: https://www.anaconda.com/download/ )\n",
    "\n",
    "**This notebook requires Python 3.6 (or higher)** (jupyter Notebook or Jupyter Lab using Chrome)\n",
    "\n",
    "### Installing Modules (conda)\n",
    "\n",
    "Install the main modules using conda, dependencies will be resolved (gdal etc.)\n",
    "\n",
    "    conda install shapely\n",
    "    conda install fiona\n",
    "    conda install rasterio\n",
    "    conda install geopandas\n",
    "    conda install folium -c conda-forge\n",
    "    conda install ipyleaflet -c conda-forge "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## 2. Shapely\n",
    "\n",
    "Shapely is a BSD-licensed Python package for **manipulation** and **analysis** of **planar geometric objects**. \n",
    "\n",
    "* Shapely is **not** concerned with data formats or coordinate systems.\n",
    "* Shapely is based on the widely deployed GEOS (the engine of PostGIS) and JTS (from which GEOS is ported) libraries.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Simple Feature Access\n",
    "\n",
    "**Simple Feature Access** is both an Open Geospatial Consortium (OGC) and International Organization for Standardization (ISO) standard **ISO 19125** that specifies a common storage and access model of mostly two-dimensional geometries (point, line, polygon, multi-point, multi-line, etc.) used by geographic information systems.\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Shapely supports the following **Simple Feature Access** Features:\n",
    "\n",
    "* Point\n",
    "* LineString\n",
    "* LinearRing        \n",
    "* Polygon\n",
    "* MultiLineString\n",
    "* MultiPoint\n",
    "* MultiPolygon\n",
    "\n",
    "![Features](img/Features.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Examples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "from shapely.geometry import Polygon\n",
    "\n",
    "polygon1 = Polygon([(30, 10), (40, 40), (20, 35), (10, 20), (30, 10)])\n",
    "\n",
    "print(f\"Polygon area: {polygon1.area}, polygon length: {polygon1.length}\") "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "You can display shapely objects directly in the Jupyter Notebook (usually for \"debug\" purposes)\n",
    "\n",
    "It is also possible to display this in matplotlib/descartes, we will see later how to bring this directly on a map)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "polygon1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "polygon2 = Polygon([(20,20),(80,30),(50,40),(20,20)])\n",
    "polygon2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "polygon2.union(polygon1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "polygon2.intersection(polygon1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "polygon2.symmetric_difference(polygon1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "result = polygon2.symmetric_difference(polygon1)\n",
    "\n",
    "print(f\"Polygon area: {result.area}, polygon length: {result.length}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Well Known Text (WKT)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "result.wkt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "s = result.wkt\n",
    "type(s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "import shapely.wkt\n",
    "\n",
    "mypolygon = shapely.wkt.loads(s)\n",
    "mypolygon"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Binary operations on shapes:\n",
    "\n",
    "- **contains** (Returns True if the interior of the object intersects the interior of the other but does not contain it, and the dimension of the intersection is less than the dimension of the one or the other.)\n",
    "- **intersects** (Returns True if the boundary and interior of the object intersect in any way with those of the other.)\n",
    "- **witin** (Returns True if the object’s boundary and interior intersect only with the interior of the other (not its boundary or exterior).\n",
    "- **touches** (Returns True if the objects have at least one point in common and their interiors do not intersect with any part of the other.)\n",
    "- **crosses** (Returns True if the interior of the object intersects the interior of the other but does not contain it, and the dimension of the intersection is less than the dimension of the one or the other.)\n",
    "- **equals** (Returns True if the set-theoretic boundary, interior, and exterior of the object coincide with those of the other.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "polygon1.intersects(polygon2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "polygon1.within(polygon2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "polygon1.equals(polygon1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## 3. Fiona\n",
    "\n",
    "Spatial information is not only geometry. Together with it we always have metadata (properties). Lets look at a first dataset.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Reading Vector Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "import fiona\n",
    "\n",
    "c = fiona.open('data/ne_10m_airports/ne_10m_airports.shp', 'r')\n",
    "\n",
    "airport = next(iter(c))\n",
    "airport"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "airport['properties']['name']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "airport['geometry']['type']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "airport['geometry']['coordinates']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "c.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Two ways to read vector data\n",
    "\n",
    "Lets read all data. There are basically two ways:\n",
    "\n",
    "a) load everything to memory: (if dataset isn't too large...)\n",
    "\n",
    "    alldata = list(c)\n",
    "    \n",
    "b) iterate through all data: (one by one):\n",
    "\n",
    "    for element in c:\n",
    "        ...\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "with fiona.open('data/ne_10m_airports/ne_10m_airports.shp', 'r') as c:\n",
    "    for airport in c:\n",
    "        if airport['properties']['iata_code'] == \"MSQ\":\n",
    "            print(airport['properties']['name'])\n",
    "            print(airport['geometry']['coordinates'])\n",
    "            print(airport['properties']['wikipedia'])        \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Coordinate System ? (Coordinate Reference System?)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "with fiona.open('data/ne_10m_airports/ne_10m_airports.shp', 'r') as c:\n",
    "    print(c.crs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "EPSG: 4326 look it up at: http://epsg.io/\n",
    "\n",
    "( http://epsg.io/4326 )\n",
    "\n",
    "\n",
    "For transformations to other spatial reference systems, use **pyproj**\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Lets look at another dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "import fiona\n",
    "\n",
    "c = fiona.open('data/ne_10m_admin_0_countries/ne_10m_admin_0_countries.shp', 'r', encoding=\"utf-8\")\n",
    "\n",
    "country = next(iter(c))\n",
    "#country"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "print(country['properties']['NAME'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "print(country['properties']['NAME_ZH'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "print(country['properties']['NAME_RU'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "print(country['properties']['CONTINENT'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "print(country['properties']['POP_EST'])\n",
    "print(country['properties']['POP_YEAR'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Example: Find the \"Belarus\" Polygon and store it"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "with fiona.open('data/ne_10m_admin_0_countries/ne_10m_admin_0_countries.shp', 'r') as c:\n",
    "    for country in c:\n",
    "        if country['properties']['NAME'] == \"Belarus\":\n",
    "            print(country['properties']['POP_EST'])\n",
    "            print(country['properties']['POP_YEAR'])   \n",
    "            print(country['geometry']['type'])\n",
    "            # print(country['geometry']['coordinates']) # you don't want to print all coordinates!\n",
    "            geometry = country['geometry']\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Example: find Airports inside the Belarus polygon"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "from shapely.geometry import shape\n",
    "from shapely.geometry import MultiPolygon\n",
    "from shapely.geometry import Polygon\n",
    "\n",
    "by = Polygon(shape(geometry))\n",
    "by"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "from shapely.geometry import Point\n",
    "\n",
    "with fiona.open('data/ne_10m_airports/ne_10m_airports.shp', 'r') as c:\n",
    "    for airport in c:      \n",
    "        position = Point(airport['geometry']['coordinates'])\n",
    "        if position.within(by):\n",
    "            print(airport['properties']['iata_code'], airport['properties']['name'], airport['geometry']['coordinates'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## 4. Rasterio"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "### Reading data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "import rasterio\n",
    "\n",
    "dataset = rasterio.open('data/BlueMarble.tif', 'r')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "dataset.name"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "dataset.mode"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "dataset.count # number of raster bands, in our case 3 for r,g,b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "dataset.indexes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "dataset.width, dataset.height"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "dataset.crs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "dataset.bounds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "dataset.transform  # affine transformation pixel to crs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "dataset.transform * (0, 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "~dataset.transform # inverse affine transformation crs to pixel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "~dataset.transform * (0,0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "px,py = ~dataset.transform * (27.487117, 53.869146) # Minsk coordinates to pixel coordinates! \n",
    "print(px,py)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "print(dataset.bounds.left)\n",
    "print(dataset.bounds.bottom)\n",
    "print(dataset.bounds.right)\n",
    "print(dataset.bounds.top)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "r = dataset.read(1)\n",
    "g = dataset.read(2)\n",
    "b = dataset.read(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "rgb = np.dstack((r,g,b))  # stack r,g,b so we can display it...\n",
    "rgb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(15,9))\n",
    "ax.imshow(rgb, interpolation='nearest')\n",
    "ax.plot(px,py, 'ro'); # IBB Minsk"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## 5. Geopandas\n",
    "\n",
    "GeoPandas is an open source project to make working with geospatial data in python easier. GeoPandas **extends** the datatypes used by pandas to **allow spatial operations** on **geometric types**. Geometric operations are performed by shapely. Geopandas further depends on fiona for file access and descartes and matplotlib for plotting.\n",
    "(geopandas.org)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "\n",
    "df = pd.read_csv('data/cities5k.csv', encoding=\"utf-8\", sep=\",\", header=None, low_memory=False)\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "df2 = df[[1,4,5,14]]\n",
    "df2.columns = [\"name\", \"lat\", \"lng\", \"population\"]\n",
    "df2.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "df2.query(\"name == 'Minsk'\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Geometry Column\n",
    "\n",
    "Now lets create a geopandas dataframe. Basically just create a new column \"geometry\" - using Shapely to create Points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "import geopandas as gpd\n",
    "from shapely.geometry import Point\n",
    "\n",
    "geometry = [Point(pos) for pos in zip(df2['lng'], df2['lat'])]\n",
    "gdf = gpd.GeoDataFrame(df2, geometry=geometry)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "gdf.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "gdf = gdf.drop(['lat', 'lng'], axis=1) # remove redundant data, we only want \"geometry\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "gdf.plot(color='green', markersize=1, figsize=(15,9))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# Export to Shapefile\n",
    "# gdf.to_file(\"cities.shp\", driver=\"Shapefile\", encoding=\"utf-8\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "cities = gdf\n",
    "cities.head()  # store as \"cities\" for later"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# another query:\n",
    "\n",
    "bigcities = cities[cities.population > 10_000_000]\n",
    "\n",
    "bigcities = bigcities.sort_values(['population'], ascending=False)\n",
    "bigcities"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Dataset with Polygons\n",
    "\n",
    "Now we load the admin 0 country files again:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "countries = gpd.read_file(\"data/ne_10m_admin_0_countries/ne_10m_admin_0_countries.shp\", encoding=\"utf-8\")\n",
    "countries.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "countries2 = countries[[\"NAME\", \"POP_EST\", \"geometry\"]]\n",
    "countries2.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "by = countries2[countries2['NAME'] == \"Belarus\"]\n",
    "by"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "by.plot(figsize=(15,9))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "restoftheworld = countries2[countries2['NAME'] != \"Belarus\"]\n",
    "\n",
    "ax = by.plot(figsize=(15,9), color=\"green\")\n",
    "restoftheworld.plot(ax=ax, color=\"black\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "#type(de.geometry) # GeoSeries\n",
    "type(by.iloc[0].geometry) # shapely.geometry.multipolygon.MultiPolygon or Polygon"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Query positions within Polygons"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# Get All Cities within Belarus Polygon\n",
    "\n",
    "bycities = cities[cities.within(by.iloc[0].geometry)]\n",
    "bycities.head(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## 6. Folium"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "import folium\n",
    "\n",
    "\n",
    "m = folium.Map(location=[49.001575, 8.3832086], zoom_start=17)\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### GeoPandas & Folium\n",
    "\n",
    "our previous polygon \"by\" can be converted to GeoJson"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "import folium\n",
    "\n",
    "center = [49.001575, 8.3832086] # ZKM in Karlsruhe!\n",
    "map_zkm = folium.Map(center, zoom_start=5)   \n",
    "\n",
    "folium.GeoJson(by).add_to(map_zkm)\n",
    "\n",
    "map_zkm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "import folium\n",
    "\n",
    "center = [49.001575, 8.3832086] # ZKM in Karlsruhe!\n",
    "map_zkm = folium.Map(center, zoom_start=5)   \n",
    "\n",
    "folium.GeoJson(by,style_function=lambda feature: {\n",
    "        'fillColor': 'green',   # you can also replace this with functions with feature as argument\n",
    "        'color': 'black',\n",
    "        'weight': 2,\n",
    "        'dashArray': '5, 5'\n",
    "    }).add_to(map_zkm)\n",
    "\n",
    "map_zkm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# remember bigcities? Let's display them on a map using markers\n",
    "\n",
    "bigcities.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "world_map = folium.Map(location=[0,0], zoom_start=2)\n",
    "\n",
    "def f(row):\n",
    "    lng, lat = row[2].x, row[2].y\n",
    "    name = row[0]\n",
    "    population = row[1]\n",
    "    folium.Marker([lat, lng], \n",
    "              popup=name + \"<br/>\" + str(int(population)),\n",
    "              icon=folium.Icon(color=\"red\", prefix=\"fa\", icon=\"home\")).add_to(world_map)\n",
    "\n",
    "bigcities.apply(f, axis=1);\n",
    "    \n",
    "world_map\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "world_map.save(\"mymap.html\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Another Example: Geodetic Line"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "center = [53.869146, 27.487117] \n",
    "zoom = 4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# Geodetic Line from Minsk to New York\n",
    "import pyproj\n",
    "g = pyproj.Geod(ellps='WGS84')\n",
    "\n",
    "startlong = 27.487117\n",
    "startlat = 53.869146\n",
    "endlong = -74.001457\n",
    "endlat = 40.7094328\n",
    "\n",
    "lonlats = g.npts(startlong, startlat, endlong, endlat, 14)\n",
    "lonlats = [(startlong, startlat)] + lonlats +  [(endlong, endlat)]    # Add Start and End"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "world_map = folium.Map(location=center, zoom_start=zoom)\n",
    "\n",
    "for pos in lonlats:\n",
    "    folium.Marker([pos[1], pos[0]], \n",
    "                  icon=folium.Icon(color=\"red\", prefix=\"fa\", icon=\"airport\")).add_to(world_map)\n",
    "\n",
    "    \n",
    "world_map"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## More Geo?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "![gp2019](img/geopython.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Don't miss EP2019!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "![ep2019](img/ep2019.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Questions\n",
    "\n",
    "\n",
    "Download this Notebook here:  https://github.com/martinchristen/PyConBy2019\n",
    "\n",
    "\n",
    "mailto:martin.christen@fhnw.ch<br/>\n",
    "Twitter: @MartinChristen<br/>\n",
    "LinkedIn: https://www.linkedin.com/in/martinchristen/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}