BinderNotebook.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# EuroSciPy 2018\n",
    "\n",
    "Martin Christen, August 2018<br/>\n",
    "martin.christen@fhnw.ch\n",
    "\n",
    "## Handling Geospatial Data using Python, Part I: Base Modules\n",
    "\n",
    "What we learn today:\n",
    "\n",
    "* Reading Vector File Formats (Shapefiles, GeoJSON, ...)\n",
    "* Reading Raster File Formats (GeoTiff, ...)\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",
    "* Creating maps using Basemap & Folium ( https://matplotlib.org/basemap/ & https://github.com/python-visualization/folium )\n",
    "\n",
    "In this notebook the following data is used:\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",
    "* EarthQuake Data: https://earthquake.usgs.gov/earthquakes/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Setup\n",
    "\n",
    "To run this tutorial in Binder, the following modules need to be installed into the environment. This can take a while.\n",
    "\n",
    "### 1.1 Installing Modules\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 pillow\n",
    "    conda install folium -c conda-forge\n",
    "    conda install basemap -c conda-forge"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
		"! conda install shapely -y\n",
		"! conda install fiona -y\n",
		"! conda install rasterio -y\n",
		"! conda install geopandas -y\n",
		"! conda install pillow -y\n",
		"! conda install folium -y -c conda-forge\n",
		"! conda install basemap -y -c conda-forge"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Shapely\n",
    "\n",
    "http://toblerity.org/shapely\n",
    "\n",
    "http://toblerity.org/shapely/manual.html\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",
    "\n",
    "### 2.1 Simple Feature Access\n",
    "\n",
    "http://www.opengeospatial.org/standards/sfa\n",
    "\n",
    "https://en.wikipedia.org/wiki/Simple_Features\n",
    "\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",
    "Shapely supports the following Features:\n",
    "\n",
    "* Point\n",
    "* LineString\n",
    "* LinearRing          \n",
    "* Polygon\n",
    "* MultiLineString\n",
    "* MultiPoint\n",
    "* MultiPolygon\n",
    "\n",
    "\n",
    "### 2.2 Spatial Data Model\n",
    "\n",
    "The fundamental types of geometric objects implemented by Shapely are points, curves, and surfaces. Each is associated with three sets of (possibly infinite) points in the plane. The interior, boundary, and exterior sets of a feature are mutually exclusive and their union coincides with the entire plane\n",
    "\n",
    "* A Point has an interior set of exactly one point, a boundary set of exactly no points, and an exterior set of all other points. A Point has a topological dimension of 0.\n",
    "\n",
    "* A Curve has an interior set consisting of the infinitely many points along its length (imagine a Point dragged in space), a boundary set consisting of its two end points, and an exterior set of all other points. A Curve has a topological dimension of 1.\n",
    "\n",
    "* A Surface has an interior set consisting of the infinitely many points within (imagine a Curve dragged in space to cover an area), a boundary set consisting of one or more Curves, and an exterior set of all other points including those within holes that might exist in the surface. A Surface has a topological dimension of 2."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.3 Examples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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": {},
   "source": [
    "You can display shapely objects directly in the Jupyter Notebook. It is more or less a \"debug\" output.\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": {},
   "outputs": [],
   "source": [
    "polygon1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "polygon2 = Polygon([(20,20),(80,30),(50,40),(20,20)])\n",
    "polygon2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "polygon2.union(polygon1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "polygon2.intersection(polygon1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "polygon2.symmetric_difference(polygon1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "result = polygon2.symmetric_difference(polygon1)\n",
    "\n",
    "print(f\"Polygon area: {result.area}, polygon length: {result.length}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "result.wkt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s = result.wkt\n",
    "type(s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import shapely.wkt\n",
    "\n",
    "mypolygon = shapely.wkt.loads(s)\n",
    "mypolygon"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are also several binary operations available:\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": {},
   "outputs": [],
   "source": [
    "polygon1.intersects(polygon2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "polygon1.within(polygon2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "polygon1.equals(polygon1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Basemap & Spatial Reference Systems\n",
    "\n",
    "Example: https://www.openstreetmap.org/\n",
    "\n",
    "### 3.1 Coordinate Transformations using pyProj"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pyproj"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wgs84=pyproj.Proj(\"+init=EPSG:4326\")\n",
    "webmercator=pyproj.Proj(\"+init=EPSG:3857\")\n",
    "\n",
    "# EPSG Codes: http://epsg.io/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r0 = pyproj.transform(wgs84, webmercator, 11.1155568, 46.064887) \n",
    "\n",
    "print(r0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r1 = pyproj.transform(webmercator, wgs84, 1237378.122859708, 5790753.511617366)\n",
    "\n",
    "print(r1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2 Basemap Module"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from mpl_toolkits.basemap import Basemap\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(9, 15))\n",
    "\n",
    "m = Basemap(projection='cyl')     # llcrnrlat=-90,urcrnrlat=90, llcrnrlon=-180,urcrnrlon=180\n",
    "\n",
    "m.drawcoastlines(color=\"black\")\n",
    "m.fillcontinents(color='coral',lake_color='aqua')\n",
    "m.drawmapboundary(fill_color='aqua')\n",
    "\n",
    "plt.title(\"Equidistant Cylindrical Projection\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are a lot of different map projections available:\n",
    "\n",
    "https://matplotlib.org/basemap/users/mapsetup.html\n",
    "\n",
    "\n",
    "Map Projections can be:\n",
    "\n",
    "* conformal - preserve angles (locally)\n",
    "* equal-area - preserve area\n",
    "* equidistant - preserve distance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(9, 15))\n",
    "\n",
    "m = Basemap(projection='cyl')     # llcrnrlat=-90,urcrnrlat=90, llcrnrlon=-180,urcrnrlon=180\n",
    "\n",
    "m.bluemarble() # https://visibleearth.nasa.gov/\n",
    "\n",
    "plt.title(\"Equidistant Cylindrical Projection (Blue Marble)\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(9, 15))\n",
    "\n",
    "m = Basemap(projection='cyl')     # llcrnrlat=-90,urcrnrlat=90, llcrnrlon=-180,urcrnrlon=180\n",
    "\n",
    "m.shadedrelief() # http://www.shadedrelief.com/\n",
    "\n",
    "plt.title(\"Equidistant Cylindrical Projection (Shaded Relief)\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(9, 15))\n",
    "\n",
    "m = Basemap(projection='cyl')     # llcrnrlat=-90,urcrnrlat=90, llcrnrlon=-180,urcrnrlon=180\n",
    "\n",
    "m.etopo() # https://www.ngdc.noaa.gov/mgg/global/global.html\n",
    "# ETOPO1 is a 1 arc-minute global relief model of Earth's surface that integrates \n",
    "# land topography and ocean bathymetry\n",
    "\n",
    "plt.title(\"Equidistant Cylindrical Projection (etopo)\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(9, 15))\n",
    "\n",
    "m = Basemap(projection='cyl', llcrnrlat=30,urcrnrlat=50, llcrnrlon=0,urcrnrlon=20)\n",
    "\n",
    "m.etopo()\n",
    "\n",
    "# plot using matplotlib...\n",
    "m.plot(11.1155568, 46.064887, \"ro\", markersize=8) # Trento\n",
    "plt.text(11.1155568, 46.064887, \"Trento\", fontsize=16, color=\"white\")\n",
    "\n",
    "plt.title(\"Plot\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Other Projections"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(9, 9))\n",
    "m = Basemap(projection='merc',llcrnrlat=-85.15,urcrnrlat=85.15,llcrnrlon=-180,urcrnrlon=180,lat_ts=20,resolution='c')\n",
    "\n",
    "m.drawcoastlines()\n",
    "m.fillcontinents(color='coral',lake_color='aqua')\n",
    "m.drawmapboundary(fill_color='aqua')\n",
    "\n",
    "plt.title(\"Mercator Projection (Gerhard Mercator 1512-1594)\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Geodesic Line: Trento to New York\n",
    "\n",
    "https://en.wikipedia.org/wiki/Geodesic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g = pyproj.Geod(ellps='WGS84')\n",
    "\n",
    "startlong = 11.1155568\n",
    "startlat = 46.064887\n",
    "\n",
    "endlong = -74.001457\n",
    "endlat = 40.7094328\n",
    "\n",
    "lonlats = g.npts(startlong, startlat, endlong, endlat, 14)\n",
    "\n",
    "lonlats = [(startlong, startlat)] + lonlats +  [(endlong, endlat)]    # Append start and end...\n",
    "\n",
    "print(lonlats)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from mpl_toolkits.basemap import Basemap\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "plt.figure(figsize=(15, 9))\n",
    "m = Basemap(projection='cyl')     # llcrnrlat=-90,urcrnrlat=90, llcrnrlon=-180,urcrnrlon=180\n",
    "\n",
    "m.drawcoastlines()\n",
    "m.fillcontinents(color='coral',lake_color='aqua')\n",
    "m.drawmapboundary(fill_color='aqua')\n",
    "\n",
    "for coord in lonlats:\n",
    "    m.plot(coord[0], coord[1], 'ro', markersize=8)\n",
    "\n",
    "#m.drawgreatcircle(startlong,startlat,endlong,endlat,linewidth=2,color='b')\n",
    "   \n",
    "    \n",
    "plt.title(\"Trento to New York\")\n",
    "plt.show();\n",
    "\n",
    "print(m.proj4string)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(9, 9))\n",
    "# lon_0, lat_0 are the center point of the projection.\n",
    "# resolution = 'l' means use low resolution coastlines.\n",
    "m = Basemap(projection='ortho',lon_0=-30,lat_0=40,resolution='l')\n",
    "m.drawcoastlines()\n",
    "m.fillcontinents(color='coral',lake_color='aqua')\n",
    "\n",
    "m.drawmapboundary(fill_color='aqua')\n",
    "\n",
    "for coord in lonlats:\n",
    "    x,y = m(coord[0], coord[1])\n",
    "    m.plot(x,y, 'ro', markersize=8)\n",
    "\n",
    "m.drawgreatcircle(startlong,startlat,endlong,endlat,linewidth=2,color='b')\n",
    "        \n",
    "plt.title(\"Full Disk Orthographic Projection - Trento to New York\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Fiona\n",
    "\n",
    "https://github.com/Toblerity/Fiona\n",
    "\n",
    "Spatial information is not only geometry. Together with it we always have metadata (properties). Lets look at a first dataset.\n",
    "\n",
    "\n",
    "### 4.1 Reading Vector Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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": {},
   "outputs": [],
   "source": [
    "airport['properties']['name']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "airport['geometry']['type']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "airport['geometry']['coordinates']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "c.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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": {},
   "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'] == \"VCE\":\n",
    "            print(airport['properties']['name'])\n",
    "            print(airport['geometry']['coordinates'])\n",
    "            print(airport['properties']['wikipedia'])        \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Coordinate System ? (Coordinate Reference System?)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with fiona.open('data/ne_10m_airports/ne_10m_airports.shp', 'r') as c:\n",
    "    print(c.crs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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": {},
   "source": [
    "Lets look at another dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import fiona\n",
    "\n",
    "c = fiona.open('data/ne_10m_admin_0_countries/ne_10m_admin_0_countries.shp', 'r')\n",
    "\n",
    "country = next(iter(c))\n",
    "country;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(country['properties']['NAME'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(country['properties']['NAME_ZH'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(country['properties']['CONTINENT'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(country['properties']['POP_EST'])\n",
    "print(country['properties']['POP_YEAR'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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'] == \"Italy\":\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": {},
   "source": [
    "### 4.2 Example: find Airports in Italy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from shapely.geometry import shape\n",
    "from shapely.geometry import MultiPolygon\n",
    "\n",
    "it = MultiPolygon(shape(geometry))\n",
    "it"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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(it):\n",
    "            print(airport['properties']['iata_code'], airport['properties']['name'], airport['geometry']['coordinates'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.3 Example: find Populated Places in Italy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import fiona\n",
    "\n",
    "with fiona.open('data/ne_10m_populated_places_simple/ne_10m_populated_places_simple.shp', 'r') as c:\n",
    "    for place in c:\n",
    "        geom = place[\"geometry\"]\n",
    "        position = Point(place['geometry']['coordinates'])\n",
    "        name = place[\"properties\"][\"name\"]\n",
    "        pop = int(place[\"properties\"][\"pop_max\"])\n",
    "        if pop>100000 and position.within(it):\n",
    "            print(name, pop, position)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Rasterio\n",
    "\n",
    "https://github.com/mapbox/rasterio\n",
    "\n",
    "Based on GDAL: https://www.gdal.org/\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.1 Reading data\n",
    "\n",
    "    rasterio.open(fp, mode='r', driver=None, width=None, height=None, count=None, crs=None, transform=None, dtype=None, nodata=None, **kwargs)\n",
    "    \n",
    "    \n",
    "#### Parameters:\t\n",
    "**fp**: string or file\n",
    "  A filename or URL, or file object opened in binary mode.\n",
    "\n",
    "**mode**: string\n",
    "“r” (read), “r+” (read/write), or “w” (write)\n",
    "\n",
    "**driver**: string\n",
    "Driver code specifying the format name (e.g. “GTiff” or “JPEG”). See GDAL docs at http://www.gdal.org/formats_list.html (optional, required for writing).\n",
    "\n",
    "**width**: int\n",
    "Number of pixels per line (optional, required for write).\n",
    "\n",
    "**height**: int\n",
    "Number of lines (optional, required for write).\n",
    "\n",
    "**count**: int > 0\n",
    "Count of bands (optional, required for write).\n",
    "\n",
    "**dtype**: rasterio.dtype\n",
    "the data type for bands such as rasterio.ubyte for 8-bit bands or rasterio.uint16 for 16-bit bands (optional, required for write)\n",
    "\n",
    "**crs**: dict or string\n",
    "Coordinate reference system (optional, recommended for write).\n",
    "\n",
    "**transform**: Affine instance\n",
    "Affine transformation mapping the pixel space to geographic space (optional, recommended for writing).\n",
    "\n",
    "**nodata**: number\n",
    "Defines pixel value to be interpreted as null/nodata (optional, recommended for write, will be broadcast to all bands)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import rasterio\n",
    "\n",
    "dataset = rasterio.open('data/BlueMarble.tif', 'r')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset.name"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset.mode"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset.count # number of raster bands, in our case 3 for r,g,b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset.indexes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset.width"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset.height"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset.crs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset.affine  # affine transformation pixel to crs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset.affine * (0, 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "~dataset.affine # inverse affine transformation crs to pixel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "~dataset.affine * (0,0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "px,py = ~dataset.affine * (11.1155568, 46.064887) # Trento\n",
    "print(px,py)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset.bounds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(dataset.bounds.left)\n",
    "print(dataset.bounds.bottom)\n",
    "print(dataset.bounds.right)\n",
    "print(dataset.bounds.top)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.2 Plot using Matplotlib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r = dataset.read(1)\n",
    "g = dataset.read(2)\n",
    "b = dataset.read(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(15,9))\n",
    "ax.imshow(rgb, interpolation='nearest')\n",
    "ax.plot(px,py, 'ro'); "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.3 Creating Overviews\n",
    "\n",
    "Overviews are reduced resolution versions of the dataset that can speed up rendering when you don't need full resolution. \n",
    "\n",
    "By precomputing the upsampled pixels, rendering can be **significantly faster** when zoomed out."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import shutil\n",
    "\n",
    "path = shutil.copy('data/BlueMarble.tif', 'data/BlueMarbleResampling.tif')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import rasterio\n",
    "\n",
    "src = rasterio.open(\"data/BlueMarbleResampling.tif\", 'r')\n",
    "for i in src.indexes:\n",
    "    print(src.overviews(i))\n",
    "src.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from rasterio.enums import Resampling\n",
    "\n",
    "dst = rasterio.open(\"data/BlueMarbleResampling.tif\", 'r+')\n",
    "dst.build_overviews([2, 4, 8, 16], Resampling.average)\n",
    "dst.update_tags(ns='rio_overview', resampling='average') # optional\n",
    "dst.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "src = rasterio.open(\"data/BlueMarbleResampling.tif\", 'r')\n",
    "for i in src.indexes:\n",
    "    print(src.overviews(i))\n",
    "src.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "src = rasterio.open(\"data/BlueMarbleResampling.tif\", 'r')\n",
    "\n",
    "r = src.read(1, out_shape=(1, int(src.height / 4), int(src.width / 4)))\n",
    "g = src.read(2, out_shape=(1, int(src.height / 4), int(src.width / 4)))\n",
    "b = src.read(3, out_shape=(1, int(src.height / 4), int(src.width / 4)))\n",
    "\n",
    "rgb = np.dstack((r,g,b))  # stack r,g,b so we can display it...\n",
    "rgb\n",
    "\n",
    "src.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(15,9))\n",
    "ax.imshow(rgb, interpolation='nearest');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.4 Windowed read"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset = rasterio.open('data/BlueMarbleResampling.tif')\n",
    "\n",
    "w = ((256,512),(1024,1280))\n",
    "\n",
    "r = dataset.read(1, window=w)\n",
    "g = dataset.read(2, window=w)\n",
    "b = dataset.read(3, window=w)\n",
    "\n",
    "rgb = np.dstack((r,g,b))\n",
    "\n",
    "dataset.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(15,9))\n",
    "ax.imshow(rgb, interpolation='nearest');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6. Folium\n",
    "\n",
    "http://python-visualization.github.io/folium/\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import folium"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = folium.Map(location=[47.37825,8.5367835], zoom_start=14) \n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Different Maps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = folium.Map(location=[47,8],\n",
    "           tiles='Stamen Toner',\n",
    "           zoom_start=6)\n",
    "m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = folium.Map(location=[47,8],\n",
    "           tiles='Mapbox Bright',\n",
    "           zoom_start=6)\n",
    "m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = folium.Map(location=[47,8],\n",
    "           tiles='Mapbox Control Room',\n",
    "           zoom_start=6)\n",
    "m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = folium.Map(location=[47,8],\n",
    "           tiles='openstreetmap',\n",
    "           zoom_start=6)\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "More examples (for some you need API keys)\n",
    "\n",
    "http://deparkes.co.uk/2016/06/10/folium-map-tiles/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Adding Markers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = folium.Map(location=[47.37825,8.5367835], zoom_start=16)\n",
    "\n",
    "folium.Marker([47.37695,8.5387885], popup=\"Hotel Schweizerhof\").add_to(m)\n",
    "folium.Marker([47.376386,8.5386506], popup=\"Hotel St. Gotthard\").add_to(m)\n",
    "m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "help(folium.Icon)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = folium.Map(location=[47.37825,8.5367835], zoom_start=16)\n",
    "\n",
    "folium.Marker([47.37695,8.5387885], \n",
    "              popup=\"Hotel Schweizerhof\",\n",
    "              icon=folium.Icon(color=\"red\", prefix=\"fa\", icon=\"hotel\")).add_to(m)\n",
    "\n",
    "folium.Marker([47.376386,8.5386506], \n",
    "              popup=\"Hotel St. Gotthard\",\n",
    "              icon=folium.Icon(color=\"green\", prefix=\"fa\", icon=\"hotel\")).add_to(m)\n",
    "\n",
    "folium.Marker([47.376192, 8.540005], \n",
    "              popup=\"Hotel Townhouse\", \n",
    "              icon=folium.Icon(color=\"blue\", prefix=\"fa\", icon=\"hotel\")).add_to(m)\n",
    "\n",
    "m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = folium.Map(location=[0,0], zoom_start=2)\n",
    "\n",
    "red = \"red\"\n",
    "green = \"green\"\n",
    "\n",
    "# Remember: Populated Places in Italy:\n",
    "with fiona.open('data/ne_10m_populated_places_simple/ne_10m_populated_places_simple.shp', 'r') as c:\n",
    "    for place in c:\n",
    "        geom = place[\"geometry\"]\n",
    "        position = Point(place['geometry']['coordinates'])\n",
    "        name = place[\"properties\"][\"name\"]\n",
    "        pop = int(place[\"properties\"][\"pop_max\"])\n",
    "        if pop>100000 and position.within(it):\n",
    "            lng, lat = place['geometry']['coordinates']\n",
    "            if name == \"Trento\":\n",
    "                color = \"green\"\n",
    "            else:\n",
    "                color = \"red\"\n",
    "                \n",
    "            folium.Marker([lat, lng], \n",
    "                            popup=name + \"<br/>\" + str(pop),\n",
    "                           icon=folium.Icon(color=color, prefix=\"fa\", icon=\"building\")).add_to(m)\n",
    "\n",
    "\n",
    "m\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Adding GeoJson Shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "center = [46.064887, 11.1155568] \n",
    "m = folium.Map(center, zoom_start=6)  \n",
    "\n",
    "red = \"red\"\n",
    "green = \"green\"\n",
    "\n",
    "# Remember: Populated Places in Italy:\n",
    "with fiona.open('data/ne_10m_populated_places_simple/ne_10m_populated_places_simple.shp', 'r') as c:\n",
    "    for place in c:\n",
    "        geom = place[\"geometry\"]\n",
    "        position = Point(place['geometry']['coordinates'])\n",
    "        name = place[\"properties\"][\"name\"]\n",
    "        pop = int(place[\"properties\"][\"pop_max\"])\n",
    "        if pop>100000 and position.within(it):\n",
    "            lng, lat = place['geometry']['coordinates']\n",
    "            if name == \"Trento\":\n",
    "                color = \"green\"\n",
    "            else:\n",
    "                color = \"red\"\n",
    "                \n",
    "            folium.Marker([lat, lng], \n",
    "                            popup=name + \"<br/>\" + str(pop),\n",
    "                           icon=folium.Icon(color=color, prefix=\"fa\", icon=\"building\")).add_to(m)\n",
    "\n",
    "\n",
    "folium.GeoJson(it,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(m)\n",
    "\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 7. Live Data Example: Earthquakes\n",
    "\n",
    "https://earthquake.usgs.gov/earthquakes/feed/v1.0/geojson.php"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import requests\n",
    "\n",
    "url = \"https://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/2.5_week.geojson\"\n",
    "#url = \"https://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/significant_month.geojson\"\n",
    "#url = \"https://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/2.5_month.geojson\"\n",
    "\n",
    "data = requests.get(url)\n",
    "file = open(\"earthquakes.geojson\",\"wb\")\n",
    "file.write(data.content)\n",
    "file.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import fiona\n",
    "\n",
    "with fiona.open('earthquakes.geojson', 'r') as c:\n",
    "    for earthquake in c:\n",
    "        geom = earthquake[\"geometry\"]\n",
    "        position = Point(earthquake['geometry']['coordinates'])\n",
    "        magnitude = earthquake[\"properties\"][\"mag\"]\n",
    "        title = earthquake[\"properties\"][\"title\"]\n",
    "        \n",
    "        if magnitude >4.5:\n",
    "            lng = position.x\n",
    "            lat = position.y\n",
    "            elv = position.z\n",
    "            print(title, lng, lat, elv)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import html\n",
    "\n",
    "center = [0, 0] \n",
    "m = folium.Map(center, zoom_start=2)\n",
    "\n",
    "\n",
    "with fiona.open('earthquakes.geojson', 'r') as c:\n",
    "    for earthquake in c:\n",
    "        geom = earthquake[\"geometry\"]\n",
    "        position = Point(earthquake['geometry']['coordinates'])\n",
    "        magnitude = earthquake[\"properties\"][\"mag\"]\n",
    "        title = earthquake[\"properties\"][\"title\"]\n",
    "        \n",
    "        if magnitude >4.5:\n",
    "            lng = position.x\n",
    "            lat = position.y\n",
    "            elv = position.z\n",
    "            folium.CircleMarker([lat, lng], \n",
    "                                radius=magnitude,\n",
    "                                popup=html.escape(title, quote=True),\n",
    "                                fill=True,\n",
    "                                color = \"#ff0000\",\n",
    "                                fill_color='#ff0000',\n",
    "                                fill_opacity=1.0).add_to(m)\n",
    "                                \n",
    "m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}