{"word_count": 698, "character_count": 10509}
Allez, un peu de scraping aujourd'hui (promis, il y aura d'autres billets d'ici la fin de la semaine sur le sujet). Je vais m'inspirer des codes de Romain de la formation Data Science pour l’Actuariat. Le but est de récupérer des données de températures, de précipitations, de vitesse de vent, journalières, en France. On commence par charger quelques librairies,library(plyr)library(stringr)library(OpenStreetMap)library(leaflet)library(shiny)library(rsconnect)library(mapview)library(png)library(magick)library(yaml)On va récupérer les deux dernières années, mais comme toujours, le code peut facilement s'adapter.annee = c(2016:2017)mois = c('01','02','03','04','05','06','07','08','09','10','11','12')aaaamm = sort(sub(pattern=' ',replacement = '', x=outer(annee,mois,paste)))myCols = c(NA,NA,"NULL","NULL","NULL",NA,NA,"factor","NULL","NULL","NULL",NA,"NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL",NA,NA,NA,NA,NA,"NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL")for (mois in aaaamm) {  fichier = paste0("synop_",mois,".csv.gz")  if (!file.exists(fichier)) {    page = paste0('https://donneespubliques.meteofrance.fr/donnees_libres/Txt/Synop/Archive/synop.',mois,'.csv.gz')    download.file(page, fichier, quiet = TRUE, cacheOK = FALSE)  }}  for (mois in aaaamm) {  fichier = paste0("synop_",mois,".csv.gz")  if (mois=='201601') {    data = read.table(gzfile(fichier),header=TRUE, sep=";",quote='"',fileEncoding="UTF-8",colClasses = myCols)  } else {    data = rbind(data,read.table(gzfile(fichier),header=TRUE, sep=";",quote='"',fileEncoding="UTF-8",colClasses = myCols))  }}On supprime une des colonnes qui ne sert a riendata[,"X"] = NULLet on bricole un peu (sinon on a des soucis avec minuit, qui est décalé de 24 heures)data$temp = as.numeric(substr(data$date,1,14))decalage_minuit = function(dateheure){  r = dateheure  if (dateheure-trunc(dateheure/1000000)*1000000==0){    r = dateheure + 10000  }  return(r)}data$temp = sapply(data$temp,decalage_minuit)data$jourheure_obs = strptime(data$temp, format = '%Y%m%d%H%M%S', 'UTC')data$date = as.POSIXct((strptime(data$jourheure_obs,format = "%Y-%m-%d")))data$heure = as.numeric(substr(data$jourheure_obs, start=12, stop=13))data[,"temp"] = NULLdata$t = as.numeric(as.character(data$t))-273.15oui, on doit convertir en degrés Celsius. On peut aussi travailler sur les vitesses de ventdata$dd = as.numeric(as.character(data$dd))data$ff = as.numeric(as.character(data$ff))/1000*3600data$ff = round(data$ff/5)*5que l'on va convertir en km/h, et enfin les précipitationsdata$rr1 = pmax(0,as.numeric(as.character(data$rr1)))data$rr3 = pmax(0,as.numeric(as.character(data$rr3)))data$rr6 = pmax(0,as.numeric(as.character(data$rr6)))data$rr12 = pmax(0,as.numeric(as.character(data$rr12)))data$rr24 = pmax(0,as.numeric(as.character(data$rr24)))On peut d'ailleurs commencer par une simple visualisation de séries de températuresdata_temperature = subset(data, !is.na(data$t))temp_min = aggregate(x=data_temperature$t,by=list(date=data_temperature$date),FUN=min)names(temp_min)[2] = "temp_min"temp_max = aggregate(x=data_temperature$t,by=list(date=data_temperature$date),FUN=max)names(temp_max)[2] = "temp_max"temp_moy = aggregate(x=data_temperature$t,by=list(date=data_temperature$date),FUN=mean)names(temp_moy)[2] = "temp_moy"temp_stat = cbind(temp_min,temp_max[],temp_moy)rm(temp_min,temp_max,temp_moy)temp_stat[,3] = NULLtemp_stat[,4] = NULLavec le maximum observé, par jour, sur toutes les stations, le minimum, et la moyenne des stationsplot(x=temp_stat$date,y=temp_stat$temp_moy,'l',ylim=c(-15,40),xlab='Date',ylab='Température (°C)')title('Températures françaises métropolitaines de 01/2016 à 12/2017')points(x=temp_stat$date,y=temp_stat$temp_min,col="blue",'l')points(x=temp_stat$date,y=temp_stat$temp_max,col="red",'l')Pour faire des jolies cartes, il nous faut des informations sur les stations météo.stations = read.table('https://donneespubliques.meteofrance.fr/donnees_libres/Txt/Synop/postesSynop.csv',header=TRUE, sep=";",quote='"',fileEncoding="UTF-8")stations = rename(stations, replace=c("ID"="numer_sta"))stations = rename(stations, replace=c("Nom"="nom"))stations = rename(stations, replace=c("Latitude"="latitude"))stations = rename(stations, replace=c("Longitude"="longitude"))stations = rename(stations, replace=c("Altitude"="altitude"))data = merge(x=data,y=stations)Pour la carte, on va se restreindre aux données en France métropolitainedata = subset(data, data$latitude &#62;= 39)data = subset(data, data$latitude &#60;= 53)data = subset(data, data$longitude &#62;= -5.3)data = subset(data, data$longitude &#60;= 9.8)coordonnees_stations = unique(data[, c("nom","latitude","longitude")])Commencons par visualiser le maximum observee le 19 juillet 2016date_sel = "2016-07-19"date_sel_OK = paste0(substr(date_sel,9,10),"/",substr(date_sel,6,7),"/",substr(date_sel,1,4))donnees_carte_0 = subset(data, data$date==as.POSIXct(date_sel))temp = aggregate(x=donnees_carte_0$t,by=list(numer_sta=donnees_carte_0$numer_sta),FUN="max")donnees_carte_0 = merge(x=donnees_carte_0,y=temp)donnees_carte_0 = rename(donnees_carte_0, replace=c("x"="TEMP"))donnees_carte_0$TEMP = round(donnees_carte_0$TEMP)donnees_carte = unique(donnees_carte_0[,c('numer_sta','nom','longitude','latitude','TEMP')])get_france = get_map(c(lon=2.25,lat=46), zoom=5, col='bw')colfunc = colorRampPalette(c("darkblue","royalblue","cyan","lightblue","orange","red","darkred"))liste_couleurs = data.frame(TEMP=as.numeric(seq(-49,50,1)),couleur=as.character(colfunc(100)))donnees_carte = merge(x=donnees_carte,y=liste_couleurs)france_tmax = ggmap(get_france) +   ggtitle(paste0('Températures maximales du ', date_sel_OK)) +  scale_x_continuous(limits = c(-5, 10), expand = c(0, 0)) +  scale_y_continuous(limits = c(41.3, 51.1), expand = c(0, 0))france_tmax = france_tmax + geom_label(data=donnees_carte,aes(x=longitude,y=latitude,label=TEMP),col=as.character(donnees_carte$couleur))france_tmaxOn peut tenter aussi une visualisation des précipitations. On se place le 13 octobre 2016 (oui, il y a eu des inondations dans la region de Montpelier ce jour la)date_sel = "2016-10-13"colfunc = colorRampPalette(c('white','darkblue'))data$precip_en_mm = data$rr3data_precipitation = subset(data, !is.na(data$precip_en_mm))precipitations = aggregate(x=data_precipitation$precip_en_mm,by=list(date=data_precipitation$date,station=data_precipitation$nom),FUN=sum)names(precipitations)[3] = "pp"precipitations$pp = pmax(precipitations$pp,0)donnees_carte = subset(precipitations, precipitations$date==date_sel)qt_prec = quantile(round(donnees_carte$pp),seq(0, 1, 0.05))donnees_carte$Couleur_prec = colfunc(21)[(findInterval(round(donnees_carte$pp), qt_prec, all.inside=TRUE))]donnees_carte = rename(donnees_carte, replace=c("station"="nom"))donnees_carte = merge(x=donnees_carte,y=coordonnees_stations)On peut tenter un leaflet ici (ça sera plus joli)france = leaflet() %&#62;% addTiles() %&#62;% fitBounds(lng1=-4.412, lat1=41.92, lng2=9.485, lat2=50.57)france = france %&#62;% addCircles(lng=donnees_carte$longitude, lat=donnees_carte$latitude, color=donnees_carte$Couleur_prec, opacity = 1, fillColor=donnees_carte$Couleur_prec, fillOpacity = 1, radius=donnees_carte$pp*500)franceComme souvent, j'ai des soucis pour integrer du leaflet en wordpress, alors on se contentera dans le billet d'une copie d'ecran (mais un copier/coller du code suffit).Et pour finir, la meme chose avec les vitesses de vent. On se place le 6 mars 2017 (la encore, la date n'est pas choisie au hasard, c’était la tempete Zeus)date_sel = "2017-03-06"date_sel_OK = paste0(substr(date_sel,9,10),"/",substr(date_sel,6,7),"/",substr(date_sel,1,4))heure_sel = 18donnees_carte_0 = subset(data, data$date==as.POSIXct(date_sel))donnees_carte_0 = subset(donnees_carte_0, donnees_carte_0$heure==heure_sel)donnees_carte = donnees_carte_0Pour les precipitations, des cercles plus ou moins grands suffisaient. La, on va mettre des flecheslogo = image_trim(image_read("https://image.freepik.com/icones-gratuites/fleche-noire-vers-le-haut_318-30934.jpg"))if (file.exists('flecheblanche90.png')==FALSE){  for (angle in seq(10,360, by=10)) {    temp = image_rotate(logo,angle)    image_write(temp,paste0('flecheblanche',angle,'.png'))    temp = image_colorize(image=temp,opacity=40,color='green')    image_write(temp,paste0('flecheverte',angle,'.png'))    temp = image_colorize(image=temp,opacity=40,color='yellow')    image_write(temp,paste0('flechejaune',angle,'.png'))    temp = image_colorize(image=temp,opacity=40,color='orange')    image_write(temp,paste0('flecheorange',angle,'.png'))    temp = image_colorize(image=temp,opacity=40,color='red')    image_write(temp,paste0('flecherouge',angle,'.png'))    temp = image_colorize(image=temp,opacity=40,color='black')    image_write(temp,paste0('flechenoire',angle,'.png'))  }}arrowIcons = icons(  iconUrl = ifelse(donnees_carte$dd == 0,            ifelse(donnees_carte$ff&#60;39,                    'flecheverte360.png',            ifelse(donnees_carte$ff&#60;79,                    'flechejaune360.png',            ifelse(donnees_carte$ff&#60;119,                    'flecheorange360.png',            ifelse(donnees_carte$ff&#60;159,                    'flecherouge360.png',            ifelse(donnees_carte$ff&#62;159,                    'flechenoire360.png',                    'flecheblanche360.png'))))),            ifelse(donnees_carte$ff&#60;39,                   paste0('flecheverte',donnees_carte$dd,'.png'),            ifelse(donnees_carte$ff&#60;79,                    paste0('flechejaune',donnees_carte$dd,'.png'),            ifelse(donnees_carte$ff&#60;119,                    paste0('flecheorange',donnees_carte$dd,'.png'),            ifelse(donnees_carte$ff&#60;159,                    paste0('flecherouge',donnees_carte$dd,'.png'),            ifelse(donnees_carte$ff&#62;159,                    paste0('flechenoire',donnees_carte$dd,'.png'),           paste0('flecheblanche',donnees_carte$dd,'.png'))))))),  iconWidth = 30, iconHeight = 30)Cette fois, on est bon ! On peut faire la cartefrance = leaflet(options = leafletOptions()) %&#62;% addTiles() %&#62;% fitBounds(lng1=-4.412, lat1=41.92, lng2=9.485, lat2=50.57)france = addMarkers(map = france,lng=donnees_carte$longitude,lat=donnees_carte$latitude,icon=arrowIcons,popup=as.character(donnees_carte$ff))france