Tag Archives: scrap

October, grant proposal season

09/10/2018 Arthur Charpentier Leave a comment

In 2012, Danielle Herbert, Adrian Barnett, Philip Clarke and Nicholas Graves published an article entitled “on the time spent preparing grant proposals: an observational study of Australian researchers“, whose conclusions had been included in Nature under a more explicit title, “Australia’s grant system wastes time” ! In this study, they included 3700 grant applications sent to the National Health and Medical Research Council, and showed that each application represented 37 working days: “Extrapolating this to all 3,727 submitted proposals gives an estimated 550 working years of researchers’ time (95% confidence interval, 513-589)“. But in these times when I have to write my funding application, I find that losing 37 days of work is huge. Because it’s become the norm! And somehow, it’s sad.

Forget about the crazy idea that I would rather, in fact, spend more time doing my research. In fact, the thought I had this morning was that it is rather sad that in the Faculty of Science, mathematicians are asked to spend a considerable amount of time, comparable to that required of physicists or chemists, for often smaller amounts of funding… And I thought it could be easily verified. We start by retrieving the discipline codes

url="http://www.nserc-crsng.gc.ca/NSERC-CRSNG/FundingDecisions-DecisionsFinancement/ResearchGrants-SubventionsDeRecherche/ResultsGSC-ResultatsCSS_eng.asp"
download.file(url,destfile = "GSC.html")
library(XML)
tables=readHTMLTable("GSC.html")
GSC=tables[[1]]$V1
GSC=as.character(GSC[-(1:2)])
namesGSC=tables[[1]]$V2
namesGSC=as.character(namesGSC[-(1:2)])

We’re going to need a small function, to remove the $ and other symbols that pollute the data (and prevent them from being treated as numbers)

library(stringr)
Correction = function(x) as.numeric(gsub('[$,]', '', x))

We will now read the 12 pages, and harvest (we will just take the 2017 data, but we could go back a few years before)

grants= function(gsc){
     url=paste("http://www.nserc-crsng.gc.ca/NSERC-CRSNG/FundingDecisions-DecisionsFinancement/ResearchGrants-SubventionsDeRecherche/ResultsGSCDetail-ResultatsCSSDetails_eng.asp?Year=2017&amp;GSC=",gsc,sep="")
    download.file(url,destfile = "GSC.html")
    library(XML)
    tables=readHTMLTable("GSC.html")
    X=as.character(tables[[1]]$"Awarded Amount")
    A=as.numeric(Vectorize(Correction)(X))
return(c(median(A),mean(A),as.numeric(quantile(A,(1:99)/100))))
}
M=Vectorize(grants)(GSC[1:12])

The average amounts of individual grants can be compared,

barplot(M[2,])

In mathematics, the average grant amount is $24400. If we normalize by this quantity, we obtain

barplot(M[2,]/M[2,8])

In other words, the average amount of a (individual) grant in chemistry (to pay for students, conferences, etc.) is twice that in mathematics, 60% higher in physics than in maths…

We can also look at the median values (rather than the averages)

barplot(M[1,])

Here again, it is in mathematics that it is the weakest….

barplot(M[1,]/M[1,8])

in comparable proportions. If we think that the time spent writing should be proportional to the amount allocated, we should spend half as much time in math as in chemistry.

Cumulative functions can also be ploted,

plot(M[3:101,8],(1:99)/100,type="s",xlim=range(M))
lines(M[3:101,5],(1:99)/100,type="s",col="red")
lines(M[3:101,4],(1:99)/100,type="s",col="blue")

with math in black, physics in red, and chemistry in blue. What is surprising is the bottom part: a “bad” researcher in chemistry or physics will earn more than the median researcher in mathematics…

Now that my intuition is confirmed, I have to go back, writing my proposal… and explain to my coauthors that I have to postpone some research projects because, well, you know…

Computer

Scraper, ou pas ?

19/09/2018 Arthur Charpentier Leave a comment

Ce matin, je mettais en ligne un billet scraper la base d’incendies de forêts expliquant comment scraper la base Promethee, en remplaçant automatique une formulaire. Dans la soirée, Alek Vladimir Nevski sur Twitter me faisait remarquer que c’était intéressant, mais peut être inutile (sur cet exemple particulier en tous cas).

Sur Firefox, il existe un “Moniteur Réseau“, qui s’ouvre en cliquant sur le clavier Ctrl + Maj + E ( Command + Option + E quand on est est sur un Mac). Ça ouvre un espace en bas de la fenêtre

Si maintenant je lance a la main ma requête, on voit ce qui est fait

Une methode GET apparait. En cliquant dessus, on a toute l’information nécessaire a droite

Dans l’onglet ‘Reponse’, on voit l’information sur le fichier json créé,

on a plus d’information bien entendu, si on agrandit

et surtout si on va dans l’onglet ‘En-tetes’

on a l’URL complet de la requête ! On peut ensuite copier le lien, et l’ouvrir

On a ainsi a accès a presque toute la table. Quand on regarde l’URL, seuls 20 incendies sont renvoyés, a cause du &nbrLigne=20& mais on peut tenter de le changer. Alek Vladimir Nevski propose de mettre &nbrLigne=10000& et il affirmait que ça marchait. J’avoue avoir essayé, sans succès.

Mais je retiens le tuyau en tous cas. Parfois, on peut faire simple !

Computer

Scraper la base d’incendies de forêts

19/09/2018 Arthur Charpentier 2 Comments

Aller, deux billets pour parler de scraping, ou harvesting serait peut être plus juste ? Disons qu’on va faire un code automatique pour récupérer des informations en ligne. Plus particulièrement, on va s’intéresser au cas ou les données sont obtenues par une requête manuelle, et le site n’a pas d’API pour générer un joli fichier csv. Pierre-Marie, de la formation Data Science pour l’Actuariat a propose un exemple dont je vais m’inspirer (avec une autre methode pour scraper car j’ai eu des soucis pour installer le logiciel).

En 1973 a été lancée en France Prométhée, une base de données sur les incendies de forêts de la région méditerranéenne (couvrant 15 départements du Sud-Est). Oui, Προμηθεύς est ce titan qui a voulu jouer avec le feu… On peut consulter la base sur http://promethee.com/default/incendies,

En remplissant les champs, on peut faire une extraction suivant tel ou tel critère. Par exemple, en 2016, en Corse, il y a eu 6 incendies de nature “accidentelle”,

On le voit, le lien de la page ne change pas… il faut donc trouver une methode pour remplir manuellement les champs, et récupérer les informations. En particulier, on voit qu’on peut avoir accès aux informations relatives aux différents sinistres en cliquant sur la loupe

Idéalement, on voudrait un code R, automatique, pour faire ça ! Une première solutions, proposée par Pierre-Marie, est de passer par RSelenium. Le package R (non maintenu sur le CRAN depuis quelques semaines) est base sur Selenium, “Selenium automates browsers” comme dit la publicité. Le code pour scraper une page comme par

require("devtools")
devtools::install_github("ropensci/RSelenium")
require(RSelenium)
remDr = remoteDriver(remoteServerAddr = "localhost" 
                      , port = 4445L
                      , browserName = "firefox")

Le soucis est que ce code ne tourne pas sur mon mac. Il y a de nombreuses pages qui parlent des spécificités de mac (en particulier le soucis du port par défaut 4445, qu’il faut changer). J’ai tout essayé, sans succès.

Une autre piste est de passer par un autre package

require("devtools")
devtools::install_github("ropensci/wdman")

wdman, pour “Webdriver Manager“, propose aussi de le faire, en s’appuyant sur diverses technologies, selenium, chromedriver, phantomJS binary ou encore geckodriver. Moyennant l’installation de chrome, la seconde option est intéressante

library(wdman)
cDrv = chrome()
eCaps = list(chromeOptions = list(
  args = c('--headless', '--disable-gpu', '--window-size=1280,800')
))
remDr = remoteDriver(browserName = "chrome", port = 4567, 
                     extraCapabilities = eCaps)
remDr$open()

qui marche ! On peut lui demander d’aller sur la base promethee, et de faire une copie d’écran

remDr$navigate("http://www.promethee.com/default/incendies")
remDr$screenshot(display = TRUE)

On va ensuite remplir les champs, pour notre requête, sur la page. Par exemple, pour la date de début, on peut remonter au début des années 80

webElem = remDr$findElement(using = 'name', 'dtAlerteDeb')
webElem$clearElement()
webElem$sendKeysToElement(list('01/01/1981'))

et terminer fin 2017

webElem = remDr$findElement(using = 'name', 'dtAlerteFin')
webElem$clearElement()
webElem$sendKeysToElement(list('31/12/2017'))

On se focalise seulement sur la Corse

