Chargement et préparation des données

library(sp)
library(SpatialPosition)

# Téléchargement des données
guetta <- url(description = "http://dansmonlabo.com/files/pays_v3.csv")
ctry <- read.csv(guetta)
# Sélection de certaines colonnes
ctry <- ctry[,c("code_pays_iso_a3","count")]
# Suppression des lignes non renseignées
ctry <- ctry[!is.na(ctry$count),]
# Sélection des pays cités plus d'une fois
ctry <- ctry[ctry$count > 0,]
head(ctry)
##   code_pays_iso_a3 count
## 1              AFG    12
## 2              ZAF     3
## 3              ALB     1
## 4              DZA    22
## 5              DEU    95
## 9              SAU    30
# Import du fond de carte des pays du monde, des capitales
load(url('http://wukan.ums-riate.fr/bernardguetto/worlddata.RData'))
par (mar = c(0,0,0,0))
plot(graticule, col = "grey90", border = NA)
plot(world, col = "grey60", border = 'grey70', add= T)
plot(cap, pch = 20, col = "#920000", add = T)

Calcul des potentiels à partir des centroïdes des pays

# Collage du nombre de citations dans le df de world
world@data <- data.frame(world@data, 
                         count = ctry[match(x = world@data$ISO3, 
                                            table = ctry$code_pays_iso_a3),2])
# Création d'une grille de points
gridpts <- CreateGrid(w = world,resolution = 200000)

# Selection des points dans les pays (pour alleger les calculs)
gridptsctry <- gridPart <- gridpts[!is.na(over(gridpts, world)$ISO3),]

# Création de la matrice de distances entre centroides de pays
mat <- CreateDistMatrix(knownpts = world, unknownpts = gridptsctry)

# Calcul des potentiels 
pot <- stewart(knownpts = world,
               unknownpts = gridptsctry, 
               matdist = mat, 
               varname = 'count', 
               typefct = "exponential", 
               span = 750000,
               beta = 3)

# Ajouts valeurs des potentiels à la grille complète
gridpts@data <- data.frame(gridpts@data, 
                           OUTPUT = pot@data[match(x = gridpts@data$ID, 
                                                   table = pot@data$ID ),
                                             "OUTPUT"])
gridpts@data[is.na(gridpts@data$OUTPUT),"OUTPUT"] <- 0

# Transformation de la grille de points en raster
rasterpot <- rasterStewart(x = gridpts)

# Création d'un SpatialPolygonDataFrame des contours
bks <- c(10,25,50,100,150,200,250,300,350,400)
contourpot <- contourStewart(x = rasterpot, breaks = bks, type = "poly")
## rgeos version: 0.3-8, (SVN revision 460)
##  GEOS runtime version: 3.4.2-CAPI-1.8.2 r3921 
##  Polygon checking: TRUE
# Gestion des couleurs
cols <-c("#B6EFB6", "#A5D8B1", "#95C2AC", 
         "#85ACA7", "#7190A1", "#5B739B", 
         "#435490", "#253175", "#080E5B")
contourpot@data$cols <- contourpot@data$levels
levels(contourpot@data$cols)[order(as.numeric(levels(contourpot@data$cols)), 
                                   decreasing = F)] <- cols

# Affichage de la carte
# png('centroids.png', width = 715, height = 350)
par (mar = c(0,0,0,0))
plot(graticule, col = "grey90", border = NA)
plot(world, col = "grey60", border = 'grey70', add= T)
plot(contourpot, add=T,border = "grey50", 
     col = as.character(contourpot@data$cols) , lwd= 0.25)
mtext("Quelles zones évoque B. Guetta?", side = 3, -1, cex = 1.2)
mtext('Carte : T. Giraud, 2015, Données : Y. Guéguan', side = 1, line = -1, cex = 0.8)