Category Archives: Maps

Interactive Maps for John Snow’s Cholera Data

This week, in Istanbul, for the second training on data science, we’ve been discussing classification and regression models, but also visualisation. Including maps. And we did have a brief introduction to the  leaflet package,

devtools::install_github("rstudio/leaflet")
require(leaflet)

To see what can be done with that package, we will use one more time the John Snow’s cholera dataset, discussed in previous posts (one to get a visualisation on a google map background, and the second one on an openstreetmap background),

library(sp)
library(rgdal)
library(maptools)
setwd("/cholera/")
deaths <- readShapePoints("Cholera_Deaths")
df_deaths <- data.frame(deaths@coords)
coordinates(df_deaths)=~coords.x1+coords.x2
proj4string(df_deaths)=CRS("+init=epsg:27700") 
df_deaths = spTransform(df_deaths,CRS("+proj=longlat +datum=WGS84"))
df=data.frame(df_deaths@coords)
lng=df$coords.x1
lat=df$coords.x2

Once installed the leaflet package, we can use the package at the RStudio console (which is what we will do here), or within R Markdown documents, and within Shiny applications. But because of restriction we got on this blog (rules of hypotheses.org) So there will be only copies of my screen. But if you run the code, in RStudio you will get interactvive maps in the viewer window.

First step. To load a map, centered initially in London, use

m = leaflet()%>% addTiles() 
m %>% fitBounds(-.141,  51.511, -.133, 51.516)

In the viewer window of RStudio, it is just like on OpenStreetMap, e.g. we can zoom-in, or zoom-out (with the standard + and – in the top left corner)

And we can add additional material, such as the location of the deaths from cholera (since we now have the same coordinate representation system here)

rd=.5
op=.8
clr="blue"
m = leaflet() %>% addTiles()
m %>% addCircles(lng,lat, radius = rd,opacity=op,col=clr)

We can also add some heatmap.

X=cbind(lng,lat)
kde2d <- bkde2D(X, bandwidth=c(bw.ucv(X[,1]),bw.ucv(X[,2])))

But there is no heatmap function (so far) so we have to do it manually,

x=kde2d$x1
y=kde2d$x2
z=kde2d$fhat
CL=contourLines(x , y , z)

We have now a list that contains lists of polygons corresponding to isodensity curves. To visualise of of then, use

m = leaflet() %>% addTiles() 
m %>% addPolygons(CL[[5]]$x,CL[[5]]$y,fillColor = "red", stroke = FALSE)

Of course, we can get at the same time the points and the polygon

m = leaflet() %>% addTiles() 
m %>% addCircles(lng,lat, radius = rd,opacity=op,col=clr) %>%
  addPolygons(CL[[5]]$x,CL[[5]]$y,fillColor = "red", stroke = FALSE)

Continue reading Interactive Maps for John Snow’s Cholera Data

Langue et Cartographie

La plupart des cartes du monde standard de R sont en anglais. L’autre jour, des étudiants souhaitaient visualiser des données tirées d’une base où les noms des pays sont en anglais. Pour obtenir une correspondance entre des noms anglais et des noms français, on peut utiliser la base suivante

> library(gdata)
> library(xlsx)
> download.file("http://www.stat.gouv.qc.ca/statistiques/divisions-territoriales/pays-liste-isq-web.xls","corresp")
>  xls_corresp <- read.xls("corresp",sheet=1,encoding="latin1")

On a ici

>  df_corresp <- data.frame(
+ FR=xls_corresp$X.5,
+ EN=xls_corresp$X.11)
> df_corresp[5:10,]
                    FR                 EN
5  Belgique-Luxembourg Belgium-Luxembourg
6    Îles du Pacifique    Pacific Islands
7          Afghanistan        Afghanistan
8       Afrique du Sud       South Africa
9           Îles Åland      Åland Islands
10             Albanie            Albania

Pour avoir une correspondance, entre les noms sous R, et ceux dans la base à notre disposition, il faut manipuler un peu les chaînes de caractères,

>  df_corresp$FR = as.character(df_corresp$FR)
>  df_corresp$FR = iconv(df_corresp$FR, to="ASCII//TRANSLIT") 
>  df_corresp$FR = tolower(df_corresp$FR)
>  remove_minus = function(s) paste(unlist(strsplit(s, split='-',fixed=TRUE)),collapse="")
>  remove_space = function(s) paste(unlist(strsplit(s, split=' ',fixed=TRUE)),collapse="")
>  df_corresp$FR = sapply(df_corresp$FR,remove_minus)
>  df_corresp$FR = sapply(df_corresp$FR,remove_space)

> df_corresp$EN = as.character(df_corresp$EN)
> df_corresp$EN = iconv(df_corresp$EN, to="ASCII//TRANSLIT") 
> df_corresp$EN = tolower(df_corresp$EN)
> df_corresp$EN = sapply(df_corresp$EN,remove_minus)
> df_corresp$EN = sapply(df_corresp$EN,remove_space)
> split_dots = function(s) strsplit(s, split=':',fixed=TRUE)[[1]][1]

Si on regarde les pays que l’on a pu convertir le nom en anglais pour avoir une correspondance avec le nom du pays dans la base de R,

> library(maps)
>  world<-map(database="world")
>  world$pays_EN <- world$names  
>  world$pays_EN <- tolower(world$pays_EN)
>  world$pays_EN = sapply(world$pays_EN,remove_space) 
>  world$pays_EN = sapply(world$pays_EN,remove_minus) 
>  world$pays_EN = sapply(world$pays_EN,split_dots) 
>  world$pays_FR <- df_corresp$FR[match(world$pays_EN, df_corresp$EN)]

on obtient le graphique suivant

>  color <- !is.na(world$pays_FR)
>  map(database="world", fill=TRUE, col=color)