webElem = remDr$findElement(using = 'name', 'codeg')
webElem$sendKeysToElement(list('CORSE'))

et on veut 50 observations par page

webElem = remDr$findElement(using = 'name', 'nbrLigne')
webElem$sendKeysToElement(list('50'))

Ah oui, et on va temporiser entre chaque requête

webElem = remDr$findElement(using = 'name', 'btnSubmit')
webElem$submitElement()
Sys.sleep(15)
remDr$screenshot(TRUE)

(c’est une règle assez classique quand on scrap: on évite d’aller trop vite, pour pas se faire repérer…). On va ensuite commencer a regarder ce que l’on demande

text_base = read_html(remDr$getPageSource()[[1]])
list_tr = html_nodes(text_base, "tr")
a = html_nodes(list_tr[length(list_tr)], "a")
page_max = as.numeric(html_text(a[10]))

On a environ 620 pages ! avec au moins 10 secondes pour changer de pages, et en rajoutant 15 secondes (on en parlait avant), ca fait quelques heures pour scraper.

Que veut-on récupérer en grattant la page ? On peut récupérer les coordonnées de l’incendie (latitude/longitude), les caractéristiques de l’incendie (dates, lieu, surface), et retourner le tout dans une table

scrapping = function(feux, page) {
  list_a = html_nodes(read_html(feux), "a")
  if (length(list_a) &gt;= 2) {
    att = html_attrs(list_a[2])
    longitude = as.numeric(str_extract(att, "[0-9]\\.([0-9])*"))
    latitude = as.numeric(str_extract(att, "([0-9]){2}\\..([0-9])*"))
  } else {
    longitude = 0
    latitude = 0
  }
  list_span = html_nodes(read_html(feux), "span")
  date = html_text(list_span[2])
  commune = html_text(list_span[4])
  surface = as.numeric(html_text(list_span[5]))  
  list_div = html_nodes(read_html(feux), "div")
  if (length(list_a) &gt;= 2) {
    infos1 = html_text(list_div[3])
    infos2 = html_text(list_div[5])
  } else {
    infos1 = html_text(list_div[2])
    infos2 = html_text(list_div[4])
  }
  alerte = gsub("First alert : |First intervention", "", str_extract(infos1, "First alert : (.)*First intervention"))
  deb_inter = gsub("First intervention : |End of intervention", "", str_extract(infos1, "First intervention : (.)*End of intervention"))
  fin_inter = gsub("End of intervention : |Origins of the alert", "", str_extract(infos1, "End of intervention : (.)*Origins of the alert"))
  code_INSEE = gsub("INSEE code : |Locality", "", str_extract(infos2, "INSEE code : (.)*Locality"))
  lieudit = str_trim(gsub("Locality : ", "", str_extract(infos2, "Locality : (.)*")))
  liste_feux = data.frame(commune, date, code_INSEE, longitude, latitude, lieudit, surface, alerte, deb_inter, fin_inter, page)
  return(liste_feux)
}

Reste a faire tourner sur les 620 pages. Pour chaque ligne, on récupère les informations !

df = data.frame()
for (page in 1:page_max) {
  text_base = read_html(remDr$getPageSource()[[1]])
  list_tr = html_nodes(text_base, "tr")
  if (length(html_nodes(read_html(as.character(list_tr[14])), "a")) == 8) {
    deb = 15
  } else {
    deb = 14
  }
  for (i in deb:(deb+49)) {
    feux = as.character(list_tr[i])
    df = rbind(df, scrapping(feux, page))
  }
  if (page == 1) {
    webElem = remDr$findElement(using = 'xpath', '//*[contains(concat( " ", @class, " " ), concat( " ", "page_img", " " ))]')
  } else {
    webElem = remDr$findElement(using = 'xpath', '//*[contains(concat( " ", @class, " " ), concat( " ", "page_img", " " ))][3]')
  }
  webElem$clickElement()
  heure = (floor((page_max-page)*15/60) - floor((page_max-page)*15/60) %% 60) / 60
  minute = floor((page_max-page)*15/60) %% 60
  pb = progress_bar$new(total = 100)
  for (i in 1:100) {
    if (i==1) {cat("page", page, "sur", page_max," -  (reste", heure, "h", minute, "min)",  "\n")}
    pb$tick()
    Sys.sleep(15 / 100)
  }
  remDr$screenshot(TRUE)
}
feux = df

