IPython par l'exemple dans le cadre des données cartographiques (principalement)
Thomas Gratier - 18/10/2014 - WebGeoDataVore http://webgeodatavore.com
Ipython est défini sur Wikipédia comme "un terminal interactif, ou shell, pour le langage de programmation Python qui propose des fonctionnalités telles que l'introspection, une syntaxe additionnelle, la complétion et un historique riche"
Introduction
IPython permet de manière interactive de manipuler Python d'une manière beaucoup plus dynamique les modules comme les données que la console Python par défaut.
3 versions à retenir:
- Ligne de commande
- avec Notebook (qui inclue un serveur web et permet de manipuler Python depuis le navigateur)
- avec QT (client lourd)
Installer IPython (avec le notebook
Sous Windows (non recommandé, vous allez avoir des problèmes sur des besoins avancés à cause de modules en C)
- Récupérez https://bootstrap.pypa.io/ez_setup.py puis l'exécuter
- Installer pip avec
easy_install pip - Installer IPython avec
pip install ipython[all]
Sous Linux comme Mac
On suppose que vous savez utiliser pip, il faut faire en ligne de commande
pip install ipython[all]
Les trucs sympas
La liste des commandes magiques
Les commandes qui font que le IPython, c'est bon
%quickref
%%file filename.txt
Hello World!
!cat filename.txt
%load filename.txt
Hello World!
Charger des extensions
Partager ou découvrir des notebooks
Le principe: si vous trouver un fichier de type ipynb (l'extension des notebooks) sur Github, vous pouvez voir le résultat qu'il donne quand il est exécuté http://nbviewer.ipython.org/github/ipython/ipython/blob/1.x/examples/notebooks/Part%205%20-%20Rich%20Display%20System.ipynb
Attention, les exemples avec bibliothèques exotiques ou avec rendus web passeront incognito.
Visualiser et analyser : différentes phases
Analyser comme visualiser de la donnée implique plusieurs phases :
- Acquérir de la donnée
- La transformer (formats)
- L'analyser
- Visualiser
Acquérir
Depuis des données publiques
Où chercher?
- http://data.gouv.fr
- http://publicdata.eu
- portails locaux comme http://data.nantes.fr
- INSEE http://insee.fr
- OpenStreetMap http://openstreetmap.org
# Exemple de données des stations météo (source empruntée à Mathieu Leplatre)
#!wget ftp://ftp.wmo.ch/wmo-ddbs/VolA_New/Pub9volA141110x.flatfile
# Ouverture du navigateur sur une des pages de données du portail OpenData de Nantes
# "Localisation et caractéristiques des colonnes d'apport volontaire aériennes"
import webbrowser
webbrowser.open('http://data.nantes.fr/donnees/detail/localisation-des-colonnes-aeriennes-de-nantes-metropole/', new=2)
Depuis d'autres sources
Grâce à des adresses ou depuis des coordonnées:
GeoPy http://geopy.readthedocs.org
Plusieurs sources de données possibles:
- OpenStreetMap Nominatim
- ESRI ArcGIS
- Google Geocoding API (V3)
- Baidu Maps
- Bing Maps API
- Yahoo! PlaceFinder
- GeoNames
- MapQuest
- OpenMapQuest
- OpenCage
- SmartyStreets
- geocoder.us
- GeocodeFarm
# Geocodage depuis une adresse
from geopy.geocoders import Nominatim
geolocator = Nominatim()
location = geolocator.geocode("11 Impasse Juton, Nantes")
print(location.address)
print(location.latitude, location.longitude)
print(location.raw)
# Reverse geocoding
location = geolocator.reverse("47.21614, -1.54981")
print(location.address)
print(location.latitude, location.longitude)
print(location.raw)
Grâce à de l'aspiration de sites avec du WebScrapping (attention, limite légale à la pratique)
- lxml http://lxml.de
- BeautifulSoup http://www.crummy.com/software/BeautifulSoup/
- Scrapy http://scrapy.org
- PyQuery https://pythonhosted.org/pyquery/
- CasperJS http://casperjs.org
BeautifulSoup
from BeautifulSoup import BeautifulSoup
import requests
r = requests.get('http://leafletjs.com/plugins.html')
soup = BeautifulSoup(r.content)
tables = soup.findAll('table')
for table in tables:
trs = table.findAll('tr')
for i,row in enumerate(trs):
if i:
print i,row.findAll('td')[0].text
Transformer / stocker
Spécial carto
- GDAL / OGR http://gdal.org
- Pyproj https://code.google.com/p/pyproj/
- Fiona http://toblerity.org/fiona/
- Rasterio https://github.com/mapbox/rasterio
- Shapely http://toblerity.org/shapely/
Génériques mais pouvant être utilisé avec de la cartographie ou bien seuls
- Pandas http://pandas.pydata.org(+ GeoPandas http://geopandas.org)
- NumPy http://www.numpy.org
- SGBD génériques avec SQL (PostgreSQL/PostGIS, Spatialite,...)
OGR / GDAL
# OGR2OGR (utilitaire en ligne de commande) depuis Ipython
!rm voronoi_stations.geojson
!ogr2ogr -f "GeoJSON" -lco ENCODING=UTF-8 voronoi_stations.geojson voronoi_stations.shp
Fiona
- Surcouche à OGR, la partie vecteur de OGR / GDAL, la bibliothèque de manipulation de formats de données géographiques
import fiona
import json
import logging
from pyproj import Proj, transform
features = []
wgs84=Proj("+init=EPSG:4326")
with fiona.collection("COLONNES_AERIENNES_NM.shp", "r") as source:
p_in = Proj(source.crs)
for feat in source:
try:
assert feat['geometry']['type'] == "Point"
new_coords = []
feat['geometry']['coordinates'] = transform(p_in, wgs84, *feat['geometry']['coordinates'])
features.append(feat)
except Exception, e:
logging.exception("Error processing feature %s:", feat['id'])
my_layer = {
"type": "FeatureCollection",
"features": features
}
with open("colonnes_aeriennes.geojson", "w") as f:
f.write(json.dumps(my_layer))
Rasterio
Nous le mentionnons mais il est très spécifique.
Il est principalement destiné aux traitements d'images:
*il manipule des bandes
- permet de déformer des images pour de la reprojection
- fait aussi de la conversion de formats
- peut être converti en array NumPy: tous les processus de traitement d'images de NumPy peuvent être utilisés
Son avantage principal: moins de "boilerplate" qu'avec GDAL
Shapely
C'est un module en Python basé sur Geos, une bibliothèque C++ dont la fonction principale est de pouvoir effectuer des opérations spatiales.
Quand on entend opération spatiale, il faut comprendre que dans les SIG, trois éléments de base:
- le point,
- la ligne,
- le polygone
Selon le cas, on a:
- besoin de construire des géométries
- créer des nouvelles géométries (par croisement)
- vérifier des relations entre des objets géométriques.
## Enable embedding Matplotlib in IPython
%matplotlib inline
# From http://py2gis.blogspot.com/2010/10/import-polygon-by-numpy.html
from osgeo import ogr
from matplotlib import pyplot
from numpy import asarray
from shapely.wkb import loads
#plot polygon only
source = ogr.Open("voronoi_stations.shp")
layer = source.GetLayerByName("voronoi_stations")
#print layer
fig = pyplot.figure(1, figsize=(8,4), dpi=150)
ax = fig.add_subplot(111)
def plot_poly(g):
a = asarray(g.exterior)
pyplot.plot(a[:,0], a[:,1])
while 1:
feature = layer.GetNextFeature()
if not feature:
break
geom = loads(feature.GetGeometryRef().ExportToWkb())
#print geom
if geom.geom_type == 'MultiPolygon':
for g in geom:
plot_poly(g)
if geom.geom_type == 'Polygon':
plot_poly(geom)
#pyplot.show()
Pandas
Pivot principal pour lire et écrire des données multidimensionnelles. Il inclue par exemple la gestion des séries temporelles. Vous pouvez penser "cube de données". Il se base sur des "data frames", un concept qui vient de R, un language orienté vers les statistiques.
Il fonctionne très bien en s'articulant avec des fichiers à plats, HDF5 comme des bases de données.
Pour la partie cartographie, il a une extension dénommée GeoPandas qui s'appuie principalement sur Fiona et Shapely. Pour ceux qui ont des besoins cartographiques alors qu'ils viennent plutôt du monde du calcul numérique, cela évite de repartir de zéro.
Analyser
Attention, même si nous inventorions une bonne partie des outils, notre métier nous amène rarement à dépasser de l'analyse descriptive. Si nous allons au delà avec par exemple de l'analyse spatiale, nous avons tendance à utiliser d'autres outils plus avancés (par exemple des bindings Python sur des logiciels SIG bureautiques genre QGIS ou Grass) ou bien à passer sur PostGIS (la cartouche spatiale du SGBD PostgreSQL)
Spécial cartographie
- Pysal https://pysal.readthedocs.org en particulier la partie "esda.mapclassify — Choropleth map classification"
Composant quasi agnostique vis à vis des cartes
- CSVKit https://csvkit.readthedocs.org
- NLTK http://www.nltk.org
- SciPy http://www.scipy.org
- SciKit-Learn http://scikit-learn.org
Visualiser
Trois sorties principales:
Sur Matplotlib ou des dérivés
- Matplotlib http://matplotlib.org
- Descartes https://bitbucket.org/sgillies/descartes/
- Basemap http://matplotlib.org/basemap/
- CartoPy http://scitools.org.uk/cartopy/
Sur des contenus web
Cas "bâtard":
- Mapnik
Matplotlib
Déjà vu conjointement avec Shapely
Un exemple ci-dessous dans un contexte non cartographique mais fortement lié car basé sur la démographie (tiré de https://commons.wikimedia.org/wiki/File:Iran-age-pyramid.svg)
# coding: utf-8
#!/usr/bin/env python
import matplotlib.pyplot as plt
import numpy as np
# first column is total number of people
# second column is percentage of the population
data = np.array([
[6232552, 8.29],
[5657791, 7.53],
[5671435, 7.55],
[6607043, 8.79],
[8414497, 11.20],
[8672654, 11.54],
[6971924, 9.28],
[5571018, 7.41],
[4906749, 6.53],
[4030481, 5.36],
[3527408, 4.69],
[2680119, 3.57],
[1862907, 2.48],
[1343731, 1.79],
[1119968, 1.49],
[913531, 1.22],
[919539, 1.22]])
labels = [
'0-4',
'5-9',
'10-14',
'15-19',
'20-24',
'25-29',
'30-34',
'35-39',
'40-44',
'45-49',
'50-54',
'55-59',
'60-64',
'65-69',
'70-74',
'75-79',
'80 and over']
# Prepare axis
fig = plt.figure(figsize=(9, 5))
ax = fig.add_subplot(111)
y_pos = np.arange(len(labels))
percentages = data[:, 1]
ax.barh(y_pos, percentages, align='center')
plt.yticks(y_pos, labels)
plt.axis('tight')
plt.title('Iranian Demography 2011')
plt.xlabel('Percentage')
plt.savefig('iran-demography.svg', transparent=True)
# Avec NetworkX, on génère une image
!python napoleon_russian_campaign.py
!cat napoleon_russian_campaign.py
from IPython.display import Image
Image(filename='napoleon_russian_campaign.png')
Basemap
Il nécessite de reprendre les informations des stations téléchargées antérieurement: elles ont des valeurs de type degrés, minutes, secondes (DMS). Les outils utilisés attendent des degrés décimaux.
import csv
def dms2dec(value):
"""
Degres Minutes Seconds to Decimal degres
"""
degres, minutes, seconds = value.split()
seconds, direction = seconds[:-1], seconds[-1]
dec = float(degres) + float(minutes)/60 + float(seconds)/3600
if direction in ('S', 'W'):
return -dec
return dec
csv_file = 'Pub9volA141110x.flatfile'
with open('new_file.csv', 'wb') as csvfile:
reader = csv.DictReader(open(csv_file, 'rb'), delimiter="\t")
writer = csv.writer(csvfile, delimiter=';',
quotechar='"', quoting=csv.QUOTE_MINIMAL)
for line in reader:
lng = dms2dec(line.pop('Longitude'))
lat = dms2dec(line.pop('Latitude'))
wmoid = int(line.pop('StationId'))
name = line.pop('StationName').title()
writer.writerow([wmoid, name, lng, lat])
import numpy as np
# On ne garde que les colonnes numériques
lat, lon = np.loadtxt('new_file.csv', delimiter=';', unpack=True, usecols=[2, 3])
print lat
from mpl_toolkits.basemap import Basemap
from matplotlib import pyplot
#basemap_graph = Basemap(projection='robin', lat_0=45.0, lon_0=0.0)
basemap_graph = Basemap(projection='ortho',lon_0=0,lat_0=45,resolution='l')
## Enable embedding Matplotlib in IPython
%matplotlib inline
pyplot.figure(figsize=(16,16))
basemap_graph.drawcoastlines()
basemap_graph.drawcountries()
basemap_graph.fillcontinents()
basemap_graph.drawmeridians(np.arange(-180,180,20))
basemap_graph.drawparallels(np.arange(-90,90,20))
x, y = basemap_graph(lat, lon)
basemap_graph.plot(x, y, 'o', markersize=4, color='red')
pyplot.show()
CartoPy
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from cartopy.io.shapereader import Reader
from cartopy.feature import ShapelyFeature
fname = 'voronoi_stations.shp'
plt.figure(figsize=(18,18))
ax = plt.axes(projection=ccrs.PlateCarree())
ax.stock_img()
shape_feature = ShapelyFeature(Reader(fname).geometries(),
ccrs.PlateCarree(), facecolor='none')
ax.add_feature(shape_feature)
plt.show()
Folium
Création simple de cartes web depuis la ligne de commandes
Pratique pour l'exploration de données rapide
import folium
from IPython.display import HTML
map_1 = folium.Map(location=[47.2199, -1.5499], zoom_start=11,
tiles='Stamen Toner')
map_1.geo_json(geo_path='colonnes_aeriennes.geojson')
map_1.create_map(path='terrain.html')
# Open Browser in a new tab
import webbrowser
webbrowser.open('terrain.html', new=2)
from IPython.display import HTML
leaflet_file = open("leaflet_geojson.html").read()
HTML(leaflet_file)
%%javascript
console.log(map._layers)
from pygments_magic import *
%highlight_html leaflet_geojson.html
Mapnik
Un exemple avec libellé et contours départementaux
#!/usr/bin/env python
# -*- coding: UTF-8 -*-
# get extent from shp (for mapnik)
from ogr_extent import *
# Retrieve extent and projection using gdal
shp_extent , proj4 = extent_and_proj('departement.shp', sourcetype = 'ESRI Shapefile')
import cairo
from mapnik import Style, Rule, Color, Filter, LineSymbolizer, PolygonSymbolizer, TextSymbolizer, label_placement, Shapefile, SQLite, Layer, Map, render, Shapefile, Expression, save_map
map_output = 'france'
m = Map(300, 300, proj4)
m.background = Color('white')
t = TextSymbolizer(Expression('[CODE_DEPT]'), 'DejaVu Sans Book', 8, Color('black'))
f = Expression("[CODE_DEPT]<>'75' and [CODE_DEPT]<>'92' and [CODE_DEPT]<>'93' and [CODE_DEPT]<>'94'")
t.allow_overlap = 1
t.label_placement = label_placement.POINT_PLACEMENT
s1 = Style()
r1 = Rule()
r1.symbols.append(t)
r1.filter = f
s1.rules.append(r1)
m.append_style('Text', s1)
s = Style()
r = Rule()
r.symbols.append(PolygonSymbolizer(Color('#f2eff9')))
r.symbols.append(LineSymbolizer(Color('rgb(50%,50%,50%)'), 0.1))
s.rules.append(r)
m.append_style('My Style', s)
lyr = Layer('france', proj4)
import os
shp = Shapefile(base='.',file='departement')
lyr.datasource = shp
#lyr.datasource = SQLite(base=os.getcwd(), file = sqlitedatabase, table = tablename, geometry_field = 'Geometry', key_field = 'ID_GEOFLA', extent = shp_extent, wkb_format = 'spatialite')
lyr.styles.append('My Style')
lyr.styles.append('Text')
m.layers.append(lyr)
m.zoom_to_box(lyr.envelope())
file_formats = {'svg': cairo.SVGSurface,
}
for type_format in file_formats:
print '// -- Rendering %s -----------------------------' % type_format
file_open = open('%s.%s' % (map_output, type_format), 'w')
surface = file_formats[type_format](file_open.name, m.width, m.height)
render(m, surface)
save_map(m,"tmp_map.xml")
surface.finish()
from IPython.display import SVG
SVG(filename='france.svg')
# Pour aller plus loin, quelques bibliothèques non citées
from IPython.display import SVG
SVG(filename='pydata_stack_columns.svg')
Côté carto à tester (non testé mais avec un existant)
- Bokeh permet de faire des cartes http://bokeh.pydata.org/docs/gallery/texas.html
- Vincent (basé sur D3) http://vincent.readthedocs.org/en/latest/quickstart.html#map-data-binding (peut être utilisé par Folium)
- Pygal http://pygal.org/chart_types/#idworldmap-charts
Autres (pour mémoire/information):