In [36]:
%matplotlib inline

import matplotlib
import numpy as np
import matplotlib.pyplot as plt
matplotlib.rcParams['figure.figsize'] = (15.0, 8.0)
from matplotlib.pyplot import cm 
from IPython.display import Image, display
import grass.script as gscript

Compute difference 3D raster and query in predefined places:


In [37]:
gscript.run_command('g.region', raster3d='map3d_3760', overwrite=True)
gscript.run_command('r3.mapcalc', expression='map3d_37601_diff = map3d_3760 - map3d_3761', overwrite=True)
gscript.run_command('r3.what', overwrite=True, input='map3d_3761', points='POI',
                    output='POI_before.csv', separator='comma', null_value='0')
gscript.run_command('r3.what', overwrite=True, input='map3d_3760', points='POI',
                    output='POI_after.csv', separator='comma', null_value='0')
gscript.run_command('r3.what', overwrite=True, input='map3d_37601_diff', points='POI',
                    output='POI_diff.csv', separator='comma', null_value='0')


Out[37]:
0

In [38]:
data_diff = np.loadtxt('POI_diff.csv', delimiter=',')
data_before = np.loadtxt('POI_before.csv', delimiter=',')
data_after = np.loadtxt('POI_after.csv', delimiter=',')

Load values and convert to number of people per 100 m$^2$ per 1 hour


In [39]:
POI = gscript.read_command('v.out.ascii', input='POI', columns='cat,name', format='point').strip()
POI_data_diff = {}
POI_data_before = {}
POI_data_after = {}
for line in POI.splitlines():
    x, y, cat, name = line.split('|')
    x = float(x)
    y = float(y)
    for line in data_diff:
        if abs(line[0] - x) < 1e-3 and abs(line[1] - y) < 1e-3:
            POI_data_diff[name] = 8 *line[2:] # 8 - from 5m x 5m x 30 min to 10m x 10m x 1h
            
    for line in data_before:
        if abs(line[0] - x) < 1e-3 and abs(line[1] - y) < 1e-3:
            # 8 - from 5m x 5m x 30 min to 10m x 10m x 1h
            POI_data_before[name] = 8 *line[2:]
            
    for line in data_after:
        if abs(line[0] - x) < 1e-3 and abs(line[1] - y) < 1e-3:
            # 8 - from 5m x 5m x 30 min to 10m x 10m x 1h
            POI_data_after[name] = 8 *line[2:]

In [40]:
time = np.arange(345, 1260, 30) / 60.
color=iter(cm.Set1(np.linspace(0, 1, len(POI_data_diff))))
keys =  POI_data_diff.keys()

In [41]:
color=iter(cm.Set1(np.linspace(0, 1, len(POI_data_before))))
for key in keys:
    c = next(color)
    plt.plot(time, POI_data_before[key], c=c, linewidth=3, label=key)
plt.legend(loc=0)
plt.title("Pedestrian density before reconstruction")
plt.ylabel("# persons per 100m$^2$ per 1 hour")
plt.xlabel("Time of the day in hours")
plt.ylim(0, 0.3)
plt.show()



In [42]:
color=iter(cm.Set1(np.linspace(0, 1, len(POI_data_after))))
for key in keys:
    c = next(color)
    plt.plot(time, POI_data_after[key], c=c, linewidth=3, label=key)
plt.legend(loc=0)
plt.title("Pedestrian density after reconstruction")
plt.ylabel("# persons per 100m$^2$ per 1 hour")
plt.xlabel("Time of the day in hours")
plt.ylim(0, 0.3)
plt.show()



In [43]:
color=iter(cm.Set1(np.linspace(0, 1, len(POI_data_after))))
for key in keys:
    c = next(color)
    plt.plot(time, POI_data_diff[key], c=c, linewidth=3, label=key)
plt.legend(loc=0)
plt.title("Change in pedestrian density")
plt.ylabel("# persons per 100m$^2$ per 1 hour")
plt.xlabel("Time of the day in hours")
plt.show()



In [44]:
display(Image(filename='POI_before.png'))
display(Image(filename='POI_after.png'))