Ça y est, on a une base, avec presque 31 000 lignes ! on va ensuite ranger un peu, et clarifier

colnames(feux)[colnames(feux)=="longitude"] = "x"
colnames(feux)[colnames(feux)=="latitude"] = "y"
feux[, date := as.Date(date, "%d/%m/%Y")]
feux[, date_alerte := as.Date(substr(alerte, 1, 10), "%d/%m/%Y")]
feux$heure_alerte[nchar(feux$alerte) &gt; 10] = substr(feux$alerte[nchar(feux$alerte) &gt; 10], 12, 16)
feux$heure_alerte[!(nchar(feux$alerte) &gt; 10)] = ""
feux[, date_deb_inter := substr(deb_inter, 1, 10)]
feux[, date_deb_inter := as.Date(date_deb_inter, "%d/%m/%Y")]
feux$heure_deb_inter[nchar(feux$deb_inter) &gt; 10] = substr(feux$deb_inter[nchar(feux$deb_inter) &gt; 10], 12, 16)
feux$heure_deb_inter[!(nchar(feux$deb_inter) &gt; 10)] = ""
feux$heure_deb_inter[feux$heure_deb_inter == "00:00"] = ""
feux[, date_fin_inter := substr(fin_inter, 1, 10)]
feux[, date_fin_inter := as.Date(date_fin_inter, "%d/%m/%Y")]
feux$heure_fin_inter[nchar(feux$fin_inter) &gt; 10] = substr(feux$fin_inter[nchar(feux$fin_inter) &gt; 10], 12, 16)
feux$heure_fin_inter[!(nchar(feux$fin_inter) &gt; 10)] = ""
feux$heure_fin_inter[feux$heure_fin_inter == "00:00"] = ""
feux$annee = year(feux$date_alerte)
feux$month = month(feux$date_alerte)

Pour voir un peu nos données, rien de tel que des statistiques descriptives. On peut commencer avec du barplot du nombre de feux de forêts par année

agg_annee = aggregate(alerte~annee, FUN=length, data=feux)
colnames(agg_annee)[colnames(agg_annee)=="alerte"] = "Nombre"
nb = ggplot(agg_annee) +
  geom_bar(aes(x=annee, y=Nombre, fill="Feux de forêts"), stat="identity", width = 0.8) +
  scale_colour_manual(name=NULL, values=c("orange")) +
  ggtitle("Nombre de feux par an") +
  theme(plot.title = element_text(hjust = 0.5, size=14),
        legend.position = "none",
        axis.text.x = element_text(angle=50, vjust=0.5, size=10),
        axis.text.y = element_text(size=12),
        axis.title.x = element_blank(),
        axis.title.y = element_blank()) +
  scale_fill_discrete("") +
  scale_x_discrete(limit=seq(1981, 2017, by = 2))

ou un barplot de la surface de feux de forêts par anee

agg_surface = aggregate(surface~annee, FUN=sum, data=feux)
surf = ggplot(agg_surface) +
  geom_bar(aes(x=annee, y=surface, fill="Feux de forêts"), stat="identity", width = 0.8) +
  scale_colour_manual(name=NULL, values=c("orange")) +
  ggtitle("Surface totale incendiée par an (en ha)") +
  theme(plot.title = element_text(hjust = 0.5, size=14),
        legend.position = "none",
        axis.text.x = element_text(angle=50, vjust=0.5, size=10),
        axis.text.y = element_text(size=12),
        axis.title.x = element_blank(),
        axis.title.y = element_blank()) +
  scale_fill_discrete("") +
  scale_x_discrete(limit=seq(1981, 2017, by = 2)) +
  scale_y_discrete(limit=seq(0, 25000, by = 5000))

Et on fait le graphique

require(grid)
grid.newpage()
pushViewport(viewport(layout = grid.layout(1, 2)))
define_region = function(row, col){
  viewport(layout.pos.row = row, layout.pos.col = col)
} 
print(nb, vp = define_region(1, 1))
print(surf, vp = define_region(1, 2))

Passons aux choses serieuses, avec des cartes (la partie la plus serieuse, on en conviendra, c’etait la phase initiale de scraping). On commence par des informations relatives a la Corse

url="https://raw.githubusercontent.com/data-PM/Promethee/master/communes_corse-shp.zip"
temp="corse.zip"
download.file(url, "corse.zip")
unzip("corse.zip")
communes_corse = sp::spTransform(shapefile("/Users/arthur/communes_corse-shp/communes_corse.shp"), CRS("+proj=longlat +datum=WGS84"))

