# Extracting information from a picture, round 2

Yesterday, I published a post on extracting information from a picture, but it did not work as expected. I claimed that it was because of the original graph I had. More precisely, the was based on some weird projection, and I could not reconcile. So I decide to cheat a little bit, by creating my own map, Colors are ugly, I know. But I got them using

```u = seq(0,1,length=30) couleurs = rgb(u,rev(u),0,1)```

The picture is

```url = "https://freakonometrics.hypotheses.org/files/2018/12/chomage3.png" library(pixmap) library(png) IMG = readPNG(url)```

I used those colors because it would make things easy when extracting reds and greens…

```ROUGE=t(IMG[,,1])[x1:x2,] ROUGE=ROUGE[,y2:y1] library(scales) image(x1:x2,y1:y2,ROUGE,col=alpha(colour=rgb(1,0,0,1), alpha = seq(0,1,by=.01))) VERT=t(IMG[,,2])[x1:x2,] VERT=VERT[,y2:y1] image(x1:x2,y1:y2,VERT,col=alpha(colour=rgb(0,1,0,1), alpha = seq(0,1,by=.01)))``` Let us see if the contour of France can be overlaid

```library(maptools) library(PBSmapping) download.file("http://biogeo.ucdavis.edu/data/gadm2.8/rds/FRA_adm0.rds","FRA_adm0.rds") FR=readRDS("FRA_adm0.rds") library(maptools) PP = SpatialPolygons2PolySet(FR) par(mfrow=c(1,1)) PP=PP[(PP\$X&lt;=8.25)&amp;(PP\$Y&gt;=42.2),] u=(x1:x2)-x1 v=(y1:y2)-y1 ax=min(PP\$X) bx=max(PP\$X)-min(PP\$X) ay=min(PP\$Y) by=max(PP\$Y)-min(PP\$Y) PP\$X=(PP\$X-ax)/bx*max(u) PP\$Y=(PP\$Y-ay)/by*max(v) image(u,v,ROUGE,col=alpha(colour=rgb(1,0,0,1), alpha = seq(0,1,by=.01))) points(PP\$X,PP\$Y)```

We have a perfect match, don’t we…? Let us now use a shapefile based on départements,

```download.file("http://biogeo.ucdavis.edu/data/gadm2.8/rds/FRA_adm2.rds","FRA_adm2.rds") FR2=readRDS("FRA_adm2.rds") library(maptools) PP = SpatialPolygons2PolySet(FR2) image(u,v,ROUGE,col=alpha(colour=rgb(1,0,0,1), alpha = seq(0,1,by=.01))) k=35 pX=(PP\$X[PP\$PID==k]-ax)/bx*max(u) pY=(PP\$Y[PP\$PID==k]-ay)/by*max(v) points(pX,pY)nge(pX)```

For instance, the thirty-fifth polygon is the following Let us extract the color inside that polygon

```u=1:nrow(ROUGE) v=1:ncol(ROUGE)```

The code would be

```pX=(PP\$X[PP\$PID==k]-ax)/bx*max(u) pY=(PP\$Y[PP\$PID==k]-ay)/by*max(v) E=expand.grid(u,v) M=matrix(point.in.polygon(E[,1],E[,2],pX,pY)&gt;0,length(u),length(v)) image(u,v,ROUGE*M,col=alpha(colour=rgb(1,0,0,1), alpha = seq(0,1,by=.01))) points(pX,pY)``` Now, for each département, I extract the average value of red, and the average value of green,

```extract_info = function(k){ pX=(PP\$X[PP\$PID==k]-ax)/bx*max(u) pY=(PP\$Y[PP\$PID==k]-ay)/by*max(v) E=expand.grid(u,v) M=matrix(point.in.polygon(E[,1],E[,2],pX,pY)&gt;0,length(u),length(v)) nom=FR2[FR2\$OBJECTID ==k,c("NAME_2","CCA_2")] return(c(as.numeric(nom\$CCA_2),sum(ROUGE[M==1])/sum(M),sum(VERT[M==1])/sum(M))) } donnees = Vectorize(extract_info)(1:95) x2=donnees[1,] y2=donnees[2,]/(donnees[2,]+donnees[3,]) df2=data.frame(dpt=x2,extract=y2) x1=as.numeric(as.character(baseChomage\$no)) y1=baseChomage\$chomagePremierTrimestre2017 df1=data.frame(dpt=x1,obs=y1) df=merge(df1,df2) plot(df\$obs,df\$extract)```