Les seuls pays pour lesquels on n’a pas de correspondance sont les États-Unis d’Amérique (usa dans la base de R), la Russie (ussr), les Congos (avec la République Démocratique et l’autre), et la Côte d’Ivoire. En bricolant un peu sur ces 4 pays, on pourra avoir une correspondance entre les noms utilisés sous R, et les noms en français.

John Snow, and Google Maps

In my previous post, I discussed how to use OpenStreetMaps (and standard plotting functions of R) to visualize John Snow’s dataset. But it is also possible to use Google Maps (and ggplot2 types of graphs).

library(ggmap)
get_london <- get_map(c(-.137,51.513), zoom=17)
london <- ggmap(get_london)

Again, the tricky part comes from the fact that the coordinate representation system, here, is not the same as the one used on Robin Wilson’s blog.

> library(maptools)
> setwd("/cholera/")
> deaths <- readShapePoints("Cholera_Deaths")
> head(deaths@coords)
coords.x1 coords.x2
0  529308.7  181031.4
1  529312.2  181025.2
2  529314.4  181020.3
3  529317.4  181014.3
4  529320.7  181007.9
5  529336.7  181006.0
> X <- deaths@coords

or, use d X_deaths.RData. So now, we have to change it

df_deaths <- data.frame(X)
library(sp)
library(rgdal)
coordinates(df_deaths)=~coords.x1+coords.x2
proj4string(df_deaths)=CRS("+init=epsg:27700") 
df_deaths = spTransform(df_deaths,CRS("+proj=longlat +datum=WGS84"))

Here, we have the same coordinate system as the one used in Google Maps. Now, we can add a layer, with the points,

london + geom_point(aes(x=coords.x1, y=coords.x2),data=data.frame(df_deaths@coords),col="red")

Again, it is possible to add the density, as an additional layer,

london + geom_point(aes(x=coords.x1, y=coords.x2), 
data=data.frame(df_deaths@coords),col="red")+
geom_density2d(data = data.frame(df_deaths@coords), 
aes(x = coords.x1, y=coords.x2), size = 0.3) + 
stat_density2d(data = data.frame(df_deaths@coords), 
aes(x = coords.x1, y=coords.x2,fill = ..level.., alpha = ..level..),size = 0.01, bins = 16, geom = "polygon") + scale_fill_gradient(low = "green", high = "red",guide = FALSE) + 
scale_alpha(range = c(0, 0.3), guide = FALSE)

 

John Snow, and OpenStreetMap

While I was working for a training on data visualization, I wanted to get a nice visual for John Snow’s cholera dataset. This dataset can actually be found in a great package of famous historical datasets.

library(HistData)
data(Snow.deaths)
data(Snow.streets)

One can easily visualize the deaths, on a simplified map, with the streets (here simple grey segments, see Vincent Arel-Bundock’s post)

plot(Snow.deaths[,c("x","y")], col="red", pch=19, cex=.7,xlab="", ylab="", xlim=c(3,20), ylim=c(3,20))
slist <- split(Snow.streets[,c("x","y")],as.factor(Snow.streets[,"street"]))
invisible(lapply(slist, lines, col="grey"))

Continue reading John Snow, and OpenStreetMap

Subjective Ways of Cutting a Continuous Variables

You have probably seen @coulmont‘s maps. If you haven’t, you should probably go and spend some time on his blog (but please, come back afterwards, I still have a story to tell you). Consider for instance the maps we obtained for a post published in Monkey Cage, a few months ago,

The codes were discussed on a blog post (I spent some time on the econometric model, not really on the map, by that time).

My mentor in cartography, Reka (aka @visionscarto) taught me that maps were always subjective. And indeed.

Continue reading Subjective Ways of Cutting a Continuous Variables

Shapefiles from Isodensity Curves

Recently, with @3wen, we wanted to play with isodensity curves. The problem is that it is difficult to get – numerically – the equation of the contour (even if we can easily plot it). Consider the following surface (just for fun, in order to illustrate the idea)

> f=function(x,y) x*y+(1-x)*(1-y)
> u=v=seq(0,1,length=21)
> v=seq(0,1,length=11)
> f=outer(u,v,f)
> persp(u,v,f,theta=angle,phi=10,box=TRUE,
+ shade=TRUE,ticktype="detailed",xlab="",
+ ylab="",zlab="",col="yellow")

For instance, assume that we want to locate areas where the density exceed 0.7 (here in the lower left corner, SW, and the upper right corner, NE)

> image(u,v,f)
> contour(u,v,f,add=TRUE,levels=.7)

Continue reading Shapefiles from Isodensity Curves

Kernel Density Estimation with Ripley’s Circumferential Correction

The revised version of the paper Kernel Density Estimation with Ripley’s Circumferential Correction is now online, on hal.archives-ouvertes.fr/.

In this paper, we investigate (and extend) Ripley’s circumference method to correct bias of density estimation of edges (or frontiers) of regions. The idea of the method was theoretical and difficult to implement. We provide a simple technique — based of properties of Gaussian kernels — to efficiently compute weights to correct border bias on frontiers of the region of interest, with an automatic selection of an optimal radius for the method. We illustrate the use of that technique to visualize hot spots of car accidents and campsite locations, as well as location of bike thefts.

There are new applications, and new graphs, too

Most of the codes can be found on https://github.com/ripleyCorr/Kernel_density_ripley (as well as datasets).

Generating Hurricanes with a Markov Spatial Process

The National Hurricane Center (NHC) collects datasets with all  storms in North Atlantic, the North Atlantic Hurricane Database (HURDAT). For all sorms, we have the location of the storm, every six jours (at midnight, six a.m., noon and six p.m.). Note that we have also the date, the maximal wind speed – on a 6 hour window – and the pressure in the eye of the storm.

