{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Geospatial analysis with Python\n",
"\n",
"Valencia - [PyConES 2015](http://2015.es.pycon.org/en/) - 2015/11/22\n",
"\n",
"- Pedro-Juan Ferrer · [@vehrka](http://twitter.com/vehrka)\n",
"- Jorge Sanz - [@xurxosanz](http://twitter.com/xurxosanz)\n",
"- Geoinquietos Valencia - [@geoinquietosvlc](http://twitter.com/geoinquietosvlc)\n",
"\n",
"Slides and repo\n",
"\n",
"- http://bit.ly/pycones2015-geo\n",
"- https://github.com/geoinquietosvlc/2015.es.pycon\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# 1. Introduction to Spatial Data\n",
"\n",
"the *Geographic Data problem*\n",
"\n",
"or\n",
"\n",
"**Spacial is Special**"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## GEO Problem?\n",
"\n",
"\n",
"\n",
"You have to deal with data that has special rules.\n",
"\n",
"Rules based on a reality not told in schools"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## Earth is not a sphere\n",
"\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## Earth is more like a potato\n",
"\n",
"\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## A potato\n",
"\n",
"Believe me, I know what I'm talking about\n",
"\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## Geoids are physical\n",
"\n",
"Although its the true shape of the Earth, we can't measure over it.\n",
"\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## We need mathematics\n",
"\n",
"An **Ellipsoid**\n",
"\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## Ellipsoid vs Geoid\n",
"\n",
"On average, they are quite similar\n",
"\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## Datum is the name of the game \n",
"\n",
"The Ellipsoid and a couple more of things is what we call:\n",
"\n",
"### the **DATUM**\n",
"\n",
"(And there are not one but several Datums, understanding this takes several courses of Geodesy, trust me on this)\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## Earth is not a sphere\n",
"\n",
"*Remember*\n",
"\n",
"> The **GEO** information lays over something *mathematical* that has its own rules.\n",
"\n",
"> **THE DATUM**"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Maps are flat\n",
"\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## Making things flat\n",
"\n",
"But you can't turn **flat** something spherical\n",
"\n",
"(without breaking it)\n",
"\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## Breaking Flat\n",
"\n",
"You have to choose what you want to *break*:\n",
"\n",
"* **Areas**\n",
"* **Angles**\n",
"* **Distances**\n",
"\n",
"In the best case you can choose two of the three.\n",
"\n",
"(I'm skipping several courses on Cartography with this slide, trust me on this)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## Choose your weapon!: The Projection\n",
"\n",
"Cartographers have tricks for breaking things **mathematically**\n",
"\n",
"[](https://xkcd.com/977/)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## Maps are flat\n",
"\n",
"*Remember*\n",
"\n",
"> The **GEO** information uses a *mathematical* trick to make things **flat** orderly.\n",
"\n",
"> **THE PROJECTION**"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## SRS and CRS\n",
"\n",
"Toghether a *DATUM* and a *PROJECTION* make a **CRS**\n",
"\n",
"The most famous catalog of *CRS* is **EPSG**\n",
"\n",
"* EPSG:4326\n",
"* EPSG:3857 or EPSG:900913\n",
"* EPSG:4258\n",
"* EPSG:25830, EPSG:25831"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Geo data types\n",
"\n",
"These one are easy\n",
"\n",
"Because you are more used to it"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## Raster and Vector data\n",
"\n",
"\n",
"\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## Data model\n",
"\n",
"How many geometry types do you know?\n",
"\n",
"... Point, Line, Surface ..."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## Well, it's complicated...\n",
"\n",
"OGC **Simple** Feature Access\n",
"\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## lon/lat vs lat/lon\n",
"\n",
"There's a little bit of fuss\n",
"\n",
"In theory it should be lon/lat (x,y)\n",
"\n",
"But we are used to lat/lon"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"**Always check things ... twice**"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Recap\n",
"\n",
"* Know your **Datum**\n",
"* Know your **Projection**\n",
"* Know your **Data Type**\n",
"* Know your **Model**"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"**KNOW YOUR (geo) DATA!!!**"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# 2. Writing and reading data"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## Fiona\n",
"- http://toblerity.org/fiona/\n",
"- Modern pythonic wrapper around [OGR](http://www.gdal.org/index.html)\n",
"- Simple and dependable\n",
"- Integrates with [pyproj](http://pypi.python.org/pypi/pyproj/), [Rtree](http://pypi.python.org/pypi/Rtree/) and [Shapely](http://pypi.python.org/pypi/Shapely/)"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"ARCGEN r\n",
"AeronavFAA r\n",
"BNA raw\n",
"DGN raw\n",
"DXF raw\n",
"ESRI Shapefile raw\n",
"FileGDB raw\n",
"GMT raw\n",
"GPKG rw\n",
"GPSTrackMaker raw\n",
"GPX raw\n",
"GeoJSON rw\n",
"Idrisi r\n",
"MapInfo File raw\n",
"OpenFileGDB r\n",
"PCIDSK r\n",
"PDS r\n",
"SEGY r\n",
"SUA r\n"
]
}
],
"source": [
"from fiona.collection import supported_drivers\n",
"for frmt in sorted(supported_drivers):\n",
" print(\"{:20}{}\".format(frmt,supported_drivers[frmt]))"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## Exercise\n",
"\n",
"1. Import a CSV and write into a geospatial file\n",
"2. Read that file and explore contents\n",
"3. Report"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Writing a shapefile from a CSV"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"1. reading a CSV from CartoDB [SQL API](http://docs.cartodb.com/cartodb-platform/sql-api/) (querying [this table](https://team.cartodb.com/u/jsanz/tables/twitter_pycones_pycones2015_pycones15/public))\n",
"2. creating a schema\n",
"3. creating features and writing them"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"\n",
"```sql\n",
"select \n",
" ST_X(the_geom) as lon, \n",
" ST_Y(the_geom) as lat, \n",
" cartodb_id, actor_preferredusername, body, postedtime \n",
"from \n",
" jsanz.twitter_pycones_pycones2015_pycones15\n",
"```"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": true,
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [],
"source": [
"import requests\n",
"url = 'https://jsanz.cartodb.com:443/api/v2/sql?q=select ST_X(the_geom) as lon, ST_Y(the_geom) as lat, cartodb_id, actor_preferredusername, body, postedtime from jsanz.twitter_pycones_pycones2015_pycones15&format=csv'\n",
"csv_file = '/tmp/tweets.csv'\n",
"\n",
"with open(csv_file,'w') as csvfile:\n",
" req = requests.get(url)\n",
" csvfile.write(req.text)"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": true,
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [],
"source": [
"target = '/tmp/tweets.shp'\n",
"epsg = 4258 # http://epsg.io/4258\n",
"driver = \"ESRI Shapefile\"\n",
"schema = {\n",
" \"geometry\": \"Point\",\n",
" \"properties\": {\n",
" (\"cartodb_id\", \"int\"),\n",
" (\"lon\",\"float\"),\n",
" (\"lat\",\"float\"),\n",
" (\"author\",\"str\"),\n",
" (\"body\",\"str\"),\n",
" (\"postedtime\",\"str\")\n",
" }\n",
"}"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": true,
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [],
"source": [
"import fiona, csv\n",
"from fiona.crs import from_epsg"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "-"
}
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"WARNING:Fiona:CPLE_AppDefined in b\"One or several characters couldn't be converted correctly from UTF-8 to ISO-8859-1.\\nThis warning will not be emitted anymore.\"\n"
]
}
],
"source": [
"output = fiona.open(target, \"w\", driver=driver, \n",
" crs=from_epsg(epsg), schema=schema)\n",
"with open(csv_file,'r') as csvfile:\n",
" csvreader = csv.reader(csvfile,delimiter=',',quotechar='\"')\n",
" next(csvreader) #skip the header\n",
" for line in csvreader:\n",
" try:\n",
" x = float(line[0])\n",
" y = float(line[1])\n",
" \n",
" feature = {\n",
" \"geometry\" : {\n",
" \"coordinates\" : (x, y), \"type\" : \"Point\"\n",
" },\n",
" \"properties\" : {\n",
" \"cartodb_id\" : int(line[2]),\n",
" \"lon\" : x,\n",
" \"lat\" : y,\n",
" \"author\" : line[3],\n",
" \"body\" : line[4],\n",
" \"postedtime\" : line[5]\n",
" }\n",
" }\n",
" output.write(feature)\n",
" except ValueError:\n",
" pass\n",
" try:\n",
" output.close()\n",
" except RuntimeError:\n",
" pass"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Reading and exploring data"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"1. Open a Shapefile\n",
"2. Getting information of the resource\n",
"3. Looping over features"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": true,
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [],
"source": [
"source = fiona.open(target, 'r')"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"text/plain": [
"'bounds close closed crs crs_wkt driver enabled_drivers encoding env filter flush guard_driver_mode items iterator keys meta mode name next path schema session validate_record validate_record_geometry values write writerecords'"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"\" \".join([atr for atr in dir(source) if atr[0] != '_'])"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"1031 tweets\r\n",
"\n",
"bounds: (-123.26204, -43.24895, 151.77647, 55.95206)\r\n",
"\n",
"CRS: {'no_defs': True, 'proj': 'longlat', 'ellps': 'GRS80'}\n"
]
}
],
"source": [
"print(\"{} tweets\\r\\n\".format(len(source)))\n",
"print(\"bounds: {}\\r\\n\".format(source.bounds))\n",
"print(\"CRS: {}\".format(source.crs))"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"ID - Author - Coords \r\n",
"*********************************************\n",
"187 - rmajadas - (-3.69063, 40.42526)\n",
" 39 - CValdeMontes - (2.15899, 41.38879)\n",
"249 - python_granada - (-3.60667, 37.18817)\n",
"247 - sdelquin - (-16.25462, 28.46824)\n",
"105 - ipedrazas - (-0.12574, 51.50853)\n",
"206 - rafbermudez - (-3.69063, 40.42526)\n",
"192 - rafbermudez - (-3.69063, 40.42526)\n",
" 9 - seattle113 - (-4.52406, 42.00955)\n",
" 35 - d1eg0_garc1a - (2.15899, 41.38879)\n",
" 99 - ipedrazas - (-0.12574, 51.50853)\n"
]
}
],
"source": [
"print(\"{:3} - {:15} - {:^15}\\r\\n{:*^45}\".format(\"ID\",\"Author\",\"Coords\",\"\"))\n",
"for f in source[:10]:\n",
" print(\"{:3} - {:15} - {}\"\n",
" .format(f['properties']['cartodb_id'],\n",
" f['properties']['author'], \n",
" f['geometry']['coordinates']))"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Displaying the imported data"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"There are some libraries to plot geographical data\n",
"\n",
"- [Descartes](https://pypi.python.org/pypi/descartes/)\n",
"- [Cartopy](http://scitools.org.uk/cartopy/)\n",
"- [Folium](https://folium.readthedocs.org/)"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"collapsed": true,
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [],
"source": [
"import folium\n",
"basemap = r'http://{s}.basemaps.cartocdn.com/light_all/{z}/{x}/{y}.png'\n",
"map = folium.Map(location=[39.5,-2.5], \n",
" zoom_start=6, width=960, height=600,\n",
" tiles=basemap, attr='OpenStreetMap and Twitter')"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"for f in source:\n",
" x,y = f['geometry']['coordinates']\n",
" map.simple_marker(\n",
" [y,x], #lat/lon!!!\n",
" popup=f['properties']['body'])"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
""
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"map"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
""
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# 3. Processing"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## Shapely\n",
"\n",
"- http://toblerity.org/shapely/\n",
"- 2D geometry processing\n",
"- Pythonic wrapper for [GEOS](https://trac.osgeo.org/geos/)\n",
"- Agnostic of coordinate systems or data formats\n",
"\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## Exercise\n",
"\n",
"1. Define a point at PyConES venue\n",
"2. Create a buffer of a radius of 100 km around it\n",
"3. Find tweets inside that buffer\n",
"4. Report"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Create a point"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"We can use the [Well Known Text](https://en.wikipedia.org/wiki/Well-known_text?oldformat=true) format to define a Shapely geometry.\n",
"\n",
""
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"collapsed": true,
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [],
"source": [
"from shapely.wkt import loads\n",
"pycones = loads(\"POINT (-0.346713 39.482767)\")"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"image/svg+xml": [
""
],
"text/plain": [
""
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"pycones"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
""
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [
{
"data": {
"text/plain": [
"'almost_equals area array_interface array_interface_base boundary bounds buffer centroid contains convex_hull coords covers crosses ctypes difference disjoint distance empty envelope equals equals_exact geom_type geometryType has_z impl interpolate intersection intersects is_closed is_empty is_ring is_simple is_valid length overlaps project relate relate_pattern representative_point simplify svg symmetric_difference to_wkb to_wkt touches type union within wkb wkb_hex wkt x xy y z'"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"\" \".join([ atr for atr in dir(pycones) if atr[0] != '_'])"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"pycones is valid\r\n",
"\n",
"WKT: POINT (-0.346713 39.482767)\r\n",
"\n",
"SVG: \r\n",
"\n"
]
}
],
"source": [
"print('pycones {}\\r\\n'.format('is valid' \n",
" if pycones.is_valid else 'is not valid'))\n",
"print('WKT: {}\\r\\n'.format(pycones.wkt))\n",
"print('SVG: {}\\r\\n'.format(pycones.svg()))"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Create a buffer"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"- We cannot create a buffer of 100km around a geodetic point\n",
"- We need to compute the buffer on a projected CRS\n",
"\n",
" 1. Project the point to UTM coordinates\n",
" 2. Compute the buffer\n",
" 3. Project the buffer to lat/lon coordinates"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"Shapely provides a method to allow using [pyproj](http://jswhit.github.io/pyproj/) with any geometry"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {
"collapsed": true,
"slideshow": {
"slide_type": "-"
}
},
"outputs": [],
"source": [
"import pyproj\n",
"from functools import partial"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"project = partial(\n",
" pyproj.transform,\n",
" pyproj.Proj(init='epsg:4258'),\n",
" pyproj.Proj(init='epsg:25830')\n",
")\n",
"project_inv = partial(\n",
" pyproj.transform,\n",
" pyproj.Proj(init='epsg:25830'),\n",
" pyproj.Proj(init='epsg:4326')\n",
")"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"POINT (728199.1076832865 4373712.985082146)\n"
]
}
],
"source": [
"import shapely.ops\n",
"pycones_25830 = shapely.ops.transform(project,pycones)\n",
"print(pycones_25830)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"Once projected we are ready to define a Shapely Point object, compute the buffer and project it back to lon/lat coordinates"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {
"collapsed": true,
"slideshow": {
"slide_type": "-"
}
},
"outputs": [],
"source": [
"from shapely.geometry import Point"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"p = Point(pycones_25830)\n",
"pycones_buffer_25830 = p.buffer(100000)\n",
"pycones_buffer = shapely.ops.transform(project_inv,pycones_buffer_25830)"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"(-1.5089067784607482,\n",
" 38.582592551915646,\n",
" 0.8139462109021063,\n",
" 40.38277685271952)"
]
},
"execution_count": 24,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"pycones_buffer.bounds"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Intersect tweets"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"All ready to read the Shapefile, check if every feature is intersected by the buffer and fill a list of tuples with the distance and the tweet."
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {
"collapsed": true,
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [],
"source": [
"import fiona\n",
"from shapely.geometry import shape"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"118 tweets at less than 100km\n"
]
}
],
"source": [
"tweets = []\n",
"with fiona.open('/tmp/tweets.shp','r') as source:\n",
" for f in source:\n",
" geometry = shape(f['geometry'])\n",
" if pycones_buffer.intersects(geometry):\n",
" geometry_25830 = shapely.ops.transform(project,geometry)\n",
" distance = geometry_25830.distance(pycones_25830)\n",
" tweets.append((distance,f))\n",
"print(\"{} tweets at less than 100km\".format(len(tweets)))"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"**WARNING**: We need to transform to projected coordinates to get a distance in meters!"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"Let's sort the results and print the closest 10 tweets to PyConES venue"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {
"collapsed": true,
"slideshow": {
"slide_type": "-"
}
},
"outputs": [],
"source": [
"tweets = sorted(tweets, key=lambda tweet: tweet[0])"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "-"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Author - Distance - ID\r\n",
"*************************************\n",
"ch_doig - 30.17 - 666\n",
"ch_doig - 58.11 - 583\n",
"ch_doig - 60.24 - 601\n",
"manuerux - 379.37 - 718\n",
"luiyo - 514.44 - 779\n",
"xurxosanz - 514.44 - 780\n",
"vero4ka_ru - 514.44 - 783\n",
"sdelquin - 917.01 - 806\n",
"andres_sanchis - 3009.98 - 209\n",
"OFNblog - 3009.98 - 262\n"
]
}
],
"source": [
"print(\"{:17} - {:>10} - {:>4}\\r\\n{:*^37}\".format(\"Author\",\"Distance\",\"ID\",\"\"))\n",
"for distance,tweet in tweets[:10]:\n",
" print(\"{:17} - {:10.2f} - {:4}\"\n",
" .format(tweet['properties']['author'],\n",
" distance,tweet['properties']['cartodb_id']))"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Displaying results"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"Let's use again folium but rendering also the buffer and with different colour for tweets inside."
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"basemap = r'http://{s}.basemaps.cartocdn.com/light_all/{z}/{x}/{y}.png'\n",
"map = folium.Map(location=[39.5,-2.5], \n",
" zoom_start=6, width=960, height=600,\n",
" tiles=basemap, attr='OpenStreetMap and Twitter')\n",
"buffer_coords = [ [lat,lon] for lon,lat in pycones_buffer.boundary.coords]\n",
"map.line(locations= buffer_coords)\n",
"with fiona.open('/tmp/tweets.shp','r') as source:\n",
" for f in source:\n",
" geometry = shape(f['geometry'])\n",
" color = 'red' if pycones_buffer.intersects(geometry) else 'blue'\n",
" x,y = f['geometry']['coordinates']\n",
" map.simple_marker([y,x], marker_color=color,\n",
" popup=f['properties']['body'])\n"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
""
]
},
"execution_count": 30,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"map"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
""
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"# Other libraries\n",
"\n",
"- [GeoPandas](http://geopandas.org/index.html)\n",
"- [PySAL](http://pysal.readthedocs.org/)\n",
"\n",
"From major Open Source desktop GIS projects:\n",
"\n",
"- [PyQGIS](http://docs.qgis.org/testing/en/docs/pyqgis_developer_cookbook/)\n",
"- [gvSIG Scripting](http://docs.gvsig.org/plone/projects/gvsig-desktop/docs/user/gvsig-desktop-2-0-scripting)\n",
"- [GRASS Scripting](http://grasswiki.osgeo.org/wiki/GRASS_Python_Scripting_Library)\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# That's all folks!\n",
"\n",
"## Thanks! Questions?\n",
"\n",
"- Pedro-Juan Ferrer · [@vehrka](http://twitter.com/vehrka)\n",
"- Jorge Sanz - [@xurxosanz](http://twitter.com/xurxosanz)\n",
"- Geoinquietos Valencia - [@geoinquietosvlc](http://twitter.com/geoinquietosvlc)\n",
"\n",
"Slides and repo\n",
"\n",
"- http://bit.ly/pycones2015-geo\n",
"- https://github.com/geoinquietosvlc/2015.es.pycon"
]
}
],
"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.4.3+"
}
},
"nbformat": 4,
"nbformat_minor": 0
}