On the graph below, we have the original values on the x-axis (unemployement, in percent) and the “average value of red”.  Note that points are almost perfectly correlated… The accumulation can be explained because on the original map, different values could have the same color So far, I can claim that we’ve been able to extract useful information from the original picture.

Consider the case now that the original map was the following one ```url = "https://freakonometrics.hypotheses.org/files/2018/12/chomage5.png" library(pixmap) library(png) IMG = readPNG(url)```

Here, the colors are obtained from a standard palette,

```library(pals) couleurs = rev(brewer.rdylgn(30))```

Here again, we use our previous code to extract reds and greens And if we use our function

```extract_info = function(k){ pX=(PP\$X[PP\$PID==k]-ax)/bx*max(u) pY=(PP\$Y[PP\$PID==k]-ay)/by*max(v) E=expand.grid(u,v) M=matrix(point.in.polygon(E[,1],E[,2],pX,pY)&gt;0,length(u),length(v)) nom=FR2[FR2\$OBJECTID ==k,c("NAME_2","CCA_2")] return(c(as.numeric(nom\$CCA_2),sum(ROUGE[M==1])/sum(M),sum(VERT[M==1])/sum(M))) } donnees = Vectorize(extract_info)(1:95) x2=donnees[1,] y2=donnees[2,]/(donnees[2,]+donnees[3,]) df2=data.frame(dpt=x2,extract=y2) x1=as.numeric(as.character(baseChomage\$no)) y1=baseChomage\$chomagePremierTrimestre2017 df1=data.frame(dpt=x1,obs=y1) df=merge(df1,df2) plot(df\$obs,df\$extract)```

we obtain the following graph Here again, we have a strong correlation, not to say comonotonic variables (in the sense that ranks are identical). Nice, isn’t it ?

# Some sort of Otto Neurath (isotype picture) map

Yesterday evening, I was walking in Budapest, and I saw some nice map that was some sort of Otto Neurath style. It was hand-made but I thought it should be possible to do it in R, automatically.

A few years ago, Baptiste Coulmont published a nice blog post on the package osmar, that can be used to import OpenStreetMap objects (polygons, lines, etc) in R. We can start from there. More precisely, consider the city of Douai, in France, The code to read information from OpenStreetMap is the following

```library(osmar) src &lt;- osmsource_api() bb &lt;- center_bbox(3.07758808135,50.37404355, 1000, 1000) ua &lt;- get_osm(bb, source = src)```

We can extract a lot of things, like buildings, parks, churches, roads, etc. There are two kinds of objects so we will use two functions

```listek = function(vc,type="polygons"){ nat_ids &lt;- find(ua, way(tags(k %in% vc))) nat_ids &lt;- find_down(ua, way(nat_ids)) nat &lt;- subset(ua, ids = nat_ids) nat_poly &lt;- as_sp(nat, type)}   listev = function(vc,type="polygons"){ nat_ids &lt;- find(ua, way(tags(v %in% vc))) nat_ids &lt;- find_down(ua, way(nat_ids)) nat &lt;- subset(ua, ids = nat_ids) nat_poly &lt;- as_sp(nat, type)}```

For instance to get rivers, use

`W=listek(c("waterway"))`

and to get buildings

`M=listek(c("building"))`

We can also get churches

`C=listev(c("church","chapel"))`

but also train stations, airports, universities, hospitals, etc. It is also possible to get streets, or roads

```H1=listek(c("highway"),"lines") H2=listev(c("residential","pedestrian","secondary","tertiary"),"lines")```

but it will be more difficult to use afterwards, so let’s forget about those.

We can check that we have everything we need

