TP 2 - ANALYSE BIVARIEE¶

1) Analyse bivariée entre deux variables quantitatives¶

Import des librairies Python :

In [1]:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import geopandas as gpd

Ouverture du jeu de données filosofi :

In [2]:
filosofi = gpd.read_file("carreaux_1km_met.shp")
In [3]:
filosofi.head()#Inspecter les 5 premières lignes du jeu de données
Out[3]:
idcar_1km i_est_1km lcog_geo ind men men_pauv men_1ind men_5ind men_prop men_fmp ... ind_6_10 ind_11_17 ind_18_24 ind_25_39 ind_40_54 ind_55_64 ind_65_79 ind_80p ind_inc geometry
0 CRS3035RES1000mN2029000E4252000 1 2A041 3.0 1.4 0.3 0.5 0.1 0.8 0.1 ... 0.2 0.2 0.3 0.2 0.6 0.2 1.0 0.3 0.0 POLYGON ((1218038.251 6048790.697, 1217949.529...
1 CRS3035RES1000mN2029000E4254000 1 2A041 2.0 0.9 0.2 0.3 0.0 0.5 0.0 ... 0.1 0.1 0.2 0.1 0.4 0.1 0.7 0.3 0.0 POLYGON ((1220028.075 6048968.632, 1219939.34 ...
2 CRS3035RES1000mN2030000E4252000 1 2A041 6.0 2.5 0.3 0.3 0.1 1.6 0.5 ... 0.1 0.5 0.3 0.7 1.5 1.4 1.1 0.3 0.0 POLYGON ((1217949.529 6049793.871, 1217860.808...
3 CRS3035RES1000mN2030000E4253000 1 2A041 16.5 6.9 0.7 1.0 0.3 4.3 1.3 ... 0.3 1.3 1.0 2.0 4.0 3.7 3.0 0.7 0.0 POLYGON ((1218944.435 6049882.842, 1218855.707...
4 CRS3035RES1000mN2030000E4254000 1 2A041 5.0 2.4 0.5 0.9 0.1 1.4 0.1 ... 0.2 0.4 0.5 0.4 1.0 0.4 1.6 0.5 0.0 POLYGON ((1219939.34 6049971.805, 1219850.606 ...

5 rows × 32 columns

On va étudier le lien entre les deux variables choisies suivantes :

  • ind_80p (Population de plus de 80 ans)
  • log_av45 (Nombre de ménages dans un logement d'avant 1945)
In [4]:
!pip install mapclassify

Collecting mapclassify
  Downloading mapclassify-2.10.0-py3-none-any.whl.metadata (3.1 kB)
Requirement already satisfied: networkx>=3.2 in /usr/local/lib/python3.12/dist-packages (from mapclassify) (3.5)
Requirement already satisfied: numpy>=1.26 in /usr/local/lib/python3.12/dist-packages (from mapclassify) (2.0.2)
Requirement already satisfied: pandas>=2.1 in /usr/local/lib/python3.12/dist-packages (from mapclassify) (2.2.2)
Requirement already satisfied: scikit-learn>=1.4 in /usr/local/lib/python3.12/dist-packages (from mapclassify) (1.6.1)
Requirement already satisfied: scipy>=1.12 in /usr/local/lib/python3.12/dist-packages (from mapclassify) (1.16.3)
Requirement already satisfied: python-dateutil>=2.8.2 in /usr/local/lib/python3.12/dist-packages (from pandas>=2.1->mapclassify) (2.9.0.post0)
Requirement already satisfied: pytz>=2020.1 in /usr/local/lib/python3.12/dist-packages (from pandas>=2.1->mapclassify) (2025.2)
Requirement already satisfied: tzdata>=2022.7 in /usr/local/lib/python3.12/dist-packages (from pandas>=2.1->mapclassify) (2025.2)
Requirement already satisfied: joblib>=1.2.0 in /usr/local/lib/python3.12/dist-packages (from scikit-learn>=1.4->mapclassify) (1.5.2)
Requirement already satisfied: threadpoolctl>=3.1.0 in /usr/local/lib/python3.12/dist-packages (from scikit-learn>=1.4->mapclassify) (3.6.0)
Requirement already satisfied: six>=1.5 in /usr/local/lib/python3.12/dist-packages (from python-dateutil>=2.8.2->pandas>=2.1->mapclassify) (1.17.0)
Downloading mapclassify-2.10.0-py3-none-any.whl (882 kB)
   ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 882.2/882.2 kB 13.8 MB/s eta 0:00:00
Installing collected packages: mapclassify
Successfully installed mapclassify-2.10.0
In [5]:
import mapclassify
var_1 = "ind_80p"
var_2 = "log_av45"
filosofi.plot(figsize=(10,10), column=var_1, legend=True, scheme="quantiles", k=5).set_title(var_1)
filosofi.plot(figsize=(10,10), column=var_2, legend=True, scheme="quantiles", k=5).set_title(var_2)
filosofi.plot(figsize=(10,10), kind="scatter", x=var_1, y=var_2)
Out[5]:
<Axes: xlabel='ind_80p', ylabel='log_av45'>
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image

On voit que les distributions des deux variables sont très asymétriques vers la droite. On va donc refaire le scatter plot avec une échelle logarithmique. On peut aussi se questionner sur la présence de valeurs non entières sur les cartes. Elles viennent en fait des valeurs "imputées".

In [6]:
ax = filosofi.plot(figsize=(10,10), kind="scatter", x=var_1, y=var_2, c=np.array(["blue", "green", "red"])[filosofi["i_est_1km"]], s=0.1)
ax.set_xscale("log")
ax.set_yscale("log")
No description has been provided for this image

1.2) Corrélation¶

Un lien statistique semble bien exister entre les deux variables. On va donc calculer la corrélation entre les deux :

In [7]:
from scipy.stats import spearmanr, pearsonr

print("Corrélation de Pearson : ", pearsonr(filosofi[var_1], filosofi[var_2]))
print("Corrélation de Spearman : ", spearmanr(filosofi[var_1], filosofi[var_2]))

Corrélation de Pearson :  PearsonRResult(statistic=np.float64(0.7322642873072113), pvalue=np.float64(0.0))
Corrélation de Spearman :  SignificanceResult(statistic=np.float64(0.7405625506170361), pvalue=np.float64(0.0))

La corrélation est assez forte (R et $\rho$ >0.7), positive et significative (p-value très faible).

Cependant, ces deux variables sont en fait liées toutes les deux au nombre total d'habitants! On va donc diviser chacune par ind.

In [8]:
print("Corrélation de Pearson entre var_1 et ind : ", pearsonr(filosofi[var_1], filosofi["ind"]))
print("Corrélation de Pearson entre var_2 et ind : ", pearsonr(filosofi[var_2], filosofi["ind"]))

Corrélation de Pearson entre var_1 et ind :  PearsonRResult(statistic=np.float64(0.9127798594694922), pvalue=np.float64(0.0))
Corrélation de Pearson entre var_2 et ind :  PearsonRResult(statistic=np.float64(0.774901356046622), pvalue=np.float64(0.0))
In [9]:
filosofi["proportion_habitant_sup_80"] = filosofi["ind_80p"] / filosofi["ind"]
filosofi["proportion_habitant_dans_avant_1945"] = filosofi["log_av45"] / filosofi["ind"]
In [10]:
filosofi.plot(figsize=(10,10), column="proportion_habitant_sup_80", legend=True, scheme="quantiles", k=5).set_title(var_1+" par habitant")
filosofi.plot(figsize=(10,10), column="proportion_habitant_dans_avant_1945", legend=True, scheme="quantiles", k=5).set_title(var_2+" par habitant")
filosofi.plot(figsize=(10,10), kind="scatter", x="proportion_habitant_sup_80", y="proportion_habitant_dans_avant_1945")
Out[10]:
<Axes: xlabel='proportion_habitant_sup_80', ylabel='proportion_habitant_dans_avant_1945'>
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image

Les données sont écrasées vers les faibles valeurs. On va afficher le scatter plot avec une échelle logarithmique pour mieux étaler les valeurs.

In [11]:
ax = filosofi.plot(figsize=(10,10), kind="scatter", x="proportion_habitant_sup_80", y="proportion_habitant_dans_avant_1945")
ax.set_xscale("log")
ax.set_yscale("log")
No description has been provided for this image
In [12]:
print("Corrélation de Pearson : ", pearsonr(filosofi["proportion_habitant_sup_80"], filosofi["proportion_habitant_dans_avant_1945"]))
print("Corrélation de Spearman : ", spearmanr(filosofi["proportion_habitant_sup_80"], filosofi["proportion_habitant_dans_avant_1945"]))

Corrélation de Pearson :  PearsonRResult(statistic=np.float64(0.30676528673616194), pvalue=np.float64(0.0))
Corrélation de Spearman :  SignificanceResult(statistic=np.float64(0.2875337925247301), pvalue=np.float64(0.0))

On constate que la corrélation est bien moins forte ! On a encore une corrélation, encore positive et significative, mais plutôt faible (R=0.3)

1.3) Régression¶

On va essayer de prédire le niveau de vie en fonction du nombre de logements construit avant 1945.

In [13]:
var_3 = "log_av45"
var_4 = "ind_snv"
filosofi.plot(figsize=(10,10), kind="scatter", x=var_3, y=var_4)
Out[13]:
<Axes: xlabel='log_av45', ylabel='ind_snv'>
No description has been provided for this image

Il semble exister un lien statistique entre ces deux variables, et il a l'air à peu près linéaire. Faire une regression semble donc pertinent.

In [14]:
from scipy.stats import linregress

reg = linregress(filosofi[var_3], filosofi[var_4])
print("R²:", (reg.rvalue)**2, "p-value:", reg.pvalue)

R²: 0.7085580778444482 p-value: 0.0

Le R² est assez élevé : le modèle linéaire prédit bien le niveau de vie. La p-value est très faible : on peut rejetter l'indépendance des variables avec une grande certitude.

In [15]:
var_3 = "log_av45"
var_4 = "ind_snv"
ax = filosofi.plot(figsize=(10,10), kind="scatter", x=var_3, y=var_4)
#Affichage de la droite de régression
plt.plot(
    [filosofi[var_3].min(), filosofi[var_3].max()],
    [reg.intercept + reg.slope * filosofi[var_3].min(), reg.intercept + reg.slope * filosofi[var_3].max()],
    color="red"
    )
Out[15]:
[<matplotlib.lines.Line2D at 0x7ff28e11b140>]
No description has been provided for this image

2) Analyse bivariée entre deux variables qualitatives¶

Les distributions de CODE_CS et CODE_US sont-elles indépendantes ?

In [16]:
ocs_ge = gpd.read_file("OCCUPATION_SOL.shp")
In [ ]:
ocs_ge.plot(figsize=(10,10), column="CODE_CS", legend=True)
ocs_ge.plot(figsize=(10,10), column="CODE_US", legend=True)
Out[ ]:
<Axes: >
No description has been provided for this image
No description has been provided for this image

On affiche la table de contingence entre couverture et usage du sol :

In [ ]:
effectifs_observes = pd.crosstab(ocs_ge["CODE_CS"], ocs_ge["CODE_US"])
effectifs_observes
Out[ ]:
CODE_US US1.1 US1.2 US1.3 US1.4 US2 US235 US3 US4.1.1 US4.1.2 US4.1.3 US4.1.4 US4.2 US4.3 US5 US6.1 US6.2 US6.3
CODE_CS
CS1.1.1.1 13113 0 82 8 4024 48 10954 49 133 0 1 19 361 114534 346 89 0
CS1.1.1.2 597 0 8 4 1672 6 5555 33942 145 5 1 13 181 12490 6 0 0
CS1.1.2.1 4329 0 98 9 1051 32 2549 84 1956 0 0 4 98 2687 483 2 0
CS1.2.1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 362
CS1.2.2 12 0 38 11 26 0 157 5 0 0 11 0 73 157 0 0 2696
CS2.1.1.1 259 10879 40 0 483 3 2237 332 0 0 0 0 0 13739 0 0 13656
CS2.1.1.2 0 1701 1 0 3 1 41 1 0 0 0 0 0 104 0 0 141
CS2.1.1.3 0 2101 1 0 16 0 143 9 0 0 0 0 0 733 0 0 351
CS2.1.2 178 4 26 0 18 0 69 54 80 0 0 0 339 93 2 0 2917
CS2.1.3 571 0 0 0 0 0 0 0 0 0 0 0 0 5 0 0 0
CS2.2.1 19263 0 95 7 2111 8 6640 1035 221 10 0 21 598 43590 113 0 4700

On peut directement réaliser le test avec la fonction chi2_contingency de scipy.

In [ ]:
from scipy.stats import chi2_contingency
chi2, p, ddl, expected = chi2_contingency(effectifs_observes)
print("\n--- Résultats du Test du Chi2 ---")
print(f"Statistique du Chi2 (χ²): {chi2:.4f}")
print(f"p-valeur (p): {p:.4f}")
print(f"Degrés de liberté (ddl): {ddl}")

--- Résultats du Test du Chi2 ---
Statistique du Chi2 (χ²): 525035.8906
p-valeur (p): 0.0000
Degrés de liberté (ddl): 160

La p-valeur est très faible (inférieure à $10^{-5}$). On peut donc rejetter l'hypothèse d'indépendance avec une certitude importante.

On peut aussi faire toutes les étapes du Chi2 :

In [ ]:
frequences_observes = effectifs_observes/ocs_ge.shape[0]
frequences_observes
Out[ ]:
CODE_US US1.1 US1.2 US1.3 US1.4 US2 US235 US3 US4.1.1 US4.1.2 US4.1.3 US4.1.4 US4.2 US4.3 US5 US6.1 US6.2 US6.3
CODE_CS
CS1.1.1.1 0.038002 0.000000 0.000238 0.000023 0.011662 0.000139 0.031745 0.000142 0.000385 0.000000 0.000003 0.000055 0.001046 0.331926 0.001003 0.000258 0.000000
CS1.1.1.2 0.001730 0.000000 0.000023 0.000012 0.004846 0.000017 0.016099 0.098366 0.000420 0.000014 0.000003 0.000038 0.000525 0.036197 0.000017 0.000000 0.000000
CS1.1.2.1 0.012546 0.000000 0.000284 0.000026 0.003046 0.000093 0.007387 0.000243 0.005669 0.000000 0.000000 0.000012 0.000284 0.007787 0.001400 0.000006 0.000000
CS1.2.1 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.001049
CS1.2.2 0.000035 0.000000 0.000110 0.000032 0.000075 0.000000 0.000455 0.000014 0.000000 0.000000 0.000032 0.000000 0.000212 0.000455 0.000000 0.000000 0.007813
CS2.1.1.1 0.000751 0.031528 0.000116 0.000000 0.001400 0.000009 0.006483 0.000962 0.000000 0.000000 0.000000 0.000000 0.000000 0.039816 0.000000 0.000000 0.039576
CS2.1.1.2 0.000000 0.004930 0.000003 0.000000 0.000009 0.000003 0.000119 0.000003 0.000000 0.000000 0.000000 0.000000 0.000000 0.000301 0.000000 0.000000 0.000409
CS2.1.1.3 0.000000 0.006089 0.000003 0.000000 0.000046 0.000000 0.000414 0.000026 0.000000 0.000000 0.000000 0.000000 0.000000 0.002124 0.000000 0.000000 0.001017
CS2.1.2 0.000516 0.000012 0.000075 0.000000 0.000052 0.000000 0.000200 0.000156 0.000232 0.000000 0.000000 0.000000 0.000982 0.000270 0.000006 0.000000 0.008454
CS2.1.3 0.001655 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000014 0.000000 0.000000 0.000000
CS2.2.1 0.055825 0.000000 0.000275 0.000020 0.006118 0.000023 0.019243 0.002999 0.000640 0.000029 0.000000 0.000061 0.001733 0.126326 0.000327 0.000000 0.013621
In [ ]:
frequences_usages = frequences_observes.sum(axis=0)
display(frequences_usages)

frequences_couvertures = frequences_observes.sum(axis=1)
display(frequences_couvertures)
0
CODE_US
US1.1 0.111059
US1.2 0.042558
US1.3 0.001127
US1.4 0.000113
US2 0.027253
US235 0.000284
US3 0.082145
US4.1.1 0.102913
US4.1.2 0.007347
US4.1.3 0.000043
US4.1.4 0.000038
US4.2 0.000165
US4.3 0.004782
US5 0.545217
US6.1 0.002753
US6.2 0.000264
US6.3 0.071938

0
CODE_CS
CS1.1.1.1 0.416627
CS1.1.1.2 0.158306
CS1.1.2.1 0.038782
CS1.2.1 0.001049
CS1.2.2 0.009233
CS2.1.1.1 0.120640
CS2.1.1.2 0.005776
CS2.1.1.3 0.009720
CS2.1.2 0.010955
CS2.1.3 0.001669
CS2.2.1 0.227242

In [ ]:
frequences_theoriques = pd.DataFrame(
    frequences_couvertures.to_numpy().reshape(-1,1).dot(frequences_usages.to_numpy().reshape(1,-1)),
    index = frequences_couvertures.index,
    columns = frequences_usages.index
)
frequences_theoriques
Out[ ]:
CODE_US US1.1 US1.2 US1.3 US1.4 US2 US235 US3 US4.1.1 US4.1.2 US4.1.3 US4.1.4 US4.2 US4.3 US5 US6.1 US6.2 US6.3
CODE_CS
CS1.1.1.1 0.046270 0.017731 0.000470 4.708895e-05 0.011354 1.183261e-04 0.034224 0.042876 0.003061 1.811113e-05 1.569632e-05 6.882231e-05 0.001992 0.227152 0.001147 1.098742e-04 0.029972
CS1.1.1.2 0.017581 0.006737 0.000178 1.789243e-05 0.004314 4.496047e-05 0.013004 0.016292 0.001163 6.881704e-06 5.964143e-06 2.615048e-05 0.000757 0.086311 0.000436 4.174900e-05 0.011388
CS1.1.2.1 0.004307 0.001650 0.000044 4.383277e-06 0.001057 1.101439e-05 0.003186 0.003991 0.000285 1.685876e-06 1.461092e-06 6.406328e-06 0.000185 0.021144 0.000107 1.022765e-05 0.002790
CS1.2.1 0.000117 0.000045 0.000001 1.185732e-07 0.000029 2.979531e-07 0.000086 0.000108 0.000008 4.560507e-08 3.952439e-08 1.732993e-07 0.000005 0.000572 0.000003 2.766707e-07 0.000075
CS1.2.2 0.001025 0.000393 0.000010 1.043575e-06 0.000252 2.622317e-06 0.000758 0.000950 0.000068 4.013750e-07 3.478583e-07 1.525225e-06 0.000044 0.005034 0.000025 2.435008e-06 0.000664
CS2.1.1.1 0.013398 0.005134 0.000136 1.363526e-05 0.003288 3.426296e-05 0.009910 0.012415 0.000886 5.244331e-06 4.545087e-06 1.992846e-05 0.000577 0.065775 0.000332 3.181561e-05 0.008679
CS2.1.1.2 0.000641 0.000246 0.000007 6.528076e-07 0.000157 1.640388e-06 0.000474 0.000594 0.000042 2.510798e-07 2.176025e-07 9.541034e-07 0.000028 0.003149 0.000016 1.523218e-06 0.000416
CS2.1.1.3 0.001080 0.000414 0.000011 1.098603e-06 0.000265 2.760593e-06 0.000798 0.001000 0.000071 4.225398e-07 3.662011e-07 1.605651e-06 0.000046 0.005300 0.000027 2.563408e-06 0.000699
CS2.1.2 0.001217 0.000466 0.000012 1.238140e-06 0.000299 3.111223e-06 0.000900 0.001127 0.000080 4.762076e-07 4.127133e-07 1.809589e-06 0.000052 0.005973 0.000030 2.888993e-06 0.000788
CS2.1.3 0.000185 0.000071 0.000002 1.886689e-07 0.000045 4.740911e-07 0.000137 0.000172 0.000012 7.256497e-08 6.288964e-08 2.757469e-07 0.000008 0.000910 0.000005 4.402275e-07 0.000120
CS2.2.1 0.025237 0.009671 0.000256 2.568387e-05 0.006193 6.453895e-05 0.018667 0.023386 0.001669 9.878411e-06 8.561289e-06 3.753796e-05 0.001087 0.123896 0.000626 5.992902e-05 0.016347
In [ ]:
effectifs_theoriques = frequences_theoriques * ocs_ge.shape[0]
effectifs_theoriques
Out[ ]:
CODE_US US1.1 US1.2 US1.3 US1.4 US2 US235 US3 US4.1.1 US4.1.2 US4.1.3 US4.1.4 US4.2 US4.3 US5 US6.1 US6.2 US6.3
CODE_CS
CS1.1.1.1 15965.991445 6118.171921 162.068020 16.248465 3917.963143 40.829476 11809.300859 14794.852101 1056.150209 6.249410 5.416155 23.747756 687.435047 78380.927470 395.795936 37.913084 10341.939503
CS1.1.1.2 6066.612521 2324.727438 61.581135 6.173944 1488.712075 15.514014 4487.190959 5621.613623 401.306371 2.374594 2.057981 9.023457 261.205330 29782.473432 150.390948 14.405870 3929.636309
CS1.1.2.1 1486.195126 569.510345 15.086110 1.512489 364.703798 3.800614 1099.269371 1377.179561 98.311796 0.581727 0.504163 2.210561 63.989926 7296.092622 36.842685 3.529141 962.679965
CS1.2.1 40.203455 15.405974 0.408098 0.040915 9.865698 0.102811 29.736625 37.254446 2.659458 0.015736 0.013638 0.059798 1.731008 197.368520 0.996641 0.095468 26.041709
CS1.2.2 353.834828 135.589595 3.591716 0.360095 86.829047 0.904854 261.715156 327.880293 23.406171 0.138498 0.120032 0.526293 15.234786 1737.061059 8.771543 0.840222 229.195813
CS2.1.1.1 4623.175214 1771.601900 46.929053 4.704969 1134.500801 11.822743 3419.547556 4284.055504 305.823004 1.809604 1.568323 6.876494 199.056393 22696.289319 114.608226 10.978262 2994.652636
CS2.1.1.2 221.341121 84.817973 2.246795 0.225257 54.315847 0.566031 163.715727 205.105281 14.641713 0.086637 0.075086 0.329222 9.530109 1086.617292 5.487033 0.525600 143.373275
CS2.1.1.3 372.492785 142.739329 3.781110 0.379083 91.407603 0.952568 275.515578 345.169649 24.640395 0.145801 0.126361 0.554044 16.038127 1828.657499 9.234073 0.884527 241.281468
CS2.1.2 419.804034 160.869011 4.261358 0.427231 103.017513 1.073556 310.509507 389.010517 27.770034 0.164320 0.142410 0.624415 18.075170 2060.919901 10.406916 0.996873 271.927236
CS2.1.3 63.970138 24.513373 0.649350 0.065102 15.697907 0.163589 47.315734 59.277793 4.231624 0.025039 0.021701 0.095149 2.754312 314.044937 1.585816 0.151904 41.436531
CS2.2.1 8708.379332 3337.053142 88.397254 8.862450 2136.986568 22.269745 6441.182928 8069.601233 576.059225 3.408634 2.954150 12.952811 374.949791 42751.547950 215.880183 20.679049 5640.835556
In [ ]:
test_chi2 = ((effectifs_observes - effectifs_theoriques)**2/effectifs_theoriques).sum().sum()
test_chi2
Out[ ]:
np.float64(525035.8906068673)
In [ ]:
from scipy.stats import chi2

degre_de_liberte = (len(frequences_usages)-1) * (len(frequences_couvertures)-1)
p_value = 0.01

valeur_critique = chi2.ppf(1 - p_value, degre_de_liberte)

print(f"La valeur critique du chi2 pour un degré de liberté={degre_de_liberte} et une p-value={p_value} est : {valeur_critique}")

La valeur critique du chi2 pour un degré de liberté=160 et une p-value=0.01 est : 204.5300945903855

La valeur obtenue par le test du Chi² est bien supérieure à la valeur critique. On peut donc rejetter l'indépendance des deux variables, avec moins d'1% de chance de se tromper.

Analyse bivariée entre une variable qualitative et une variable quantitative¶

On va d'abord restreindre le fichier filosofi au département de l'Ain. On va utiliser le code des communes pour faire cette selection facilement

In [17]:
filosofi_ain = filosofi[filosofi["lcog_geo"].str.startswith("01")]
#Ce qu'il y a entre les crochets est une pd.Series de Booleen "True" quand le carreaux est dans l'Ain et False sinon.
#La syntaxe df[booleen_serie] permet de renvoyer le sous ensemble du dataframe initial avec tous les indices pour lesquels booleen_serie est "True"
In [18]:
filosofi_ain.plot(column="ind", legend=True)
Out[18]:
<Axes: >
No description has been provided for this image

On va réaliser l'intersection entre les deux couches

In [19]:
intersection = ocs_ge.overlay(filosofi_ain,
                              how="intersection",
                              keep_geom_type=True
                              )

Chaque polygone de cette intersection a un usage bien défini. On va désagréger la population à l'échelle de ces unités cartographiques en faisant l'hypothèse que la population est répartie de manière homogène dans les carreaux de 1 km.

In [20]:
intersection["ind_estimee"] = intersection["ind"] * (intersection.area/1000000)
In [21]:
intersection.head()
Out[21]:
ID CODE_CS CODE_US MILLESIME SOURCE OSSATURE ID_ORIGINE CODE_OR idcar_1km i_est_1km ... ind_25_39 ind_40_54 ind_55_64 ind_65_79 ind_80p ind_inc proportion_habitant_sup_80 proportion_habitant_dans_avant_1945 geometry ind_estimee
0 OCSGE0000000010019471939 CS1.1.1.1 US5 2021 calcul 0 NC NC CRS3035RES1000mN2580000E4012000 0 ... 6.0 31.0 7.0 11.0 2.1 0.0 0.023204 0.000000 POLYGON ((930611.6 6577732.26, 930609.6 657773... 0.029604
1 OCSGE0000000010019471929 CS1.1.1.1 US5 2021 calcul 0 NC NC CRS3035RES1000mN2580000E4011000 1 ... 41.8 66.0 34.5 47.7 13.9 0.0 0.047039 0.019966 POLYGON ((930397.1 6577673.9, 930388.11 657768... 0.128211
2 OCSGE0000000010019471930 CS1.1.1.1 US5 2021 calcul 0 NC NC CRS3035RES1000mN2580000E4011000 1 ... 41.8 66.0 34.5 47.7 13.9 0.0 0.047039 0.019966 POLYGON ((929937.96 6577691.9, 929926.31 65777... 0.059096
3 OCSGE0000000010019471931 CS1.1.1.1 US5 2021 calcul 0 NC NC CRS3035RES1000mN2580000E4012000 0 ... 6.0 31.0 7.0 11.0 2.1 0.0 0.023204 0.000000 POLYGON ((930752.74 6577721, 930766.5 6577715.... 0.021780
4 OCSGE0000000010019471932 CS1.1.1.1 US5 2021 calcul 0 NC NC CRS3035RES1000mN2580000E4011000 1 ... 41.8 66.0 34.5 47.7 13.9 0.0 0.047039 0.019966 POLYGON ((930382.8 6577714.3, 930394.16 657772... 0.062773

5 rows × 43 columns

In [22]:
intersection["ind_estimee"].head()
Out[22]:
ind_estimee
0 0.029604
1 0.128211
2 0.059096
3 0.021780
4 0.062773

Certains morceaux n'ont pas été intersectés. Ils ont pour l'instant une valeur "ind_estimee": NaN. Pour notre analyse, on va considérer qu'ils ont une population nulle.

In [23]:
intersection["ind_estimee"] = intersection["ind_estimee"].fillna(0)

On va regrouper au niveau des unités cartographiques initiales de l'OCS GE. Pour cela, on va faire la somme des populations estimées par ID de l'OCS GE.

In [24]:
pop_par_ocsge = intersection.groupby("ID")["ind_estimee"].sum()
In [25]:
ocs_ge = ocs_ge.set_index("ID")
ocs_ge["ind_estimee"] = pop_par_ocsge
In [26]:
ocs_ge.plot(kind="hist", column="ind_estimee", bins=100)
Out[26]:
<Axes: ylabel='Frequency'>
No description has been provided for this image

L'histogramme montre qu'une très grande majorité des polygones a une population très faible, et que la distribution est très asymétrique. On va donc faire une tranformation logarithmique

In [31]:
ocs_ge["log_ind_estimee"] = np.log10(ocs_ge["ind_estimee"])
ocs_ge.plot(kind="hist", column="log_ind_estimee", bins=100)
Out[31]:
<Axes: ylabel='Frequency'>
No description has been provided for this image
In [32]:
ocs_ge.boxplot(column="log_ind_estimee", by="CODE_US", figsize=(10,10), rot=90)
Out[32]:
<Axes: title={'center': 'log_ind_estimee'}, xlabel='CODE_US'>
No description has been provided for this image

On constate que les distributions des populations sont différentes selon les usages : il semble exister un lien entre les deux. Pour autant, le modèle qu'on a utilisé pour la désagrégation de la population ne semble pas très pertinent : l'usage avec la plus haute population médiane est l'usage "Navigable" ! D'autre part les unités spatiales ont des tailles très variées, et comparer directement leur population n'est pas si pertinent.

In [33]:
ocs_ge["densite"] = ocs_ge["ind_estimee"]/ocs_ge.area *1000000 #(habitants/km²)
ocs_ge.plot(kind="hist", column="densite", bins=100)
ocs_ge.boxplot(column="densite", by="CODE_US", figsize=(10,10), rot=90)
Out[33]:
<Axes: title={'center': 'densite'}, xlabel='CODE_US'>
No description has been provided for this image
No description has been provided for this image
In [34]:
ocs_ge["log_densite"] = np.log10(ocs_ge["densite"])
ocs_ge.plot(kind="hist", column="log_densite", bins=100)
ocs_ge.boxplot(column="log_densite", by="CODE_US", figsize=(10,10), rot=90)
Out[34]:
<Axes: title={'center': 'log_densite'}, xlabel='CODE_US'>
No description has been provided for this image
No description has been provided for this image

On va essayer de faire l'inverse : plutôt que de calculer la population des polygones de l'OCS GE, on va chercher l'usage majoritaire des carreaux.

In [ ]:
intersection["intersected_area"] = intersection.area
grouped = intersection.groupby(by=["idcar_1km", "CODE_US"])
In [ ]:
area_per_US = grouped.aggregate({"intersected_area": "sum"})
In [ ]:
area_per_US
Out[ ]:
intersected_area
idcar_1km CODE_US
CRS3035RES1000mN2510000E3979000 US1.1 199509.974876
US1.2 247506.319747
US3 148022.904047
US4.1.1 8643.642790
US5 26332.931045
... ... ...
CRS3035RES1000mN2613000E3932000 US1.2 72979.564333
US3 6874.732000
US4.1.1 5182.361158
US5 2594.768041
US6.3 94273.165015

25388 rows × 1 columns

In [ ]:
 def main_us(groupe):
    #On définit une fonction d'aggrégation qu'on va appliquer à chaque groupe de carreaux, et qui va retourner l'usage avec la plus grande surface
    return groupe["CODE_US"].iloc[np.argmax(groupe["intersected_area"])]
In [ ]:
usage_majoritaire = area_per_US.reset_index().groupby("idcar_1km").apply(main_us, include_groups=False)
In [ ]:
filosofi_ain = filosofi_ain.set_index("idcar_1km")
filosofi_ain["Usage majoritaire"] = usage_majoritaire
In [ ]:
filosofi_ain.plot(column="Usage majoritaire", legend=True)
Out[ ]:
<Axes: >
No description has been provided for this image
In [ ]:
filosofi_ain.boxplot(column="ind", by="Usage majoritaire", figsize=(10,10), rot=90)
Out[ ]:
<Axes: title={'center': 'ind'}, xlabel='Usage majoritaire'>
No description has been provided for this image

Ça semble plus crédible

La carte permet de constater une limite à cette approche : les trous correspondent à "des carreaux inhabités", mais qui ont pourtant un usage.

In [ ]: