Étude données météorologiques à Besançon d’août à octobre 2018
L.-C. Canon
27 août 2019
On a les données de 2 stations: Besançon (relevés amateurs) et Besançon-Thyse (Météo France) séparées de 2 km et toutes les deux à plus de 10 km du point d’enregistrement. En plus, meteociel.fr permet aussi d’avoir accès aux données historique.
infoclimat <- read_rds("infoclimat.rds")
infoclimat_t <- read_rds("infoclimat_t.rds")
meteociel <- read_rds("meteociel.rds")On s’intéresse à analyser la cohérence et la précision des relevés suivants: température, force du vent, direction du vent, précipitation. De plus, on inclut la pression et l’humidité qui peuvent être moins sensible à la localisation. On ignore donc la biométéo, les radiation, le point de rosée et la visibilité.
Relevés amateurs (chaque 10 minutes)
Observation chronologique
On va commencer par afficher les relevés sur une période de 3 jours au cœur de la période étudié (à partir du 1er septembre 2018) avec les données de la station amateur.
infoclimat %>%
gather(mesure, valeur, Température, pluie_heure, Humidité, vent_heure,
direction, Pression) %>%
filter(Heure > make_date(2018, 9, 1)) %>%
filter(Heure < make_date(2018, 9, 4)) %>%
ggplot(aes(x = Heure, y = valeur)) +
facet_wrap(~ mesure, scales = "free_y") +
geom_point(size = 1)## Warning: Removed 360 rows containing missing values (geom_point).
On observe des phénomènes journaliers sur la direction des vents. D’autre part, la direction du vent prend une certaine valeur lorsque sa force est nulle. Finalement, la force des vents semble prendre des valeurs discrètes.
Vent
directions <- sort(unique(infoclimat$direction))
table(tail(directions, -1) - head(directions, -1))##
## 1 2 3 4 5
## 329 5 1 2 1speed <- sort(unique(infoclimat$vent_heure))
table(tail(speed, -1) - head(speed, -1))##
## 1 2
## 7 11La direction semble précise à un degré près, tandis que la force semble l’être à 1 ou 2 km/h.
On étudie la relation entre heure de la journée et direction du vent (lors qu’il est de force non-nulle) sur la période de l’étude (du 22/8 au 30/10).
infoclimat %>%
filter(vent_heure != 0) %>%
filter(Heure > make_date(2018, 8, 22)) %>%
filter(Heure < make_date(2018, 10, 31)) %>%
mutate(Heure_jour = hms::hms(hours = hour(Heure), minutes = minute(Heure))) %>%
ggplot(aes(x = Heure_jour, y = direction, color = Heure)) +
geom_point(size = 0.3) +
scale_y_continuous(breaks = (0:4) * 90) +
scale_x_time(breaks = hms::hms(hours = (0:4) * 6))
Si le vent vient parfois du sud-ouest/ouest entre 12h et 20h, il est principalement du nord-est durant la nuit.
On peut confirmer ces observations avec une rose des vents (entre 22h et 5h).
# compass rose
infoclimat %>%
filter(vent_heure != 0) %>%
filter(Heure > make_date(2018, 8, 22)) %>%
filter(Heure < make_date(2018, 10, 31)) %>%
filter(hour(Heure) < 5 | 22 < hour(Heure)) %>%
ggplot(aes(x = direction,
fill = cut_number(pmax(0, jitter(vent_heure)), 5))) +
geom_histogram(breaks = (0:16) * (360 / 16)) +
# geom_density(aes(y = stat(count)), position = "stack") +
# expand_limits(x = 0) +
scale_x_continuous(breaks = (0:3) * 90) +
scale_fill_viridis_d("Wind speed") +
coord_polar()
On constate une dominance très forte des vents du nord-est (toute vitesse confondue).
Pluie
La dernière mesure d’intérêt est la pluie et il y a deux types de relevés pour la pluie. Commençons par trouver une période pluvieuse.
infoclimat %>%
mutate(day = date(Heure)) %>%
group_by(day) %>%
summarise(pluie = sum(pluie_heure, na.rm = TRUE)) %>%
arrange(desc(pluie)) %>%
head()## # A tibble: 6 x 2
## day pluie
## <date> <dbl>
## 1 2018-09-06 10.6
## 2 2018-09-23 7.6
## 3 2018-10-30 7.2
## 4 2018-10-29 6.6
## 5 2018-08-13 6
## 6 2018-08-25 5.6Le 6 septembre est le jour le plus pluvieux. Observons les deux relevé de pluie.
infoclimat %>%
gather(mesure, valeur, pluie_heure, pluie_inst) %>%
filter(Heure > make_date(2018, 9, 6)) %>%
filter(Heure < make_date(2018, 9, 7)) %>%
ggplot(aes(x = Heure, y = valeur)) +
facet_wrap(~ mesure, scales = "free_y", nrow = 2) +
geom_point()## Warning: Removed 251 rows containing missing values (geom_point).
Les deux types de relevés pour la pluie semblent cohérents.
(pluie_data <- infoclimat %>%
mutate(Heure_heure = ceiling_date(Heure, "1 hour")) %>%
group_by(Heure_heure) %>%
summarise(pluie_heure = max(pluie_heure, na.rm = TRUE),
pluie_inst_agg = sum(pluie_inst, na.rm = TRUE) / 6) %>%
arrange(pluie_heure, desc(pluie_inst_agg)))## # A tibble: 2,209 x 3
## Heure_heure pluie_heure pluie_inst_agg
## <dttm> <dbl> <dbl>
## 1 2018-08-25 21:00:00 -Inf 2.57
## 2 2018-08-09 02:00:00 -Inf 1.8
## 3 2018-08-09 16:00:00 -Inf 0.667
## 4 2018-08-09 13:00:00 -Inf 0.20
## 5 2018-08-01 09:00:00 -Inf 0
## 6 2018-08-01 10:00:00 -Inf 0
## 7 2018-08-01 11:00:00 -Inf 0
## 8 2018-08-01 12:00:00 -Inf 0
## 9 2018-08-01 13:00:00 -Inf 0
## 10 2018-08-01 14:00:00 -Inf 0
## # … with 2,199 more rowspluie_data %>%
filter(is.finite(pluie_heure)) %>%
summarise(corr = cor(pluie_heure, pluie_inst_agg))## # A tibble: 1 x 1
## corr
## <dbl>
## 1 0.838pluie_data %>%
ggplot(aes(x = pluie_heure, y = pluie_inst_agg, color = Heure_heure)) +
geom_point() +
stat_smooth(method = "lm")## Warning: Removed 71 rows containing non-finite values (stat_smooth).
Les deux types de relevés sont corrélés malgré quelques valeurs manquantes pour la pluie par heure. La qualité de la données pluviométriques est donc discutable mais on peut au moins espérer avoir une indication qualitative sur les jours pluvieux à défaut d’avoir des données quantitatives fiables.
Vérifions la cohérence entre l’humidité et la pluie:
infoclimat %>%
mutate(Heure_heure = ceiling_date(Heure, "1 hour")) %>%
group_by(Heure_heure) %>%
summarise(Humidité = max(Humidité, na.rm = TRUE),
pluie_heure = max(pluie_heure, na.rm = TRUE),
pluie_inst_agg = sum(pluie_inst, na.rm = TRUE) / 6) %>%
filter(pluie_heure > 0 | pluie_inst_agg > 0) %>%
ggplot(aes(x = Humidité)) +
stat_ecdf()
Les trois-quarts des journées pluvieuse présentent une humidité maximale supérieure à 90%.
Relevés Météo France (chaque heure)
Observation chronologique
On se focalise d’abord sur la même période de 3 jours.
infoclimat_t %>%
gather(mesure, valeur, Température, Pluie, Humidité, direction, Pression,
vent_heure) %>%
filter(Heure > make_date(2018, 9, 1)) %>%
filter(Heure < make_date(2018, 9, 4)) %>%
ggplot(aes(x = Heure, y = valeur)) +
facet_wrap(~ mesure, scales = "free_y") +
geom_rect(aes(xmin = coucher, ymin = -Inf, xmax = lever, ymax = Inf)) +
geom_point(size = 1)
Vent
On observe les mêmes phénomènes avec cependant aucune précipitation et des valeurs de forces de vent encore plus discrètes.
directions <- sort(unique(infoclimat_t$direction))
table(tail(directions, -1) - head(directions, -1))##
## 10
## 36speed <- sort(unique(infoclimat_t$vent_heure))
table(tail(speed, -1) - head(speed, -1))##
## 2 3 4
## 2 2 3La direction est précise à 10 degrés près, tandis que la force semble l’être à 2, voire 4 km/h.
On étudie à nouveau la relation entre heure de la journée et direction du vent.
infoclimat_t %>%
filter(vent_heure != 0) %>%
filter(Heure > make_date(2018, 8, 22)) %>%
filter(Heure < make_date(2018, 10, 31)) %>%
mutate(Heure_jour = hms::hms(hours = hour(Heure), minutes = minute(Heure))) %>%
ggplot(aes(x = Heure_jour, y = direction, color = Heure)) +
geom_jitter(size = 0.3) +
scale_y_continuous(breaks = (0:4) * 90) +
scale_x_time(breaks = hms::hms(hours = (0:4) * 6))
Les résultats sont similaires. On refait une rose des vents.
# compass rose
infoclimat_t %>%
filter(vent_heure != 0) %>%
filter(Heure > make_date(2018, 8, 22)) %>%
filter(Heure < make_date(2018, 10, 31)) %>%
filter(hour(Heure) < 5 | 22 < hour(Heure)) %>%
ggplot(aes(x = direction,
fill = cut_number(pmax(0, jitter(vent_heure)), 5))) +
geom_histogram(breaks = (0:16) * (360 / 16)) +
# geom_density(aes(y = stat(count)), position = "stack") +
# expand_limits(x = 0) +
scale_x_continuous(breaks = (0:3) * 90) +
scale_fill_viridis_d("Wind speed") +
coord_polar()
Même dominance avec des vents qui viennent un peu plus du nord, mais qui restent quand même orientés nord-est.
Pluie
Vérifions à nouveau la cohérence entre l’humidité et la pluie:
infoclimat_t %>%
filter(Pluie > 0) %>%
ggplot(aes(x = Humidité)) +
stat_ecdf()
Les deux-tiers des journées pluvieuse présentent une humidité maximale supérieure à 90%.
Cohérence avec meteociel
Comparons que les données accessible sur meteociel.fr sont les mêmes.
meteociel %>%
left_join(infoclimat_t, by = c("Heurelocale" = "Heure")) %>%
mutate(Pression = Pression.x - Pression.y,
Température = Température.x - Température.y,
Humidité = Humidité.x - Humidité.y,
Pluie = Pluie - `Précip./h`,
Vent = Vent - vent_heure) %>%
select(Pression, Température, Humidité, Pluie, Vent) %>%
summary()## Pression Température Humidité Pluie Vent
## Min. :0 Min. :0 Min. :0.00000 Min. :0 Min. :-2.0000
## 1st Qu.:0 1st Qu.:0 1st Qu.:0.00000 1st Qu.:0 1st Qu.:-1.0000
## Median :0 Median :0 Median :0.00000 Median :0 Median : 0.0000
## Mean :0 Mean :0 Mean :0.02273 Mean :0 Mean :-0.4214
## 3rd Qu.:0 3rd Qu.:0 3rd Qu.:0.00000 3rd Qu.:0 3rd Qu.: 0.0000
## Max. :0 Max. :0 Max. :1.00000 Max. :0 Max. : 2.0000
## NA's :2 NA's :2 NA's :2 NA's :7 NA's :2meteociel %>%
left_join(infoclimat_t, by = c("Heurelocale" = "Heure")) %>%
filter(Vent != 0) %>%
ggplot(aes(x = direction.x, y = direction.y)) +
geom_jitter(size = 1)## Warning: Removed 2 rows containing missing values (geom_point).
C’est correct, à quelques pourcents près.
Comparaison et cohérence entre les 2 stations météos
Comparons la proximité entre les données des deux stations météorologiques de Besançon.
meteo <- infoclimat %>%
full_join(infoclimat_t, by = "Heure")
meteo %>%
mutate(Pression = Pression.x - Pression.y,
Température = Température.x - Température.y,
Humidité = Humidité.x - Humidité.y,
Pluie = ifelse(pluie_heure == 0 & Pluie == 0, NA, pluie_heure - Pluie),
Direction = ifelse(vent_heure.x == 0 & vent_heure.y == 0, NA,
pmin((direction.x - direction.y + 360) %% 360,
(direction.y - direction.x + 360) %% 360)),
Vent = ifelse(vent_heure.x == 0 & vent_heure.y == 0, NA,
vent_heure.x - vent_heure.y)) %>%
select(Pression, Température, Humidité, Pluie, Direction, Vent) %>%
summary()## Pression Température Humidité Pluie
## Min. :-1.700 Min. :-3.800 Min. :-17.000 Min. :-6.400
## 1st Qu.:-0.875 1st Qu.:-0.500 1st Qu.: -3.000 1st Qu.:-0.350
## Median :-0.400 Median : 0.100 Median : 2.000 Median :-0.200
## Mean :-0.368 Mean : 0.335 Mean : 1.529 Mean :-0.306
## 3rd Qu.: 0.100 3rd Qu.: 1.000 3rd Qu.: 6.000 3rd Qu.: 0.000
## Max. : 1.400 Max. : 4.100 Max. : 17.000 Max. : 3.200
## NA's :11023 NA's :11023 NA's :11023 NA's :13038
## Direction Vent
## Min. : 0.00 Min. :-9.000
## 1st Qu.: 15.00 1st Qu.:-3.000
## Median : 30.00 Median :-1.000
## Mean : 38.31 Mean : 0.203
## 3rd Qu.: 52.00 3rd Qu.: 3.000
## Max. :178.00 Max. :14.000
## NA's :11353 NA's :11353Les données les plus continues et qui sont susceptibles de peu dépendre de la localité sont la pression, la température et l’humidité. On observe effectivement une grande cohérence entre les pressions et des valeurs assez proches pour la température et l’humidité bien qu’il y ait des différences de quelques degrés pour la température et jusqu’à plus de 10 % de différence pour l’humidité.
meteo %>%
drop_na(Pression.x, Pression.y, Température.x, Température.y, Humidité.x,
Humidité.y) %>%
summarise(corr_pression = cor(Pression.x, Pression.y),
corr_temp = cor(Température.x, Température.y),
corr_humi = cor(Humidité.x, Humidité.y))## # A tibble: 1 x 3
## corr_pression corr_temp corr_humi
## <dbl> <dbl> <dbl>
## 1 0.995 0.988 0.970Les corrélations confirment ces similarités.
Étant un phénomène ponctuel, la pluie demandent à être analysée plus finement. Par ailleurs, même si pour la moitié des relevés, les forces du vent sont à 3 km/h l’un de l’autre et les directions à moins de 30 °, les valeurs restantes peuvent beaucoup différer.
Pluie
Trouvons une période ou les données sur la pluie diffèrent le plus.
(pluie_coher <- meteo %>%
mutate(Heure_heure = ceiling_date(Heure, "1 hour")) %>%
group_by(Heure_heure) %>%
summarise(pluie_heure = max(0, pluie_heure, na.rm = TRUE),
pluie_inst_agg = sum(pluie_inst, na.rm = TRUE) / 6,
Pluie = max(0, Pluie, na.rm = TRUE)) %>%
mutate(day = date(Heure_heure)) %>%
group_by(day) %>%
summarise(pluie = sum(pluie_heure, na.rm = TRUE),
pluie_inst = sum(pluie_inst_agg, na.rm = TRUE),
pluie_t = sum(Pluie, na.rm = TRUE)) %>%
arrange(desc(abs(pluie - pluie_t))))## # A tibble: 93 x 4
## day pluie pluie_inst pluie_t
## <date> <dbl> <dbl> <dbl>
## 1 2018-08-09 0 2.67 21.9
## 2 2018-09-06 10.6 13.6 23.6
## 3 2018-08-13 6 4.67 12.8
## 4 2018-08-23 4 3 9
## 5 2018-08-25 5.6 8.47 9.9
## 6 2018-09-21 5.6 6.53 8.9
## 7 2018-08-14 1.8 2.17 3.9
## 8 2018-09-23 7.6 17.9 9.6
## 9 2018-09-07 0.6 0 2.5
## 10 2018-10-30 7.2 2.6 5.8
## # … with 83 more rowspluie_coher %>%
summarise(corr1 = cor(pluie, pluie_t), corr2 = cor(pluie_inst, pluie_t))## # A tibble: 1 x 2
## corr1 corr2
## <dbl> <dbl>
## 1 0.777 0.761pluie_coher %>%
gather(type, valeur, pluie, pluie_inst) %>%
ggplot(aes(x = pluie_t, y = valeur)) +
geom_point() +
facet_wrap(~ type) +
stat_smooth(method = "lm")
Le 9 août est le jour où les différences sont le plus marquées. Observons les données pluviométriques sur cette journée.
meteo %>%
gather(mesure, valeur, pluie_heure, pluie_inst, Pluie) %>%
filter(Heure > make_date(2018, 8, 9)) %>%
filter(Heure < make_date(2018, 8, 10)) %>%
ggplot(aes(x = Heure, y = valeur)) +
facet_wrap(~ mesure, scales = "free_y", nrow = 3) +
geom_point()## Warning: Removed 358 rows containing missing values (geom_point).
La différence est essentiellement dû à des données manquantes. Cependant, les précipitations peuvent sensiblement varier en fonction de l’endroit, ce qui confirme l’utilité surtout qualitative de cette donnée (journée pluvieuse ou non) avec éventuellement une indication générale de la quantité des précipitations (beaucoup, moyenne, peu).
Vent
Pour le vent, la présence de deux variables complique l’analyse. Commençons par la force.
meteo %>%
filter(!is.na(vent_heure.x) & !is.na(vent_heure.y)) %>%
filter(vent_heure.x != 0 | vent_heure.y != 0) %>%
summarise(corr = cor(vent_heure.x, vent_heure.y))## # A tibble: 1 x 1
## corr
## <dbl>
## 1 0.593meteo %>%
filter(vent_heure.x != 0 | vent_heure.y != 0) %>%
ggplot(aes(x = vent_heure.x, y = vent_heure.y)) +
geom_jitter(size = 1) +
geom_abline(slope = 1)## Warning: Removed 7955 rows containing missing values (geom_point).
La force du vent présente une forte variabilité et une faible corrélation entre les stations.
Étudions la direction observée.
meteo %>%
filter(vent_heure.x != 0 & vent_heure.y != 0) %>%
ggplot(aes(x = direction.x, y = direction.y, color = vent_heure.x)) +
geom_jitter(size = 1) +
geom_hline(yintercept = (0:8) * 360/8 + 360/16) +
geom_vline(xintercept = (0:8) * 360/8 + 360/16) +
geom_abline(slope = 1) +
scale_color_viridis_c()
On remarque quelques données proches de 360 ° assez alignées.
meteo %>%
filter(direction.x > 350) %>%
filter(200 < direction.y & direction.y < 300) %>%
select(direction.x, vent_heure.x, vent_heure.y) %>%
summary()## direction.x vent_heure.x vent_heure.y
## Min. :360 Min. : 0.000 Min. : 4.000
## 1st Qu.:360 1st Qu.: 2.000 1st Qu.: 4.000
## Median :360 Median : 5.000 Median : 7.000
## Mean :360 Mean : 6.476 Mean : 6.524
## 3rd Qu.:360 3rd Qu.: 8.000 3rd Qu.: 7.000
## Max. :360 Max. :18.000 Max. :11.000Cela ne ressemble pas à une erreur liée à des vents de faible amplitude. C’est peut-être des erreurs de mesures. Peu de données semblent impactées.
Conclusion
La pression semble extrêmement stable, suivi de la température et de l’humidité. Pour les précipitations, il faut sans doute agréger les données en quelques catégories sur une journée entière. Pour le vent, à la fois la direction et la force sont assez stables la moitié du temps, mais les différences peuvent cependant être importantes.
## R version 3.6.1 (2019-07-05)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.3 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/libopenblasp-r0.2.20.so
##
## locale:
## [1] LC_CTYPE=fr_FR.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=fr_FR.UTF-8 LC_COLLATE=fr_FR.UTF-8
## [5] LC_MONETARY=fr_FR.UTF-8 LC_MESSAGES=fr_FR.UTF-8
## [7] LC_PAPER=fr_FR.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=fr_FR.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] tidyr_0.8.3 readr_1.3.1 lubridate_1.7.4 ggplot2_3.2.0
## [5] dplyr_0.8.3
##
## loaded via a namespace (and not attached):
## [1] Rcpp_1.0.1 pillar_1.4.2 compiler_3.6.1
## [4] tools_3.6.1 zeallot_0.1.0 digest_0.6.20
## [7] evaluate_0.14 tibble_2.1.3 gtable_0.3.0
## [10] viridisLite_0.3.0 pkgconfig_2.0.2 rlang_0.4.0
## [13] cli_1.1.0 yaml_2.2.0 xfun_0.8
## [16] withr_2.1.2 stringr_1.4.0 knitr_1.23
## [19] vctrs_0.2.0 hms_0.5.0 grid_3.6.1
## [22] tidyselect_0.2.5 glue_1.3.1 R6_2.4.0
## [25] fansi_0.4.0 rmarkdown_1.14 purrr_0.3.2
## [28] magrittr_1.5 scales_1.0.0 backports_1.1.4
## [31] htmltools_0.3.6 assertthat_0.2.1 colorspace_1.4-1
## [34] labeling_0.3 utf8_1.1.4 stringi_1.4.3
## [37] lazyeval_0.2.2 munsell_0.5.0 crayon_1.3.4