```plot(M) plot(W,add=TRUE,col="blue") plot(P,add=TRUE,col="green") if(!is.null(B)) plot(B,add=TRUE,col="red") if(!is.null(C)) plot(C,add=TRUE,col="purple") if(!is.null(T)) plot(T,add=TRUE,col="red")``` Now, let us consider a rectangular grid. If there is a river in a cell, I want a river. If there is a church, I want a church, etc. Since there will be one (and only one) picture per cell, there will be priorities. But first we have to check intersections with polygons, between our grid, and the OpenStreetMap polygons.

```library(sp) library(raster) library(rgdal) library(rgeos) library(maptools) identification = function(xy,h,PLG){ b=data.frame(x=rep(c(xy-h,xy+h),each=2), y=c(c(xy-h,xy+h,xy+h,xy-h))) pb1=Polygon(b) Pb1 = list(Polygons(list(pb1), ID=1)) SPb1 = SpatialPolygons(Pb1, proj4string = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0")) UC=gUnionCascaded(PLG) return(gIntersection(SPb1,UC)) }```

and then, we identify, as follows

```whichidtf = function(xy,h){ h=.7*h label="EMPTY" if(!is.null(identification(xy,h,M))) label="HOUSE" if(!is.null(identification(xy,h,P))) label="PARK" if(!is.null(identification(xy,h,W))) label="WATER" if(!is.null(identification(xy,h,U))) label="UNIVERSITY" if(!is.null(identification(xy,h,C))) label="CHURCH" return(label) }```

Let is use colored rectangle to make sure it works

```nx=length(vx) vx=as.numeric((vx[2:nx]+vx[1:(nx-1)])/2) ny=length(vy) vy=as.numeric((vy[2:ny]+vy[1:(ny-1)])/2) plot(M,border="white") for(i in 1:(nx-1)){ for(j in 1:(ny-1)){ lb=whichidtf(c(vx[i],vy[j]),h) if(lb=="HOUSE") rect(vx[i]-h,vy[j]-h,vx[i]+h,vy[j]+h,col="grey") if(lb=="PARK") rect(vx[i]-h,vy[j]-h,vx[i]+h,vy[j]+h,col="green") if(lb=="WATER") rect(vx[i]-h,vy[j]-h,vx[i]+h,vy[j]+h,col="blue") if(lb=="CHURCH") rect(vx[i]-h,vy[j]-h,vx[i]+h,vy[j]+h,col="purple") }}``` As a first start, we us agree that it works. To use pics, I did borrow them from https://fontawesome.com/. For instance, we can have a tree

``` library(png) library(grid) download.file("http://freakonometrics.hypotheses.org/files/2018/05/tree.png","tree.png") tree &lt;- readPNG("tree.png")```

Unfortunatly, the color is not good (it is black), but that’s easy to fix using the RGB decomposition of that package

``` rev_tree=tree rev_tree[,,2]=tree[,,4]```

We can do the same for houses, churches and water actually

``` download.file("http://freakonometrics.hypotheses.org/files/2018/05/angle-double-up.png","angle-double-up.png") download.file("http://freakonometrics.hypotheses.org/files/2018/05/home.png","home.png") download.file("http://freakonometrics.hypotheses.org/files/2018/05/church.png","curch.png") water &lt;- readPNG("angle-double-up.png") rev_water=water rev_water[,,3]=water[,,4] home &lt;- readPNG("home.png") rev_home=home rev_home[,,4]=home[,,4]*.5 church &lt;- readPNG("church.png") rev_church=church rev_church[,,1]=church[,,4]*.5 rev_church[,,3]=church[,,4]*.5```

and that’s almost it. We can then add it on the map

``` plot(M,border="white") for(i in 1:(nx-1)){ for(j in 1:(ny-1)){ lb=whichidtf(c(vx[i],vy[j]),h) if(lb=="HOUSE") rasterImage(rev_home,vx[i]-h*.8,vy[j]-h*.8,vx[i]+h*.8,vy[j]+h*.8) if(lb=="PARK") rasterImage(rev_tree,vx[i]-h*.9,vy[j]-h*.8,vx[i]+h*.9,vy[j]+h*.8) if(lb=="WATER") rasterImage(rev_water,vx[i]-h*.8,vy[j]-h*.8,vx[i]+h*.8,vy[j]+h*.8) if(lb=="CHURCH") rasterImage(rev_church,vx[i]-h*.8,vy[j]-h*.8,vx[i]+h*.8,vy[j]+h*.8) }}``` Nice, isn’t it? (as least as a first draft, done during the lunch break of the R conference in Budapest, today).