It is possible to run the following function

library(XML)
extract.track=function(year=2012,p=TRUE){

Continue reading Generating Hurricanes with a Markov Spatial Process

Crowded Cities, Paris, Hong Kong and Montréal

Over the past years, I’ve been living in different cities, all of them being completely different, compared with the others. I have been living in Paris, which is a big city in Europe, with a large suburban neighborhood, too (la banlieue).

Then, I’ve been living in Hong Kong, which is a larger city, in Asia.

It was crowded. I mean, it was the feeling I had, while I was living there. And more recently, I’ve been living in Montréal, in North America. Montreal is a large city. Or to be more specific, an island,

The three cities are quite different. Paris, 2.211 million unhabitants, and 105,4 km² (density 21,057 unhabitants per km²). Montréal, 1.621 million unhabitants, and three times wider 365.1 km² (density 4,441 unhabitants per km²). Hong Kong, 7.234 million unhabitants, and again three times wider 1,104 km² (density 6,553 unhabitants per km²). In Hong Kong, there are several hill where it is not possible to build anything: on a large part of the island, the density is null.

Continue reading Crowded Cities, Paris, Hong Kong and Montréal

On retourne au pôle nord ?

Depuis plusieurs jours, j’essaye de me préparer autant que possible au débat de mardi soir, et plus j’avance dans mes lectures, plus je trouve que la question “à qui appartient le pôle nord” est complexe. Un autre exemple récent de “revendication” d’une région par d’autres pourrait être celui de l’Ukraine. Pour simplifier outrageusement un débat (là aussi), j’ai l’impression qu’une partie de l’Ukraine se sent européenne, et une autre se sent russe. Un outil pour visualiser cette appartenance peut être d’utiliser des données de langue, par exemple.

File:Ukraine cencus 2001 Ukrainian.svg

Avec des données et une visualisation sur une carte, on peut tenter d’éclairer un peu. Mais dans le cas du pôle nord, on ne peut pas vraiment utiliser de données de ce genre (pas à ma connaissance, malheureusement). Cela dit, on peut apprendre des choses passionnantes dans un ancien numéro de Courrier International sur cette même question :

Dans mon précédant billet, j’avais parlé de la discussion que j’avais eu, rapidement, avec ma fille. En emmenant mon fils à l’escrime vendredi soir, je lui expliquais la problématique du débat, et sa réponse a été sans ambiguïté: le pôle nord appartient à celui qui a mis son premier un drapeau dessus ! On retrouve là la visualisation graphique de l’affiche du débat. Maintenant je dois avouer que ce genre de revendication me gêne un peu. La symbolique de la pose du drapeau est troublante quand on y pense. Une des images à laquelle on pense immédiatement est d’ailleurs le drapeau américain, posé sur la lune (en écho au drapeau hissé lors de la bataille d’Iwo Jima, 硫黄島, en février et mars 1945). A l’époque du drapeau américain sur la lune, il avait d’ailleurs été précisé « this act is intended as a symbolic gesture of national pride in achievement and is not to be construed as a declaration of national appropriation by claim of sovereignty » dans un texte adopté par le Sénat. Donc peut être que cette réponse est trop rapide ! Je voulais donc revenir sur deux éléments clés du débat : les pays qui pourraient revendiquer le pôle nord, et le pôle nord, en tant que tel.

  • les pays autour du pôle nord

Pour l’instant, j’avais fait quelques cartes pour visualiser le pôle, par rapport aux pays qui l’entourent (essentiellement le Canada, la Russie et le Groenland). J’avais retenu comme contour pour mes cartes les limites terrestres des pays, tels que définis sur les cartes. Mais on peut aller plus loin, comme l’évoquent certains, en tenant compte de la zone de 200 milles maritimes (la notion d’Exclusive Economic Zone).

La (petite) région au centre est celle qui est découpée dans la carte de Courrier International, en début de billet. Comme on le voit sur la carte ci-dessus, en rajoutant ces zones, effectivement, nos trois nations sont très proches du pôle (le contour de glace est celui observé pendant les 12 mois de 2013),

Pour moi, l’avantage de ces zones est que c’est plus simple de regarder l’intersection entre la zone de glace, et ces zones économiques (et de calculer l’aire) que de mesurer la longueur de frontière entre un pays et la zone de glace. On va essayer en fin de billet de regarder quelle proportion de la zone de glace chevauchait une zone exclusive (je passe ici sous silence la difficulté technique de calcul d’une aire d’une région contenant le pôle nord… étrangement, c’est complexe, d’un point de vue computationel. Les plus curieux iront voir le billet dédié).

Cela dit, certains pourraient critiquer l’idée de tenir compte de cette extension (horizontale) du territoire, en tenant compte de ces zones économiques. C’est quoi la prochaine étape ? Une appropriation verticale, afin de prendre possession des nuages. Cela dit, je plaisante à peine car la possession verticale de l’espace existe, avec les espaces aériens… D’ailleurs les russes ont évoqués une appropriation verticale, mais sous terre, car la dorsale de Lomonossov  – Хребет Ломоносова – relirait le pôle directement au territoire russe.

Mais je vais arrêter là ma discussion sur les pays.

  • quel pôle nord ?

J’avais un peu laissé de côté la question, mais je pense qu’il serait quand même intéressant d’en parler au moins 2 minutes. De quoi parle-t-on quand on parle du pôle nord ? Pour l’instant, j’avais utilisé la définition usuelle du pôle nord, correspondant au point le plus septentrional de la terre (ce qui, j’en conviens, est un peu tautologique comme définition). En terme de coordonnées, il s’agit du point à 90° de latitude Nord. Il n’a pas de longitude, car c’est précisément le point au nord où tous les méridiens (et tous les fuseaux horaires) se coupent. J’avais aussi dit, dans mon précédant billet qu’il s’agissait du point de la surface terrestre, de l’hémisphère nord, situé sur l’axe de rotation de la terre. Ce point géographique n’est malheureusement pas fixe à la surface de la Terre, car l’axe de rotation de la terre varie (très faiblement).

Mais il existe un autre pôle, le pôle nord magnétique, qui est l’unique sur la surface de la terre où le champ magnétique (terrestre) pointe vers le centre de la terre (et orthogonal à la surface de la terre, s’il pointe vers le haut, c’est le pôle sud). A priori, on pourrait se dire que ce point n’aura pas grand intérêt dans la discussion, car s’approprier un point mouvant serait surprenant. Car il bouge beaucoup, en témoigne la carte ci-dessous (les coordonnées sont facilement accessibles, sur http://ngdc.noaa.gov/geomag/data/poles/ par exemple)

Un utilisant ces champs, on peut d’ailleurs définir un pôle nord gravitationnel (je me permets d’inventer le mot, par analogie au pôle nord magnétique), à l’aide non pas du champ magnétique, mais du champ gravitationnel. En anglais, on parle de polar wander. Pour faire simple, ce pôle nord sera lié à la répartition de la masse, sur terre, et des mouvements d’eau, en particulier (le reste étant assez stable à côté). Or comme l’ont montré récemment Chen, Wilson, Ries et Tapley, dans Rapid ice melting drives Earth’s pole to the east (article remarqué par Richard Lovett, en mai 2013, dans Nature) le pôle nord semble se déplace plus rapidement que prévu, à cause du réchauffement climatique.

Mais je pense que cette discussion (quoi que potentiellement passionnante) nous éloigne de ce qui devrait être évoqué mardi soir. En fait, mon point serait que, plutôt que de visualiser le pôle nord, en tant que point, il est plus intéressant de visualiser une région. Et plutôt que de considérer une région au delà du cercle polaire (j’avais utilisé des cercles dans mon précédant billet), on peut davantage utiliser la zone glacée, comme le montre la carte

Encore une fois, le soucis est que cette zone glacée varie énormément ! On observe qu’entre la moitié et les deux-tiers de la surface de la glace se situe dans les zones économiques des trois principaux pays limitrophes, la Russie et le Canada.

Et si on y regarde de plus près, on observe d’ailleurs que les surfaces occupées respectivement par la Russie et le Canada sont comparables,

avec toutefois un comportement cyclique, la présence canadienne étant au plus faible à l’automne,

J’ai l’impression que je pourrais passer des heures sur ces cartes. Mais je vais probablement m’arrêter là pour ma préparation du débat. Et même si ces débats ne laissent que peu de temps pour la prise de notes, j’essayerais de prendre un peu de temps pour revenir sur ce qui s’est dit, en milieu de semaine, sur le blog, ou sur Twitter.

Moving the North Pole to the Equator

I am still working with @3wen on visualizations of the North Pole. So far, it was not that difficult to generate maps, but we started to have problems with the ice region in the Arctic. More precisely, it was complicated to compute the area of this region (even if we can easily get a shapefile). Consider the globe,

worldmap <- ggplot() + 
geom_polygon(data = world.df, aes(x = long, y = lat, group = group)) +
scale_y_continuous(breaks = (-2:2) * 30) +
scale_x_continuous(breaks = (-4:4) * 45)

and then, add three points in the northern hemisphere, and plot the associated triangle

P1 <- worldmap + geom_polygon(data = triangle, aes(x = long, y = lat, group = group), 
fill ="blue", alpha = 0.6, col = "light blue", size = .8)+
geom_point(data = triangle, aes(x = long, y = lat, group = group),colour = "red")+

for some given projection, e.g.

coord_map("ortho", orientation=c(61, -74, 0))

This can be done with the following function

proj1=function(x=75){
triangle <- data.frame(long=c(-70,-110,-90*(x<90)+90*(x>90)),
lat=c(60,60,x*(x<90)+(90-(x-90))*(x>90)),group=1, region=1)
worldmap <- ggplot() + 
geom_polygon(data = world.df, aes(x = long, y = lat, group = group)) +
scale_y_continuous(breaks = (-2:2) * 30) +
scale_x_continuous(breaks = (-4:4) * 45)
P1 <- worldmap + geom_polygon(data = triangle, aes(x = long, y = lat, group = group), 
fill ="blue", alpha = 0.6, col = "light blue", size = .8)+
geom_point(data = triangle, aes(x = long, y = lat, group = group),colour = "red")+
coord_map("ortho", orientation=c(61, -74, 0)) 
print(P1)
}

or

I am not sure if I understand why the projection of the triangle is not convex on the graph above, but say it’s not a big deal, here. Actually, our problem is that our interest is on regions (polygons, from a geometrical point of view) that do contain the North Pole. And here, it starts to be messy. I can easily move the upper point on the other side of the globe, but the polygon is not correct,

I do understand that it should be a problem, non-trivial, but it means that it should not be that simple to compute the area of a polygon (a region) that contains the North Pole. Which is exactly what we did observe in our computation. And I believe that one heuristic interpretation is related to the following graph

My skills in geometry are extremely poor. So do not expect that I will go through the code of the function that compute the area of a polygon ! Actually, my idea is the following : if the problem is that the North Pole is in the region, let’s consider some rotation, to shift the North Pole on the Equation. The code here is, from latitudes and longitude, to get new latitudes and longitudes, after a rotation around the y-axis (the North Pole will go down, along Greenwhich meridian) is

rotation=function(Z,theta){
lon=Z[,1]/180*pi; lat=Z[,2]/180*pi
x=cos(lon)*cos(lat)
y=sin(lon)*cos(lat)
z=sin(lat)
pt1=cbind(x,y,z)
M=matrix(c(cos(theta),0,-sin(theta),0,1,0,sin(theta),0,cos(theta)),3,3)
pt2=t(M%*%t(pt1))
lat=asin(pt2[,3])*180/pi
lon=atan2(pt2[,2],pt2[,1])*180/pi
return(cbind(lon,lat))}

With a rotation from  (no change) to  (the North Pole on the equator), we get

From now on, it is possible to compute the area of any region containing the North Pole ! One should simply apply the rotation function on all datebases generated from shapefiles (and then the opposite rotation to get a proper location) ! We can then compute the centroid of the ice region, for example,

r.glace=glace
r.glace[,1:2]=rotation(glace[,1:2],pi/2)
M=matrix(NA,length(unique(glace$id)),3)
j=0
for(i in unique(glace$id)){j=j+1
Polyglace <- as(r.glace[glace$id==i,c("long","lat")],"gpc.poly")
M[j,1]=area.poly(Polyglace)
M[j,2:3]=centroid(r.glace[r.glace$id==i,c("long","lat")])
}
Z=c(weighted.mean(M[,2],M[,1]),weighted.mean(M[,3],M[,1]))
rotation(rbind(Z),-pi/2)[1,])

And we get

and below, we can visualize all the locations of the centroid of the ice region in the past 25 years

We are winter

Non, pas de billet pour annoncer la reprise de Game of Thrones dans les jours à venir. Mais un billet rapide pour rappeler qu’au Canada, l’hiver, c’est un sujet important, et sérieux.

C’est ce thème qu’avait repris le Comité Olympique Canadien avant et pendant les Jeux.

Et lors de mon excursion au Mexique la semaine passée (où nous avons retrouvé des collègues chiliens d’ailleurs), je me suis rendu compte que les Jeux Olympiques d’hiver n’intéressaient que peu de monde, en dehors des Canadiens, et des nations qui subissent l’hiver, comme les Russes, ou les Norvégiens (entre autres). Un autre point comment de toutes ces nations est d’avoir un intérêt très fort pour un endroit important d’un point de vue géopolitique, le pôle nord. Et on m’a demandé d’intervenir dans une conférence de Thierry Garcin, et Philippe Reka. J’ai toujours un plaisir fou à discuter avec Philippe, donc c’est avec plaisir que je viendrais, mais j’avoue (une nouvelle fois) avoir du mal à trouver ma place (oui, je suis un peu le sparadrap du capitaine Haddock)

Il y a tout juste un an, Philippe publiait dans le Monde Diplo un article inspirant sur la Cartographie Radicale,

Depuis des décennies, la cartographie traditionnelle revendique le statut de science exacte s’appuyant sur des données fiables. Elle se targue de fournir une image neutre et fidèle de la réalité. Mais une telle approche fait l’impasse sur l’utilisation politique et sociale de la carte, et sur son rôle tant de propagande que de contestation. Depuis le début des années 2000, et dans le désordre, émerge une pratique cartographique qui se dit « radicale » (on parle aussi de « cartographie critique » ou de « contre-cartographie »), riche combinaison revendiquée d’art, de sciences, de géographie, de politique et de militantisme social. Les cartographes orthodoxes ont vu naître avec beaucoup de méfiance les premiers projets exploratoires qui bousculaient, voire cassaient, les codes. « On peut “crobarder” des schémas réducteurs donnant l’illusion de la synthèse, et en faire des figures de style prétentieuses, écrivait un cartographe professionnel en écho à l’une des toutes premières expositions de cette mouvance. La carte géographique, ce fantastique outil de savoir et d’appréhension de l’espace et du monde, se voit ravalée au rang de simple outil de prestidigitation. »

C’est peu dire que le symposium « Art and Cartography — Cartography and Art », organisé par l’Académie des beaux-arts et l’Université technique de Vienne (Autriche) en février 2008, fut le théâtre d’intenses discussions, où s’exprimèrent des positions inconciliables. La cartographie engagée a trouvé la voie d’une forme de contestation au service de la justice sociale, dénonçant des pratiques politiques et économiques douteuses. Elle est un exercice libre de déconstruction de l’espace et des phénomènes sociaux, pour lequel les protagonistes se permettent de pervertir les conventions les plus classiques. Les premiers à s’en emparer furent des artistes, mais aussi des architectes, des urbanistes et des militants, rejoints plus tard par des géographes. Les projets, individuels ou collectifs, se sont surtout multipliés en Argentine, aux Etats-Unis, en Suisse, en Allemagne et en France, couvrant un champ thématique assez large (finance, surveillance et sécurité, marketing, environnement, organisation des territoires, pour ne donner que quelques exemples), et profitant largement de la puissance des outils de la cartographie participative et des réseaux sociaux.

Déchiffrer pour dénoncer, mettre en images des processus peu visibles qui concourent à confisquer l’espace public (voire des biens publics), à compromettre les libertés individuelles, à détourner des lois, tels sont les principaux objectifs des promoteurs de cette cartographie. Dans cet esprit, informer n’est qu’un début. L’étape suivante de la démarche reste l’action pour le changement sur le terrain. Largement informelle, cette initiative se reconnaît quelques affinités avec les mouvements de résistance tels que celui des « indignés » ou Occupy. Ainsi, le collectif d’artistes argentins Grupo de Arte Callejero ne s’est pas contenté de dénoncer les bourreaux de l’ancienne dictature en cartographiant les lieux de leur paisible retraite : il a aussi organisé des itinéraires « touristiques » qui passaient par leurs domiciles et mené campagne sur des panneaux publicitaires… Le projet « Duty Free Shop », en ce qu’il se propose de déconstruire, rendre lisible et donc fragiliser un processus d’aliénation des citoyens dans un espace public, s’inscrit de plain-pied dans ce mouvement pétillant.

En tant qu’amateur (plus qu’expert) de visualisation de données, je rejoins complètement Philippe dans sa lecture critique (je ne vois les cartes que comme un support graphique parmi beaucoup d’autres). Philippe a ainsi proposé une très jolie visualisation de l’arctique.

Et comme souvent, Philippe a raison. Pour comprendre l’arctique, il faut avoir une visualisation claire. Car comme beaucoup de monde, j’ai été éduqué avec des cartes qui ressemblaient à ça

Autrement dit, pour moi, l’arctique et le pôle nord ont longtemps été des concepts flous, un peu comme l’horizon… c’est au nord, tout là haut. Par contre, j’ai bien conscience que l’antarctique existe, par exemple, mais les cartes usuelles ne nous permettent de prendre conscience de sa superficie (un peu moins de deux fois la taille de l’Australie). D’ailleurs, si on va sur https://www.google.ca/maps/place/North+Pole/ on peut voir (en butant contre la partie supérieure de la carte) qu’on ne voit rien, si on cherche à visualiser ce pôle nord (même s’il est sous la glace)

Comme le notait ma fille, le pôle nord est une enjeu stratégique important car c’est là que se trouve le village du père Noël (c’est manifestement la croyance répandue en Amérique). Et si le père Noël n’existe pas (c’est ce que pense ma fille), peut être que le pôle nord n’existe pas non plus ? Le premier défi est donc d’avoir une visualisation aussi pertinente que possible du pôle. Je vais donc revenir un peu sur les problématique de projection. Car les Soviétiques ne se trompaient pas lorsqu’ils représentèrent l’arctique avec un point de vue qui donnait l’impression que l’arctique était russe,

Pour appréhender les projecteurs, deux ouvrages fabuleux expliquent cela très clairement. Flattening the Earth – Two Thousand Years of Map Projections de John Parr Snyder, mais aussi (et surtout) How to Lie with Maps de Mark Monmonier. Mark Monmonier a également publié un livre passionnant sur le sujet (en particulier le chapitre 6), intitulé Air Apparent – How Meteorologists Learned to Map, Predict, and Dramatize Weather

Pour visualiser le pôle nord, on peut utiliser quelques fonctions simples (basées sur un travail en cours avec @3wen), pour visualiser différentes projections. On va ici mettre l’accent sur trois pays importantes quand on parle du pôle nord, le Canada, la Russie, et le Groenland (qui est techniquement un pays consultatif), et sur trois deux régions, une située à moins de 1500km du pôle nord, et l’autre à 3000km du pôle (soit un peu en dessous du cercle arctique)

Si on veut changer de point de vue, on peut visualiser l’animation suivante

D’ailleurs, si on compare la proportion du disque occupé par chacun des pays (j’ai pris ici chacune des huit pays de l’arctique), on obtient, avec en abscisse la distance au pôle, et en ordonnée, la proportion du disque occupée par chacun des pays, avec le Canada (rouge), la Russie (bleu) et le Groenland (vert)

Techniquement, le pôle pose de réels soucis car la plupart des données sont représentées par leur latitude, et leur longitude. Les pôles sont les seuls points – dans cette représentation – à avoir une infinité de coordonnées !

File:Latitude and Longitude of the Earth.svg

La latitude est de 90°, mais toutes les longitudes sont valides. On a ainsi des soucis y compris sur les cartes construites à partir de données satellites, qui pourtant pourraient balayer uniformément la sphère.

Le pôle nord peut ainsi compter beaucoup de discontinuités, uniquement à cause de la représentation des données par la latitude et la longitude : en quelques kilomètres, on passe d’une zone rouge à une zone blanche, ce qui pose des soucis de lissage de la carte (les données n’étant manifestement pas considérées ici comme circulaire pour le lissage).

Maintenant, plus qu’un point sur la carte, on peut voir le pôle nord comme une région, et non pas comme un point sur la sphère, défini très arbitrairement. Bon, j’exagère un peu ici, car le pôle nord est défini par l’axe de rotation de la terre,

qui va induire une espère de symétrie, par exemple sur le climat. Mais le second problème, selon moi, qui va alors apparaître est que le pôle nord vu comme un territoire (un territoire de glace) est alors un territoire qui va changer avec les saisons.

Par exemple, en septembre dernier, la glace occupait la région suivante, ce qui montre que peu de pays peuvent effectivement rejoindre le pôle nord ‘à pieds’ quel que soit le moment de l’année.

On continue de travailler avec @3wen sur le sujet, mais le sujet et vaste, ambitieux… et complexe ! à suivre donc…

Equidistant points on a map

This morning, I had a comment on a recent post, regarding a graph I did upload on the blog, which was extracted from a paper now online (see http://hal.archives-ouvertes.fr/hal-00871883). Jo (from KUL, I guess I can share that piece of information) asked me

I was wondering whether you would want to share the R code for plotting figures 1 and 14? W.r.t. the former, the figure-in-figure is a nice touch; as to the latter, I am curious to know how you translated distance in km to the size parameters of the graph (par(“usr”)) for plotting the corresponding concentric circles (and the arrow indicating the radius) on top of your map.

At first, I thought I made a mistake in my plot. I mean, each time I have a question, I start to be suspicious, and I start to wonder if what I did was valid, or not. Here was the graph

Let’s make it clear: I do not draw circles here. So yes, I believe that what I did is valid. What I did is simple. First, I get the background map,

library(maps)
map("world",xlim=c(130,150),ylim=c(25,45),fill=TRUE,col="light green")

Then, I need some function to compute distance from coordinates. The functions I use are

deg2rad = function(deg) return(deg*pi/180)
DISTANCEDEG = function(long1, lat1, long2, lat2) {
R=6371; d=acos(sin(lat1)*sin(lat2) + cos(lat1)*cos(lat2) * cos(long2-long1)) * R
return(d) 
}

The center here will be Tokyo (東京),

X=139+45/60
Y=35+40/60

The idea now is simple: I generate a grid (here 501×501)

VX=seq(X-10,X+10,length=501)
VY=seq(Y-10,Y+10,length=501)
VtX=rep(VX,each=length(VY))
VtY=rep(VY,length(VX))
ZDeg=deg2rad(cbind(VtX,VtY))

I compute the distance from all those points to Tokyo, and check is the distance is larger or smaller than a given value,

L=500
D1=DISTANCEDEG(ZDeg[,1],ZDeg[,2],deg2rad(X),deg2rad(Y))<L

If the distance is smaller than 500km, then I put a blue dot on the graph,

points(VtX[D1],VtY[D1],pch=19,cex=.2,col="light blue")

Then I use the same procedure for 250km (obviously, it is more convenient to start from larger and to go to smaller distances)

L=250
D=DISTANCEDEG(ZDeg[,1],ZDeg[,2],deg2rad(X),deg2rad(Y))<L
points(VtX[D],VtY[D],pch=19,cex=.2,col="light yellow")

Then, I did draw an arrow to ilustrate the largest distance

k=which.max(VtX[D1])
arrows((VtX[D1])[k],(VtY[D1])[k],X,Y,code=3,length=.1)
text(((VtX[D1])[k]+X)/2,Y+.35,"500 km")

And now, I have the graph.

Now, the point is that it should depend on the kind of projection we use, right? So here is a function that can be used for different kind of projections (some slight changes are necessary, since the map is now centered on some point, and we cannot use standard coordinates)

library(mapproj)
mapjapan = function(pr="conic",pm=45){
map("world","japan",fill=TRUE,col="light green",projection=pr, par=pm)
MP=mapproject(data.frame(x=X,y=Y),projection="")
Xp=MP$x
Yp=MP$y
VX=seq(X-10,X+10,length=501)
VY=seq(Y-10,Y+10,length=501)
VtX=rep(VX,each=length(VY))
VtY=rep(VY,length(VX))
MP=mapproject(data.frame(x=VtX,y=VtY),projection="")
VtXp=MP$x
VtYp=MP$y
ZDeg=deg2rad(cbind(VtX,VtY))
L=500
D1=DISTANCEDEG(ZDeg[,1],ZDeg[,2],deg2rad(X),deg2rad(Y))<L
points(VtXp[D1],VtYp[D1],pch=19,cex=.2,col="light blue")
L=250
D=DISTANCEDEG(ZDeg[,1],ZDeg[,2],deg2rad(X),deg2rad(Y))<L
points(VtXp[D],VtYp[D],pch=19,cex=.2,col="light yellow")
L=100
D=DISTANCEDEG(ZDeg[,1],ZDeg[,2],deg2rad(X),deg2rad(Y))<L
points(VtXp[D],VtYp[D],pch=19,cex=.2,col="light blue")
L=50
D=DISTANCEDEG(ZDeg[,1],ZDeg[,2],deg2rad(X),deg2rad(Y))<L
points(VtXp[D],VtYp[D],pch=19,cex=.2,col="light yellow")
points(Xp,Yp,pch=19,cex=.4,col="red")
map("world","japan",projection=pr, par=pm,add=TRUE)
}

The default function here produces a map based on a conic projection,

mapjapan()

But we can also use a Bonne projection (a pseudo-conic one, named after Rigobert Bonne)

mapjapan("bonne")

or Lagrange projection,

mappjapan("lagrange",NULL)

and (as a last one), Albers projections,

mapjapan("albers",c(30,40))

Of course, much more projections are possible !

We do not see much here, right ? So let us play with a larger country to visualize something. Like Canada. And the distance to, say, Winnipeg.

The first thing to do is to define the coordinates of Winnipeg,

X=-(97+08/60)
Y=(49+53/60)

Then, we slightly change our function

mapcanada = function(pr="conic",pm=45){
map("world","canada",fill=TRUE,col="light green",projection=pr, par=pm)
MP=mapproject(data.frame(x=X,y=Y),projection="")
Xp=MP$x
Yp=MP$y
VX=seq(X-30,X+30,length=501)
VY=seq(Y-30,Y+30,length=501)
VtX=rep(VX,each=length(VY))
VtY=rep(VY,length(VX))
MP=mapproject(data.frame(x=VtX,y=VtY),projection="")
VtXp=MP$x
VtYp=MP$y
ZDeg=deg2rad(cbind(VtX,VtY))
L=2000
D1=DISTANCEDEG(ZDeg[,1],ZDeg[,2],deg2rad(X),deg2rad(Y))<L
points(VtXp[D1],VtYp[D1],pch=19,cex=.2,col="light blue")
L=1000
D=DISTANCEDEG(ZDeg[,1],ZDeg[,2],deg2rad(X),deg2rad(Y))<L
points(VtXp[D],VtYp[D],pch=19,cex=.2,col="light yellow")
L=500
D=DISTANCEDEG(ZDeg[,1],ZDeg[,2],deg2rad(X),deg2rad(Y))<L
points(VtXp[D],VtYp[D],pch=19,cex=.2,col="light blue")
L=200
D=DISTANCEDEG(ZDeg[,1],ZDeg[,2],deg2rad(X),deg2rad(Y))<L
points(VtXp[D],VtYp[D],pch=19,cex=.2,col="light yellow")
points(Xp,Yp,pch=19,cex=.4,col="red")
map("world","canada",projection=pr, par=pm,add=TRUE)
}

Now, we can have some fun

mapcanada()

mapcanada("bonne",45)

mapcanada("albers",c(30,40))

mapcanada("lagrange",NULL)

Fun, isn’t it ? Changing the projection will change the shape of equidistant curves.

Visualizing densities of spatial processes

We recently uploaded on http://hal.archives-ouvertes.fr/hal-00725090 a revised version of our work, with Ewen Gallic (a.k.a. @3wen) on Visualizing spatial processes using Ripley’s correction: an application to bodily-injury car accident location

In this paper, we investigate (and extend) Ripley’s circumference method to correct bias of density estimation of edges (or frontiers) of regions. The idea of the method was theoretical and di#cult to implement. We provide a simple technique – based of properties of Gaussian kernels – to compute e#efficiently weights to correct border bias on frontiers of the region of interest, with an automatic selection of an optimal radius for the method. An illustration on location of bodily-injury car accident (and hot spots) in the western part of France is discussed, where a lot of accident occur close to large cities, next to the sea.

Sketches of the R code can be found in the paper, to produce maps, an to describe the impact of our boundary correction. For instance, in Finistère, the distribution of car accident is the following (with a standard kernel on the left, and with correction on the right), with 186 claims (involving bodily injury)

and in Morbihan with 180 claims, observed in a specific year (2008 as far as I remember),

The code is the same as the one mentioned last year, except perhaps plotting functions. First, one needs to defi
ne a color scale and associated breaks

breaks <- seq( min( result $ZNA , na.rm = TRUE ) * 0.95 , max ( result$ZNA , na.rm = TRUE ) * 1.05 , length = 21)
col <- rev( heat . colors (20) )

to
finally plot the estimation

image . plot ( result $X, result $Y, result $ZNA , xlim = range (pol[,
1]) , ylim = range (pol[, 2]) , breaks = breaks , col = col ,
xlab = "", ylab = "", xaxt = "n", yaxt = "n", bty = "n",
zlim = range ( breaks ), horizontal = TRUE )

It is possible to add a contour, the observations, and the border of the polygon

contour ( result $X, result $Y, result $ZNA , add = TRUE , col = "grey ")
points (X[, 1], X[, 2], pch = 19, cex = 0.2 , col = " dodger blue")
polygon (pol , lwd = 2)

Now, if one wants to improve the aesthetics of the map, by adding a Google Maps base map, the
first thing to do – after loading ggmap package – is to get the base map

theMap <- get_map( location =c( left =min (pol [ ,1]) , bottom =min (pol[ ,2]) , right =max (pol [ ,1]) , 
top =max (pol [ ,2])), source =" google ", messaging =F, color ="bw")

Of course, data need to be put in the right format

getMelt <- function ( smoothed ){
res <- melt ( smoothed $ZNA)
res [ ,1] <- smoothed $X[res [ ,1]]
res [ ,2] <- smoothed $Y[res [ ,2]]
names (res) <- list ("X","Y","ZNA")
return (res )
}
smCont <- getMelt ( result )

Breaks and labels should be prepared

theLabels <- round (breaks ,2)
indLabels <- floor (seq (1, length ( theLabels ),length .out =5)) 
indLabels [ length ( indLabels )] <- length ( theLabels ) 
theLabels <- as. character ( theLabels [ indLabels ])
theLabels [ theLabels =="0"] <- " 0.00 "

Now, the map can be built

P <- ggmap ( theMap )
P <- P + geom _ point (aes(x=X, y=Y, col=ZNA), alpha =.3 , data =
smCont [!is.na( smCont $ZNA ) ,], na.rm= TRUE )

It is possible to add a contour

P <- P + geom _ contour ( data = smCont [!is.na( smCont $ZNA) ,] ,aes(x=
X, y=Y, z=ZNA ), alpha =0.5 , colour =" white ")

and colors need to be updated

P <- P + scale _ colour _ gradient ( name ="", low=" yellow ", high ="
red", breaks = breaks [ indLabels ], limits = range ( breaks ),
labels = theLabels )

To remove the axis legends and labels, the theme should be updated

P <- P + theme ( panel . grid . minor = element _ line ( colour =NA), panel
. grid . minor = element _ line ( colour =NA), panel . background =
element _ rect ( fill =NA , colour =NA), axis . text .x= element _ blank() ,
axis . text .y= element _ blank () , axis . ticks .x= element _ blank() ,
axis . ticks .y= element _ blank () , axis . title = element _ blank() , rect = element _ blank ())

The
final step, in order to draw the border of the polygon

polDF <- data . frame (pol)
colnames ( polDF ) <- list ("lon","lat")
(P <- P + geom _ polygon ( data =polDF , mapping =( aes(x=lon , y=lat)), colour =" black ", fill =NA))

Then, we’ve applied that methodology to estimate the road network density in those two regions, in order to understand if high intensity means that it is a dangerous area, or if it simply because there is a lot of traffic (more traffic, more accident),

We have been using the dataset obtained from the Geofabrik website which provides
Open-StreetMap data. Each observation is a section of a road, and contains a few points identifi
ed by their geographical coordinates that allow to draw lines. We have use those points to estimate a proxy of road intensity, with weight going from 10 (highways) to 1 (service roads).

splitroad <- function ( listroad , h = 0.0025) {
pts = NULL
weights <- types . weights [ match ( unique ( listroad $ type ), types .
weights $ type ), " weight "]
for (i in 1:( length ( listroad ) - 1)) {
d = diag (as. matrix ( dist ( listroad [[i]]))[, 2: nrow ( listroad
[[i ]]) ]))
}}
return (pts )
}

See Ewen’s blog for more details on the code, http://editerna.free.fr/blog/…. Note that Ewen did publish a poster of the paper (in French), for the http://r2013-lyon.sciencesconf.org/ conference, that will be organized in Lyon on June 27th-28th, see

All comments are welcome…