{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Working with Shapefiles and OSM Part II of II\n", "\n", "**Author**: [Erika Fille Legara](http://www.erikalegara.net)\n", "\n", "You are free to use (or change) this notebook for any purpose you'd like. However, please respect the [MIT License](https://github.com/eflegara/PythonMaps/blob/master/LICENSE) that governs its use, and for copying permission.\n", "\n", "Copyright © 2016 Erika Fille Legara\n", "\n", "---\n", "## Description\n", "\n", "This recipe is the second of a series of Python notebooks on shapefiles and OSM files that I am working on. After going through the notebook, you should be able to plot maps as in the images below. \n", "\n", "![Metro Manila Maps](./img/maps.png \"Metro Manila Maps\")\n", "\n", "\n", "First things first, import the necessary packages. The [`Basemap` package](http://matplotlib.org/basemap/) is used for drawing and plotting maps (and also for reading shapefiles). The [`shapefile` package](https://pypi.python.org/pypi/pyshp/1.1.4) provides read and write support for the ESRI Shapefile format. " ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [], "source": [ "try:\n", " from mpl_toolkits.basemap import Basemap\n", " from matplotlib.patches import Polygon\n", " from matplotlib.collections import PatchCollection\n", " from matplotlib.patches import PathPatch\n", " import shapefile\n", " import os.path\n", " import matplotlib.pyplot as plt\n", " import urllib\n", " %matplotlib inline\n", " import numpy as np\n", "except:\n", " import traceback\n", " traceback.print_exc()\n", " raise ImportError('Something failed, see above.')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Loading the shapefile for cities\n", "The shapefile used here (for the Philippines and its municipalities) is obtained from [GADM](http://www.gadm.org/), which provides other shapefiles for other countries. The next code cell is only necessary to obtain the boundaries of a given administrative region. If you already have the boundary of your region of interest, you may skip the next three code cells. " ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [], "source": [ "phl2 = shapefile.Reader(\"PHL_adm_shp/PHL_adm2\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Fields (Records) + Geometries\n", "\n", "Now, let's pull out the shapes and records and obtain the bounding-box for each city. I just want to reiterate that the code cell below is not necessary if you already have a bounding-box in mind." ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "[('Tagana-An',\n", " [125.52970123291061, 9.65027046203636, 125.70527648925791, 9.769089698791568]),\n", " ('Talitay',\n", " [124.32147979736328, 6.983870029449747, 124.41300964355484, 7.08224010467552]),\n", " ('Dumalneg',\n", " [120.79208374023438, 18.415330886840877, 120.88435363769545, 18.531850814819375]),\n", " ('Aloguinsan',\n", " [123.5244369506836, 10.143600463867188, 123.63102722167966, 10.248220443725543]),\n", " ('Bongabong',\n", " [121.24404144287121, 12.618599891662598, 121.55500030517587, 12.798540115356673])]" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "citiesRecs = phl2.shapeRecords()\n", "cities_bbox = {}\n", "\n", "for entry in citiesRecs:\n", " cities_bbox[entry.record[6]] = entry.shape.bbox\n", "\n", "cities_bbox.items()[0:5]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can see from the previous cell that for each city we have the coordinates of their bounding-box. For example, for \"Tagana-An\", the lower-left corner of the bounding-box is at latitude and longitude 9.65027046203636 and 125.52970123291061, respectively; while the upper-right corner has (9.769089698791568, 125.70527648925791). \n", "\n", "The first 20 cities in the shapefile are listed below." ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "['Tagana-An', 'Talitay', 'Dumalneg', 'Aloguinsan', 'Bongabong', 'Villaviciosa', 'Solsona', 'Infanta', 'Malabon', 'Barobo', 'Minalin', 'Caramoran', 'San Jacinto', 'Mallig', 'Santa Ignacia', 'Tabuk City', 'Pandi', 'Tongkil', 'General Macarthur', 'General Emilio Aguinaldo']\n" ] } ], "source": [ "print cities_bbox.keys()[0:20]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## OpenStreetMap\n", "[Downloading OSM data](http://wiki.openstreetmap.org/wiki/Downloading_data) is straightforward now that we have the bounding boxes." ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [], "source": [ "getfile = urllib.URLopener()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Download OSM for cities\n", "\n", "Let's get the OSM data for each of the cities in Metro Manila. " ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [], "source": [ "metro_manila = [\"Manila\", \"Quezon City\", \"Pateros\", \"Kalookan City\", \"Las Piñas\", \"Makati City\", \n", " \"Malabon\", \"Mandaluyong\", \"Marikina\", \"Muntinlupa\", \"Navotas\", \"Parañaque\", \"Pasay City\", \n", " \"Pasig City\", \"Taguig\", \"Valenzuela\"]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "I use [XAPI](http://wiki.openstreetmap.org/wiki/XAPI) to download city information. According to [OSM's wiki page](http://wiki.openstreetmap.org/wiki/Downloading_data#Huge_amounts_of_data):\n", "\n", "
XAPI and Overpass API allow to download custom data sets like arbitrary bounding boxes, elements with specific tags, public transport networks or other features. [...] The API is limited to bounding boxes of about 0.5 degree by 0.5 degree and you should avoid using it if possible. For larger areas you might try to use XAPI [...]\n", "\n", "\n", "IMPORTANT NOTE: If the code chunk results to a 504 error, the issue is on the server side (OSM overpass); not the code. You may have to wait to download again. \n", "> If a query is rejected due to too much resource consumption, this is now answered with HTTP status code 504" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false, "scrolled": true }, "outputs": [], "source": [ "urlhead = \"http://overpass.osm.rambler.ru/cgi/xapi_meta?*[bbox=\"\n", "urltail = \"]\"\n", "\n", "for k in metro_manila:\n", " url = urlhead + str(cities_bbox[k][0]) + \",\" + str(cities_bbox[k][1]) + \",\" + \\\n", " str(cities_bbox[k][2]) + \",\" + str(cities_bbox[k][3]) + urltail\n", " \n", " if not os.path.exists(\"./cities_osm/\"):\n", " print \"Creating cities_osm directory...\"\n", " os.makedirs(\"./cities_osm/\")\n", " \n", " fname = \"./cities_osm/\"+ str(k)+\".osm\"\n", " \n", " if not os.path.isfile(fname): \n", " try:\n", " getfile.retrieve(url, fname)\n", " except IOError as err:\n", " print \"Whoops! Looks like the server is busy... consuming too much resource. Let's wait a bit.\" \n", " print \"Error: \", err\n", " break" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We now have the OSM files for the cities we are interested in. Let's see what is inside an OSM file. Let's have a look at the first few lines of one OSM file. " ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "