In [59]:
import numpy as np
import pandas as pd
import geopandas as gpd
import seaborn as sns
import matplotlib.pyplot as plt
import mplleaflet as mpll
In [60]:
accidents = gpd.read_file("donnees/accidents-geobase/accidents_2018.shp")
accidents = accidents[~accidents.geometry.isnull()]
accidents.info()
In [61]:
sns.jointplot(x="LOC_X", y="LOC_Y", data=accidents)
Out[61]:
In [62]:
#accidents#
fig, ax = plt.subplots(1, figsize=(6, 6))
wgsAccidents = gpd.read_file("donnees/accidents-geobase/accidents0.shp").to_crs(epsg = 4326)
wgsAccidents.plot(ax=ax)
mpll.display(fig=fig)
Out[62]:
In [63]:
from pointpats import PointPattern
from pointpats.centrography import hull, mbr, mean_center, weighted_mean_center, manhattan_median, std_distance,euclidean_median,ellipse
pp = PointPattern(accidents[['LOC_X','LOC_Y']])
pp.summary()
pp10 = PointPattern(pp.points[::10])
In [64]:
mc = mean_center(pp.points)
print(mc)
pp.plot()
plt.plot(mc[0], mc[1], 'b^', label='Centre moyen')
plt.legend(numpoints=1)
plt.axis('equal')
Out[64]:
In [65]:
# centre moyen pondéré
wmc = weighted_mean_center(pp.points, accidents['NB_BLESSES']+accidents['NB_MORTS'])
# accidents[['NB_BLESSES','NB_MORTS']].sum(axis=1)
print(wmc)
pp.plot() #use class method "plot" to visualize point pattern
plt.plot(mc[0], mc[1], 'b^', label='Centre moyen')
plt.plot(wmc[0], wmc[1], 'gd', label='Centre moyen pondéré')
plt.legend(numpoints=1)
plt.axis('equal')
plt.figure()
pp.plot() #use class method "plot" to visualize point pattern
plt.plot(mc[0], mc[1], 'b^', label='Centre moyen', mec = 'k', ms = 10)
plt.plot(wmc[0], wmc[1], 'gd', label='Centre moyen pondéré', mec = 'k', ms = 10)
#plt.legend(numpoints=1)
plt.axis('equal')
plt.xlim(mc[0]-2000,mc[0]+2000)
plt.ylim(mc[1]-2000,mc[1]+2000)
Out[65]:
In [66]:
# Écart-type de la distance de chaque point par rapport au centre moyen
stdd = std_distance(pp.points)
print(stdd)
circle1=plt.Circle((mc[0], mc[1]),stdd,color='r')
ax = pp10.plot(get_ax=True, title='Standard Distance Circle')
ax.add_artist(circle1)
plt.plot(mc[0], mc[1], 'b^', label='Mean Center')
ax.set_aspect('equal')
plt.legend(numpoints=1)
plt.axis('equal')
Out[66]:
In [67]:
# Ellipse de dispersion
sx, sy, theta = ellipse(pp.points)
print(sx, sy, theta)
from matplotlib.patches import Ellipse
#fig, ax = plt.subplots(1, figsize = (10,10))
ax = pp10.plot(title='Ellipse de dispersion', get_ax = True)
e = Ellipse(xy=mean_center(pp.points), width=sx*2, height=sy*2, angle=-theta*180./np.pi) #angle is rotation in degrees (anti-clockwise)
ax.add_artist(e)
e.set_clip_box(ax.bbox)
e.set_facecolor([0.8,0,0])
e.set_edgecolor([1,0,0])
#ax.set_xlim(0,100)
#ax.set_ylim(0,100)
ax.set_aspect('equal')
plt.axis('equal')
plt.plot(mc[0], mc[1], 'b^', label='Mean Center')
plt.legend(numpoints=1)
Out[67]:
In [68]:
fig, ax = plt.subplots(1, figsize = (16,16))
for mois in accidents.MOIS.unique():
points = accidents.loc[accidents.MOIS == mois, ['LOC_X','LOC_Y']]
sx, sy, theta = ellipse(points)
print(sx, sy, theta)
c = mean_center(points)
e = Ellipse(xy=c, width=sx*2, height=sy*2, angle=-theta*180./np.pi) #angle is rotation in degrees (anti-clockwise)
ax.add_artist(e)
e.set_clip_box(ax.bbox)
e.set_facecolor([0.8,0,0])
e.set_alpha(0.2)
plt.plot(c[0], c[1], 'x',label='centre '+mois)
ax.set_xlim(mc[0]-1.5*sx,mc[0]+1.5*sx)
ax.set_ylim(mc[1]-1.5*sy,mc[1]+1.5*sy)
ax.set_aspect('equal')
#plt.plot(mc[0], mc[1], 'b^', label='Mean Center')
plt.legend()
Out[68]:
In [69]:
print(hull(pp.points), mbr(pp.points))
pp.plot(title='Centers', hull=True, window = True ) #plot point pattern "pp" as well as its convex hull
plt.plot(mc[0], mc[1], 'b^', label='Mean Center')
plt.plot(wmc[0], wmc[1], 'gd', label='Weighted Mean Center')
plt.legend()
Out[69]:
In [70]:
from pointpats import PoissonPointProcess, PoissonClusterPointProcess, Window, poly_from_bbox, PointPattern
import libpysal as ps
import pointpats.quadrat_statistics as qs
In [71]:
ile = ps.io.open('donnees/accidents-geobase/terre_shp.shp')
polys = [p for p in ile]
pmax=polys[0]
for p in polys:
if p.area>pmax.area:
pmax = p
window = Window(pmax.parts)
In [72]:
samples = PoissonPointProcess(window, 200, 1, conditioning=False, asPP=True)
In [73]:
pp_csr = samples.realizations[0]
pp_csr.plot(window=True, hull=True, title='Motif de points aléatoire')
pp_csr.summary()
#plt.plot(*pmax.exterior.xy)
In [74]:
pp.summary()
print(pp.lambda_mbb, pp.lambda_hull)
In [75]:
q_r = qs.QStatistic(pp,shape= "rectangle",nx = 20, ny = 10)
#plt.figure(figsize=(12,12))
q_r.plot()
print(q_r.chi2,q_r.df,q_r.chi2_pvalue)
In [76]:
q_r = qs.QStatistic(pp_csr,shape= "rectangle",nx = 10, ny = 5)
plt.figure(figsize=(12,12))
q_r.plot()
print(q_r.chi2,q_r.df,q_r.chi2_pvalue)
In [77]:
q_h = qs.QStatistic(pp_csr,shape= "hexagon",lh = 5000)
q_h.plot()
print(q_h.chi2,q_h.df,q_h.chi2_pvalue)
In [78]:
points = [[0., 0.], [0., 1.], [1., 1.], [1., 0.]]
samples2 = PoissonPointProcess(PointPattern(points).get_window(), 2000, 1, conditioning=False, asPP=True)
pp_csr2 = samples2.realizations[0]
pp_csr2.plot()
q_h = qs.QStatistic(pp_csr2,shape= "hexagon",lh = 0.1)
q_h.plot()
print(q_h.chi2,q_h.df,q_h.chi2_pvalue)
q_r = qs.QStatistic(pp_csr2,shape= "rectangle",nx = 10, ny = 10)
plt.figure(figsize=(12,12))
q_r.plot()
print(q_r.chi2,q_r.df,q_r.chi2_pvalue)
In [134]:
# kernel density estimation
sns.kdeplot(pp10.points.x, pp10.points.y, shade=True, cmap='viridis', cbar = True)
#pp.plot()
Out[134]:
In [119]:
# diagramme de Voronoi
from libpysal.cg.voronoi import voronoi, voronoi_frames, as_dataframes
regions, vertices = voronoi(pp10.points)
#vertices.shape
#plt.plot(vertices[:,0],vertices[:,1])
region_df, point_df = voronoi_frames(np.asarray(pp10.points))
#as_dataframes(regions, vertices, pp10.points)
#points = [(10.2, 5.1), (4.7, 2.2), (5.3, 5.7), (2.7, 5.3)]
#region_df, point_df = voronoi_frames(points)
fig, ax = plt.subplots(1,figsize=(9,9))
region_df.plot(ax=ax, color='blue',edgecolor='black', alpha=0.3)
point_df.plot(ax=ax, color='red')
Out[119]:
In [132]:
region_df['aires'] = 1/region_df.area
f, ax = plt.subplots(1, figsize=(16, 16))
region_df.plot(column='aires', legend = True, ax=ax)
#accidents[::10].plot(ax=ax)
#xmin, xmax = plt.xlim()
#ymin, ymax = plt.ylim()
plt.xlim(290000, 302000)
plt.ylim(5035000, 5046000)
Out[132]:
In [79]:
from libpysal.weights import Queen, Rook, KNN
In [80]:
data = gpd.read_file('donnees/accidents-geobase/accidents_2018_smod13.shp')
f, ax = plt.subplots(1, figsize=(9, 9))
data.plot(column='NACCIDENTS', legend = True, ax=ax)
data.info()
In [81]:
data['DENSITE_ACCIDENTS'] = data['NACCIDENTS']/data.geometry.area
f, ax = plt.subplots(1, figsize=(9, 9))
data.plot(column='DENSITE_ACCIDENTS', legend = True, ax=ax)
Out[81]:
In [82]:
w_rook = Rook.from_dataframe(data)
ax = data.plot(edgecolor='grey', facecolor='w', figsize = (9,9))
f,ax = w_rook.plot(data, ax=ax,
edge_kws=dict(color='r', linestyle=':', linewidth=1),
node_kws=dict(marker=''))
ax.set_axis_off()
w_rook.histogram
Out[82]:
In [83]:
w_queen = Queen.from_dataframe(data)
ax = data.plot(edgecolor='grey', facecolor='w', figsize = (9,9))
f,ax = w_queen.plot(data, ax=ax,
edge_kws=dict(color='r', linestyle=':', linewidth=1),
node_kws=dict(marker=''))
ax.set_axis_off()
w_queen.histogram
Out[83]:
In [84]:
w_knn = KNN.from_dataframe(data, k=4)
ax = data.plot(edgecolor='grey', facecolor='w', figsize = (9,9))
f,ax = w_knn.plot(data, ax=ax,
edge_kws=dict(color='r', linestyle=':', linewidth=1),
node_kws=dict(marker=''))
ax.set_axis_off()
w_knn.histogram
# from libpysal.weights.contiguity import Voronoi as Vornoi_weights
# w = Vornoi_weights(points)
Out[84]:
In [85]:
w_rook.transform = 'r'
data['NACCIDENTS_LAG'] = ps.weights.lag_spatial(w_rook, data.NACCIDENTS)
f, (ax1,ax2) = plt.subplots(2, figsize=(6,6))
data.plot(column='NACCIDENTS', legend = True, ax=ax1)
data.plot(column='NACCIDENTS_LAG', legend = True, ax=ax2)
Out[85]:
In [86]:
# moran plot
f, ax = plt.subplots(1, figsize=(9, 9))
plt.plot(data.NACCIDENTS, data.NACCIDENTS_LAG, '.', color='firebrick')
# dashed vert at mean of the last year's PCI
plt.vlines(data.NACCIDENTS.mean(), data.NACCIDENTS_LAG.min(), data.NACCIDENTS_LAG.max(), linestyle='--')
# dashed horizontal at mean of lagged PCI
plt.hlines(data.NACCIDENTS_LAG.mean(), data.NACCIDENTS.min(), data.NACCIDENTS.max(), linestyle='--')
# red line of best fit using global I as slope
b,a = np.polyfit(data.NACCIDENTS, data.NACCIDENTS_LAG, 1)
plt.plot(data.NACCIDENTS, a + b*data.NACCIDENTS, 'r')
plt.title('Moran Scatterplot')
plt.ylabel('Spatial Lag of NACCIDENTS')
plt.xlabel('NACCIDENTS')
plt.show()
In [87]:
import esda
mi = esda.moran.Moran(data.NACCIDENTS, w_rook)
print(mi.I,mi.p_sim)
In [88]:
geary = esda.geary.Geary(data.NACCIDENTS, w_rook)
print(geary.C, geary.p_sim)
In [89]:
li = esda.moran.Moran_Local(data.NACCIDENTS, w_rook)
data['MORAN'] = li.Is
f, ax = plt.subplots(1, figsize=(9, 9))
data.plot(column='MORAN', legend = True, ax=ax)
li.q
# (if permutations>0) values indicate quandrant location 1 HH, 2 LH, 3 LL, 4 HL
Out[89]:
In [90]:
from matplotlib import colors
sig = li.p_sim < 0.05
hotspots = sig * li.q==1
coldspots = sig * li.q==3
doughnut = sig * li.q==2
diamond = sig * li.q==4
In [91]:
cmap = colors.ListedColormap(['grey', 'blue'])
f, ax = plt.subplots(1, figsize=(9, 9))
data.assign(cl=hotspots*1).plot(column='cl', categorical=True, \
k=2, cmap=cmap, linewidth=0.1, ax=ax, \
edgecolor='black', legend=True)
ax.set_axis_off()
data[hotspots]
Out[91]:
In [135]:
cmap = colors.ListedColormap(['grey', 'blue'])
f, ax = plt.subplots(1, figsize=(9, 9))
data.assign(cl=coldspots*1).plot(column='cl', categorical=True, \
k=2, cmap=cmap, linewidth=0.1, ax=ax, \
edgecolor='black', legend=True)
ax.set_axis_off()
data[coldspots]
Out[135]:
In [136]:
sns.kdeplot(data.NACCIDENTS)
Out[136]:
In [137]:
import sklearn.cluster as skc
In [138]:
clusters = skc.AgglomerativeClustering(n_clusters=3,connectivity=w_rook.sparse).fit(data[['NACCIDENTS']])
data.assign(labels=clusters.labels_).plot('labels', cmap='rainbow')
Out[138]:
In [139]:
clusters = skc.AgglomerativeClustering(n_clusters=3,connectivity=w_knn.sparse).fit(data[['NACCIDENTS']])
data.assign(labels=clusters.labels_).plot('labels', cmap='rainbow')
Out[139]:
In [140]:
clusters = skc.AgglomerativeClustering(n_clusters=3,connectivity=w_rook.sparse).fit(data[['DENSITE_ACCIDENTS']])
data.assign(labels=clusters.labels_).plot('labels', cmap='rainbow')
Out[140]: