Inférence d’un réseau entre stations météorologiques

Dans ce projet, il est proposé d’analyser un jeu de données issues de mesures prises en différentes stations météorologiques durant l’année 2016 dans le nord-ouest de la France. Le but de projet est l’inférence d’un graphe de corrélation partielle entre certaines variables météorologiques à l’aide des outils d’estimation en grande dimension dans des modèles graphiques.

Données disponibles sur les stations de mesure

Sur l’année 2016, on observe des données d’apprentissage sur \(p = 262\) stations météorologiques dont on dispose des coordonnées spatiales (latitude et longitude).

Pour chaque station \(1 \leq j \leq p\), on dispose des mesures suivantes :

Variables explicatives : mesure de \(q = 6\) variables \(X_{ijt} = (X_{ijt}^{(k)})_{1 \leq k \leq q} \in \mathbb{R}^{q}\) pour la station \(j\), le jour \(i\) et l’heure \(t \in \{0,\ldots,23 \}\) (variable ordonnée). Les mesures sont

‘ff’ : Vitesse du vent (en m/s)
‘t’ : Température (en Kelvin)
‘td’ : Point de rosée (en Kelvin)
‘hu’ : Humidité (en pourcentage)
‘dd’ : Direction du vent (en degrés)
‘precip’ : Cumul de pluie sur une heure (en ml)

Travail à effectuer par groupe.

pour une ou plusieurs des 6 variables ci-dessus, vous devez créer une matrice \(X\) de taille \(n \times p\) avec \(n=366\) (nombre de jours en 2016) dont chaque entrée \(X_{ij}\) est une statistique (par exemple la moyenne sur 24 heures) issue des observations journalière d’une de ces variables le jour \(i\) et pour la station \(j\). On veillera à retirer une tendance annuelle (moyenne journalière sur l’ensemble des stations) liée à la mesure choisie pour construire les \(X_{ij}\).
inférer à partir de la matrice \(X\), un graphe de corrélation partielle entre les stations et le représenter sur une carte

On présente ci-dessous un exemple de lecture de ces données, et la représentation géographique des positions des stations météorologiques avec des liens qui correpondent aux \(k=2\) plus proches voisins de chaque station.

Dans un premier temps, vous pourrez chercher à construire un graphe basé sur la notion usuelle de corrélation, en indiquant par des arêtes entre deux stations les corrélations non-nulles (estimées à l’aide d’un test d’hypothèses multiples).

# Lecture des données
X_2016 = read.csv("X_2016_final.csv")
dim(X_2016)

## [1] 92209   146

names(X_2016)

##   [1] "ff_0"       "ff_1"       "ff_10"      "ff_11"      "ff_12"     
##   [6] "ff_13"      "ff_14"      "ff_15"      "ff_16"      "ff_17"     
##  [11] "ff_18"      "ff_19"      "ff_2"       "ff_20"      "ff_21"     
##  [16] "ff_22"      "ff_23"      "ff_3"       "ff_4"       "ff_5"      
##  [21] "ff_6"       "ff_7"       "ff_8"       "ff_9"       "t_0"       
##  [26] "t_1"        "t_10"       "t_11"       "t_12"       "t_13"      
##  [31] "t_14"       "t_15"       "t_16"       "t_17"       "t_18"      
##  [36] "t_19"       "t_2"        "t_20"       "t_21"       "t_22"      
##  [41] "t_23"       "t_3"        "t_4"        "t_5"        "t_6"       
##  [46] "t_7"        "t_8"        "t_9"        "td_0"       "td_1"      
##  [51] "td_10"      "td_11"      "td_12"      "td_13"      "td_14"     
##  [56] "td_15"      "td_16"      "td_17"      "td_18"      "td_19"     
##  [61] "td_2"       "td_20"      "td_21"      "td_22"      "td_23"     
##  [66] "td_3"       "td_4"       "td_5"       "td_6"       "td_7"      
##  [71] "td_8"       "td_9"       "hu_0"       "hu_1"       "hu_10"     
##  [76] "hu_11"      "hu_12"      "hu_13"      "hu_14"      "hu_15"     
##  [81] "hu_16"      "hu_17"      "hu_18"      "hu_19"      "hu_2"      
##  [86] "hu_20"      "hu_21"      "hu_22"      "hu_23"      "hu_3"      
##  [91] "hu_4"       "hu_5"       "hu_6"       "hu_7"       "hu_8"      
##  [96] "hu_9"       "dd_0"       "dd_1"       "dd_10"      "dd_11"     
## [101] "dd_12"      "dd_13"      "dd_14"      "dd_15"      "dd_16"     
## [106] "dd_17"      "dd_18"      "dd_19"      "dd_2"       "dd_20"     
## [111] "dd_21"      "dd_22"      "dd_23"      "dd_3"       "dd_4"      
## [116] "dd_5"       "dd_6"       "dd_7"       "dd_8"       "dd_9"      
## [121] "precip_0"   "precip_1"   "precip_10"  "precip_11"  "precip_12" 
## [126] "precip_13"  "precip_14"  "precip_15"  "precip_16"  "precip_17" 
## [131] "precip_18"  "precip_19"  "precip_2"   "precip_20"  "precip_21" 
## [136] "precip_22"  "precip_23"  "precip_3"   "precip_4"   "precip_5"  
## [141] "precip_6"   "precip_7"   "precip_8"   "precip_9"   "number_sta"
## [146] "day"