It is now very easy to read (automatically) some text that can be found in a pdf file. For instance, consider the program of the conference we had yesterday – and today – in Rennes ```> library(pdftools) > scan_pdf <- pdf_text("http://crem.univ-rennes1.fr/Documents/Docs_sem_divers/2017_03_10-11_JJD/JDD_prog.pdf") > cat(scan_pdf) Journées Jeunes Docteurs Programme du jeudi 9 mars 2017 Faculty of Economics - Rennes - Amphi Henri Krier 9h- 9h30 - Accueil 9h30-10h15 :      Présentation du CREM, de la faculté et des activités de recherche liées du ou laboratoire 10h15-10h50 :     Emmanuel LORENZON (Université de Bordeaux, GREThA) Collusion with a rent seeking agency in sponsored search auctions 10h50-11h25 :     Julien BERTHOUMIEU (Université de Bordeaux, GREThA) The Impact of “At-the-Border” and “Behind-the-Border” Policies on Cost-Reducing Research and Development Co-écrit avec Antoine Bouët```

(etc). As you can see, it is working well, even in French, where we have those weird letters (with accents). Here, it is working well because the pdf is vectorized, i.e. it was generated properly, by open office.

But sometimes, we can have only a scanned version of a letter or just a picture with some typed text. I will not mention hand-writing because it is much more complex.

The other day, my friend Fleur did show me a picture, and some very simple lines of code, ```> library('tesseract') > pic1="http://freakonometrics.hypotheses.org/files/2017/03/pic1.png" > text_fr <- ocr(pic1, engine = tesseract("fra")) > cat(text_fr) Près de 14.400 décès```

```Si [épidémie de grippe est un phénomène récurrent. celle de 2016—2017 présente plusieurs spéciﬁcités. outre sa virulence : une survenue plus précoce que d‘habitude. une activité modérée en médecine ambulatoire. mais un impact fort en milieu hospitalier.```

It looks like we’ve be able to extract typed text from a picture ! I want to check. I have to admit, first of all that installation on a linux machine is tricky: one has to install first leptonica, and then follow some guidelines to install tesseract (see also Artem‘s advices). It took me some time, but I’ve been able to install the package.

The first important step, it to train the algorithm with some texts in French (because it is in French in my picture)

```> library('tesseract') > tesseract_download("fra")```

Then, I did try with the picture that Fleur did send me (the picture was inserted in the core of the message) ```> pic2="http://freakonometrics.hypotheses.org/files/2017/03/pic2.png" > text_fr <- ocr(pic2, engine = tesseract("fra")) > cat(text_fr) Près de 14400 decès```

```s. mm….agw«… ………«…m …… a……u…u Dhs—ur; ;pmum…. ;: u…… … W»: »… w…q…na… … …… …… ………u…_ …… mm…… …… nwm/u```

```… … mm…—mg…»— sa…… su a.…….… : :mmræwesdæ ; ; m…decnﬂwtülﬂws WWW… un»… M on m…… . … … .. m...… wma: .,… … V, …… … …;………yg…gn… …… pe- le…samemeuuœwpwv m… mum```

Clearly, something went wrong here. When I got that output, I thought that I did not train properly the function. But it was not the answer. As described in that post (in French) it is necessary to have a clean picture, to read it properly And actually, if we zoom in our picture – the first one, used by Fleur, to show me that package – we have while for the second one – with a lower resolution – we have It is necessary to have a scan of a typed text with high resolution… And you have to admit that it is awesome….

The good thing is that I have to work with a judge, in France, to assess quality of experts. And since most of the reports are typed, and then scanned, I am glad to have such a function. I just have to make sure that the resolution is high enough…