On se limite aux feux de plus de 1000m2

feux_min1000m2 = feux[feux$surface &gt; 0.1, ]

Maintenant, on peut faire la carte

feux_INSEE = aggregate(x = feux_min1000m2$code_INSEE, by = list(unique.values = feux_min1000m2$code_INSEE), FUN = length)
colnames(feux_INSEE)[colnames(feux_INSEE)=="unique.values"] = "insee"
communes_corse = merge(x = communes_corse, y = feux_INSEE, by = "insee", all.x = TRUE)
communes_corse$x[is.na(communes_corse$x)] = 0
max = max(communes_corse$x)
communes_corse$x2 = communes_corse$x/max * 100
corse2 = gSimplify(communes_corse, tol = 0.001, topologyPreserve = TRUE)
corse2 = SpatialPolygonsDataFrame(corse2, communes_corse@data, match.ID = FALSE)
opar = par(mar = c(0,0,1.2,0))
plot(corse2, col = "#DAE3E6", border = "#8A0641", lwd = 0.7, bg = "#B5D0D0")

Le fond de carte (avec les communes) ressemble a ca

On peut maintenant passer a la carte choroplethe

choroLayer(spdf = corse2, df = corse2@data, var = "x", border = NA, 
           col = carto.pal("wine.pal", 6), legend.pos = "topleft", 
           add = TRUE, method = "fisher-jenks", nclass = 6, legend.title.txt = NULL)
top_villes = corse2[corse2$nom %in% c("Ajaccio", "Bastia", "Corte", "Porto-Vecchio", "Propriano", "Bonifacio", "Calvi"), ]
graphics::points(coordinates(top_villes), pch = 20, cex = 0.5)
labelLayer(top_villes, top_villes@data, txt = "nom", cex = 1, pos = 2, font = 4, offset = 0.5)
layoutLayer(title = "Nombre de feux de forêts depuis 1981 (minimum 1000 m2)",
            source = "Source : 1973-2018 Prométhée ©",
            author = NULL, scale = 10, frame = FALSE, col = "#cdd2d4", coltitle = "#8A5543")

Comme souvent avec R, plusieurs formats sont envisageables. On peut passer a du ggplot2, mais on commence par aller chercher un shapefile global de la Corse

corse = shp_zip("https://raw.githubusercontent.com/data-PM/Promethee/master/corse-shp.zip", "corse")
corse = spTransform(corse, CRS("+proj=longlat +datum=WGS84"))
corse = fortify(corse)
p=data %&gt;%
  ggplot() +
  geom_polygon(data=corse, aes(x=long, y = lat, group = group), fill="grey", alpha=0.3) +
  geom_point(data=feux_min100ha, aes(x=x, y=y, size=surface, color=surface, text=mytext, alpha=surface) ) +
  labs(title = "\nPrincipaux feux de forêts depuis 1981 \n (plus de 100 Ha)") +
  xlim(8.5, 9.6) + ylim(41.3,43.1) +
  scale_size_continuous(range=c(1,5), ) +
  scale_color_viridis(option="inferno", trans="log", name="Surface", limits=c(99, 5700), breaks=c(100, 200, 500, 1000, 2500, 5000)) +
  scale_alpha_continuous(trans="log") +
  theme_void() +
  theme(legend.title = element_text(size=10, face="bold"))
t=800
div(ggplotly(tooltip="text", width=t/1.7, height=t), align = "center")

Joli, non ?

Computer

Scraper pour avoir des infos sur les médecins sur Paris

08/10/2017 Arthur Charpentier 2 Comments

Allez, quelques codes R pour s’amuser un peu… Je me suis inspiré du projet de Marine pour la formation Data Science pour l’Actuariat. L’idée était de scraper un site qui permet de prendre un rendez-vous chez un médecin, pour voir où se trouvent les médecins (sur Paris) et de voir leurs prix. Histoire de donner un avant goût, voilà les graphiques qu’on cherche à produire

Continue reading Scraper pour avoir des infos sur les médecins sur Paris →

Freakonometrics

Tag Archives: scrap

October, grant proposal season

Scraper, ou pas ?

Scraper la base d’incendies de forêts

Scraper pour avoir des infos sur les médecins sur Paris

"sendo l'intento mio scrivere cosa utile a chi la intende…"