from cartopy import crs
import pylab as plt
import pandas as pd
import numpy as np
import os
df = pd.read_pickle('data/traffic_preprocessed.pkl')
df.head()
hour | month | id_arc_trafic | AVG(debit) | AVG(taux) | lat | lon | lndebit | time | dlo25 | dla25 | dlo50 | dla50 | dlo75 | dla75 | ||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
id_arc_trafic | ||||||||||||||||
1 | 0 | 0 | 1 | 1 | 1188.714286 | 4.638207 | 48.859838 | 2.334242 | 4.121122 | 100 | 13 | 12 | 25 | 24 | 38 | 37 |
93 | 0 | 2 | 1 | 1261.829268 | 5.031917 | 48.859838 | 2.334242 | 4.170615 | 200 | 13 | 12 | 25 | 24 | 38 | 37 | |
186 | 0 | 3 | 1 | 1468.347826 | 4.730384 | 48.859838 | 2.334242 | 4.298984 | 300 | 13 | 12 | 25 | 24 | 38 | 37 | |
79 | 0 | 4 | 1 | 1106.285714 | 5.214193 | 48.859838 | 2.334242 | 4.062314 | 400 | 13 | 12 | 25 | 24 | 38 | 37 | |
165 | 0 | 5 | 1 | 1452.818182 | 6.927379 | 48.859838 | 2.334242 | 4.289851 | 500 | 13 | 12 | 25 | 24 | 38 | 37 |
df2 = df.groupby(['hour','month','dlo25','dla25']).agg({'AVG(debit)':'sum','AVG(taux)':'mean','lat':'mean','lon':'mean','time':'mean'})
df2.columns = df2.columns.get_level_values(0)
df2.head()
lat | AVG(taux) | lon | AVG(debit) | time | ||||
---|---|---|---|---|---|---|---|---|
hour | month | dlo25 | dla25 | |||||
0 | 1 | 0 | 10 | 48.865882 | 0.957055 | 2.418085 | 304.571429 | 100 |
1 | 8 | 48.874043 | 6.193939 | 2.411269 | 6582.333333 | 100 | ||
9 | 48.870609 | 6.066761 | 2.413245 | 3409.523810 | 100 | |||
10 | 48.866809 | 1.804781 | 2.413973 | 767.428571 | 100 | |||
11 | 48.864982 | 4.341600 | 2.414442 | 4801.452381 | 100 |
# Thresholding for better visualization
trh = 15000
df2.loc[df2['AVG(debit)']>trh,'AVG(debit)']=trh
# Save extents
DLO = df.lon.max(),df.lon.min()
DLA = df.lat.max(),df.lat.min()
%pylab inline
from matplotlib.font_manager import FontProperties
from cartopy.io.img_tiles import OSM,StamenTerrain,GoogleTiles
from scipy.interpolate import griddata
from tqdm import tqdm_notebook as tqdm
font0 = FontProperties()
font0.set_family('serif')
font0.set_name('ubuntu')
Populating the interactive namespace from numpy and matplotlib
/home/astyonax/.anaconda/lib/python2.7/site-packages/IPython/core/magics/pylab.py:161: UserWarning: pylab import has clobbered these variables: ['plt'] `%matplotlib` prevents importing * from pylab and numpy "\n`%matplotlib` prevents importing * from pylab and numpy"
class StamenToner(GoogleTiles):
def _image_url(self, tile):
x, y, z = tile
url = 'http://tile.stamen.com/toner/{}/{}/{}.png'.format(z, x, y)
return url
imagery = OSM()
imagery = StamenToner()
# imagery = StamenTerrain()
# imagery = GoogleTiles()
Nota bene Plotting with the background map is slow, too slow to render each frame.
The solution is to render the background map once for all and store it to a file (traffic_frames/mapST.jpg) in my case.
Then the actual points are rendered on another figure, and saved with transparent background (see code, it's not super intuitive because you have to set multiple layers to be transparent). Since every figure and axis is created with the same properties (figsize, extent, and dpi) the actual overlay is trivial and performed on the fly by a call to convert
.
I use os.system
to KISS.
Finally, because we will use ffmpeg
to transform the resulting frames into a movie, the frame filename follows a strict naming code that is then feeded to ffmpeg
itself.
fig = plt.figure(figsize=(12, 12))
ax = plt.axes(projection=imagery.crs)
ax.set_extent((DLO[0],DLO[1],DLA[0],DLA[1]))
# # Add the imagery to the map.
ax.add_image(imagery, 13)
fig.savefig('traffic_frames/mapST.jpg',bbox_inches='tight',dpi=96)
plt.close();
month_names="January February March April May June July August September October November December".split(" ")
debit_min = df['AVG(debit)'].min()
debit_max = df['AVG(debit)'].max()
norm = plt.Normalize(vmin=debit_min,vmax=debit_max)
fig = plt.figure(figsize=(12, 12))
ax = plt.axes(projection=imagery.crs)
dx = DLO[1]-DLO[0]
dy = DLA[1]-DLA[0]
eps = 0/100.
extent = (DLO[0]-eps*dx,DLO[1]+eps*dx,DLA[0]-eps*dy,DLA[1]+eps*dy)
for time in tqdm(sorted(df.time.unique())):
ax.set_extent(extent)
dftimed = df2[df2.time==time]
try:
curr_month = dftimed.index.get_level_values('month').values.mean()
curr_month = int(curr_month)
curr_hour = dftimed.index.get_level_values('hour').values.mean()
except:
curr_month = dftimed.month.mean()
curr_month = int(curr_month)
curr_hour = dftimed.hour.mean()
xs,ys,vs = dftimed[['lon','lat','AVG(debit)']].values.T
xi_x,xi_y = np.mgrid[DLO[0]:DLO[1]:100j,DLA[0]:DLA[1]:100j]
vs = norm(vs)
vi = griddata(zip(xs,ys),vs,(xi_x,xi_y),method='linear').T
vs = plt.Normalize()(vs)
ax.imshow(vi,extent=extent,
transform=crs.PlateCarree(),cmap=plt.cm.jet,
alpha=.25,zorder=10)
ax.scatter(xs, ys, transform=crs.PlateCarree(), c=plt.cm.jet(vs), s=vs*200,lw=0,zorder=11)
ax.text(0.85,0.7,"{1:.0f}:00 {0:s}".format(month_names[curr_month-1],curr_hour),
fontsize=18, transform=ax.transAxes,va='center',ha='right', bbox=dict(facecolor='w', alpha=0.75),
fontproperties=font0,zorder=11)
ax.text(0,0,"Copyright (C) 2017 Guglielmo Saggiorato - map: Stamen Toner - data: opendata.paris.fr",fontsize=10,
fontproperties=font0,transform=ax.transAxes,va='bottom',ha='left',zorder = 12)
# make invisibility coat for the background
ax.outline_patch.set_visible(False)
ax.background_patch.set_visible(False)
fout = 'fr{0:04d}'.format(time)
fig.savefig('traffic_frames/{0:s}.png'.format(fout),bbox_inches='tight',dpi=96,transparent=True)
os.system('convert traffic_frames/mapST.jpg traffic_frames/{0:s}.png -layers merge traffic_frames/combined/{0:s}.jpg'.format(fout))
ax.cla()
# break
plt.close()
cmd = "rm {fout:s}; cat {indir:s}/*.jpg | ffmpeg -f image2pipe -r {rate:d} -i - -vf scale={xscale:d}:-1 -c:v libvpx-vp9 -b:v 2M {fout:s}"
cm = cmd.format(fout='traffic_25.webm',indir='traffic_frames/combined/',rate=10,xscale=640)
os.system(cm)
0
%%HTML
<center>
<video width="640" controls>
<source src="traffic_25.webm" type="video/webm">
</video>
</center>