# Nombre de stations avec observations en 2016
names_stations = unique(X_2016$number_sta)
length(names_stations)

## [1] 262

# Nombre de jours
length(unique(X_2016$day))

## [1] 366

# Importation des coordonnées géographiques des stations de mesures
library(readr)
stations_coordinates_all <- read_csv("stations_coordinates.csv")
names(stations_coordinates_all)

## [1] "number_sta" "lat"        "lon"        "height_sta"

# Intersection avec le nombre de stations avec observations en 2016
length(stations_coordinates_all$number_sta)

## [1] 325

pos = is.element(stations_coordinates_all$number_sta, names_stations)
sum(pos)

## [1] 262

stations_coordinates = stations_coordinates_all[pos,] 
dim(stations_coordinates)

## [1] 262   4

library(leaflet)
library(igraph)

# Recherche de kNN
library(FastKNN)
dist_mat <- as.matrix(dist(stations_coordinates[,c(2,3)]))
k <- 2 # Nbre de plus proches voisins
nrst <- lapply(1:nrow(dist_mat), function(i) k.nearest.neighbors(i, dist_mat, k = k))
# print(nrst)

# Construction d'une matrice adjacence A symmétrique
A <- matrix(0, nrow = dim(dist_mat), ncol=dim(dist_mat)) 
for(i in 1:length(nrst)) for(j in nrst[[i]]) A[i,j] = 1 
for(i in 1:length(nrst)) for(j in nrst[[i]]) A[j,i] = 1 

# Obtention du data.frame d'information sur les somments et les arretes du graphe
rownames(A) <- colnames(A) <- as.matrix(stations_coordinates[,1])
g <- graph_from_adjacency_matrix(t(A), weighted=TRUE) %>%
  set_vertex_attr("longitude", value = as.matrix(stations_coordinates["lon"])) %>%
  set_vertex_attr("latitude", value = as.matrix(stations_coordinates["lat"]))
gg <- get.data.frame(g, "both")

# Collection de tous les coornodonees à tracer 
vert <- gg$vertices
rownames(vert) <- vert$name

# Importation de la collection des coornodonees dans le systeme de SpatialLines
library(sp)
coordinates(vert) <- ~longitude+latitude

# Identifier les droites à tracer et les transferer en classe SpatialLines
line <- function(i){
  return(as(rbind(vert[vert$name == gg$edges[i, "from"], ], 
                    vert[vert$name == gg$edges[i, "to"], ]), "SpatialLines")) # arrets
}

edges <- lapply(1:nrow(gg$edges), line)
for (i in seq_along(edges)) {
  edges[[i]] <- spChFIDs(edges[[i]], as.character(i))
}
edges <- do.call(rbind, edges)

noloop <- gg$edges$from != gg$edges$to

# Visualisation
N <- dim(A)[1]
leaflet(vert[1:N,]) %>% addTiles() %>% 
  addCircleMarkers(data = vert[1:N,])%>% 
  addPolylines(data = edges[1:nrow(gg$edges),radius=0.05], 
               weight = 5*gg$edges$weight)

# Construction de la matrice des données
# Nombre de stations avec observations en 2016
names_stations = unique(X_2016$number_sta)
p = length(names_stations)
print(p)

## [1] 262

# Nombre de jours
n = length(unique(X_2016$day))
print(n)

## [1] 366

liste_jours = sort(unique(X_2016$day))
head(X_2016[,c("number_sta","day")])

# Choix d'une variable
k = 2 # Température (en Kelvin)
nb_heures = 24

# Matrice des données pour cette variable
X = matrix(0,n,p)

for (i in 1:n) {
  # print(i)
  for (j in 1:p) {
    pos = (X_2016[,"day"] == liste_jours[i]) & (X_2016[,"number_sta"] == names_stations[j])
    X[i,j] = mean(as.numeric(X_2016[pos,1+(nb_heures*(k-1)):(nb_heures*k-1)]),na.rm=TRUE)
  }
}

# Visualisation
matplot(1:n,X[,1:10]-273.15,type="l")

# Visulisation
matplot(1:n,X[,100:110]-273.15,type="l")

Inférence d’un réseau entre stations météorologiques

Université de Bordeaux - M2 MAS - UE : Statistique en grande dimension

Année 2024-2025

Données disponibles sur les stations de mesure