Cartographie avec R

Niveau 3 - Des cartes plus riches en information (reliefs, hydrographie)

Cette page est en fabrication. Merci pour votre compréhension.

1- Afficher l'altitude

Charger une image dont chaque pixel correspond à une valeur (un raster) avec la librairie raster

# Exemple - la raster des données altimétriques de Grèce

library(raster) # install.packages("raster") avant si nécessaire

# srtm : meilleure résolution que alt

srtm <- getData('SRTM', lon=23, lat=36)

# Affichage

plot(srtm)

# Décrire un raster

minValue(srtm) ; maxValue(srtm) ; srtm

Afficher une région de la carte des altitudes

# ETAPE 1 - Chargement des données administratives

library(raster)

# Contours administratifs

adm_gr <- getData('GADM', country='GRC', level=1)

# La ligne qui suit permet de limiter la zone traitée à une région,

#ici le Péloponnèse (La carte sera plus rapideà tracer)

adm <- adm_gr[match(toupper(adm_gr$NAME_1[7]),toupper(adm_gr$NAME_1)),]

plot(adm,xlim=c(23.05,23.20),ylim=c(36.422,36.53))

Les contours administratifs de la région étudiée sont affichés

# ETAPE 2 - Chargement des données altimétriques et affichage

# Chargement des données altimétriques

srtm <- getData('SRTM', lon=23, lat=36) # srtm : meilleure résolution que alt

plot(adm,xlim=c(23.05,23.20),ylim=c(36.422,36.53))

plot(srtm ,add=T)

A ce stade, on constate une résolution du raster insuffisante. Il va falloir circonscrire la région d'intérêt pour mieux afficher les données du raster.

Remarque : parfois le chargement automatique avec getData ne fonctionne pas pour les données SRTM. Il faut alors les télécharger soi-même sur le site : http://srtm.csi.cgiar.org/download.

On pourra alors charger le tif (l'image avec la fonction raster)

srtm <- raster("image_du_site_srtm.tif")

Le raster de tout le pays aura un affichage insuffisant qui justifie d'aller plus loin.

Remarque importante : je peux avoir besoin de charger plusieurs rasters d'altitude pour couvrir une grande surface. J'aurais alors besoin de fusionner mes raster de la façon suivante :

srtm1 <- getData('SRTM', lon=23, lat=36) # raster sud

srtm2 <- getData('SRTM', lon=23, lat=48) # raster nord

new_data <- raster::merge(srtm1 , srtm2)

plot(new_data)

# ETAPE 3 - Chargement des données altimétriques locales de la région étudiée

vatikasrtm <- crop(srtm,extent(c(xlim=c(23,23.3),ylim=c(36.422,36.54))))

# Affichage

plot(adm,xlim=c(23.05,23.20),ylim=c(36.422,36.53))

plot(vatikasrtm ,add=T)

plot(adm,add=T)

# On trace le polygone qui nous intéresse pour les étapes suivantes

polygon(c(23.13,23.13,23.15,23.15),

c(36.48,36.5,36.5,36.48),border="red",lwd=2)

Le raster est pleinement exploité mais si on zoome sur le rectangle rouge, on verra ci-après que sa résolution risque de se montrer insuffisante pour épouser les contours de la carte.

Par défaut, le raster risque d'être pixelisé, pas assez de données altimétriques par rapport à la finesse du contours. Certains pixels n'épousent pas les bords de la carte.

Augmenter artificiellement le nombre de pixels d'un raster en subdivisant chaque pixel

# Chargement local des données d'altitude de la région à afficher

vatikasrtm <- crop(srtm,extent(c(xlim=c(23,23.3),ylim=c(36.422,36.54))))

# Les régions vidus du raster sont remplacées par la valeur min : étape importante

vatikasrtm[is.na(vatikasrtm[])] <- minValue(vatikasrtm)

# Subdivision de chaque pixel par 18 : un pixel donne 18² pixels.

# ATTENTION : cette étape peut nécessiter plusieurs minutes

vatikasrtm<- disaggregate(vatikasrtm, fact=18)

# Applique un masque qui va effacer les pixels qui dépassent des contours administratifs

vatikasrtm <- mask(vatikasrtm ,mask= adm, inverse=F)

plot(adm,xlim=c(23.13,23.15),ylim=c(36.48,36.5))

plot(vatikasrtm ,add=T,legend=F)

plot(adm,add=T)

On voit ainsi qu'en fractionnant les pixels du raster, on a pu ajuster au mieux la carte aux contours administratifs. Toutefois, la carte reste pixelisée. Il faudrait pouvoir lister cette structure en marche d'escalier sans perdre de signal.

Lisser un raster pour éviter les phénomènes de pixelisation

# En reprenant le srtm précédent, on limite le raster d'altitude

vatikasrtm <- crop(srtm,extent(c(xlim=c(23,23.3),ylim=c(36.422,36.54))))

# Les pixels vides se voient imposés la valeur min d'altitude locale

vatikasrtm[is.na(vatikasrtm[])] <- minValue(vatikasrtm)

# Les pixels du raster sont fractionnés en 17² (on peut faire moins, plus rapide mais résolution moindre)

vatikasrtm_mask_brut <- disaggregate(vatikasrtm, fact=17)

# On crée une fenêtre de 17 sur 17 de valeurs 1 qui servira à faire une moyenne pour chaque pixel à partir de 17² valeurs

FenetreLissage=matrix(1,17,17) # Création de la fenêtre de lissage

# Lissage avec filtre mean ou median (fun)médian

vatikasrtm_mask<- focal(vatikasrtm_mask_brut,fun=mean,w=FenetreLissage,pad=T,padValue=0)

# On impose les contours administratifs comme masque pour éliminer les pixels hors de la carte

vatikasrtm_mask <- mask(vatikasrtm_mask ,mask= adm, inverse=F); #plot(vatikasrtm_mask)

# Affichage

plot(adm,xlim=c(23.05,23.20),ylim=c(36.422,36.53))

plot(vatikasrtm_mask,add=T)

plot(adm,add=T)

Cette fois-ci la carte ne présente pas de pixelisation. Elle pourra donc être imprimée sans problème.

Changer la couleur/légende de hiérarchisation du raster

Coloration par défaut

plot(vatikasrtm,

legend=F, axes=F,

xlab="", ylab="")

Exemple 1 de coloration

liste_colors <- c("#BCF5A9","#36a013","#E58544","#C6753F","#B6652F")

Exemple 2 de coloration

liste_colors <- c("white","red")

# Exemple 1 et 2

# Changer les valeurs de couleurs de liste_colors pour obtenir l'exemple 2

liste_colors <- c("#BCF5A9","#36a013","#E58544","#C6753F","#B6652F")

colfunc<-colorRampPalette(liste_colors,interpolate="spline")

colors <- colfunc(1000)

plot(adm,xlim=c(23.05,23.20),ylim=c(36.422,36.53))

plot(vatikasrtm_mask,add=T,col=colors,axes=F,xlab="",ylab="",legend=F)

plot(adm,add=T)

2- Afficher les ombrages des pentes

Par convention, l'ombrage doit venir comme si le soleil éclairait du Nord-Ouest, avec un angle de 45° environ.

# Chargement des données altimétriques

srtm <- getData('SRTM', lon=23, lat=36) # srtm : meilleure résolution que alt

# En reprenant le srtm précédent, on limite le raster d'altitude

vatikasrtm <- crop(srtm,extent(c(xlim=c(23,23.3),ylim=c(36.422,36.54))))

# Les pixels vides se voient imposés la valeur min d'altitude locale

vatikasrtm[is.na(vatikasrtm[])] <- minValue(vatikasrtm)

# Les pixels du raster sont fractionnés en 17² (on peut faire moins, plus rapide mais résolution moindre)

# On augmente ainsi la résolution de la carte

vatikasrtm_mask_brut <- disaggregate(vatikasrtm, fact=17)

# On crée une fenêtre de 17 sur 17 de valeurs 1 qui servira à faire une moyenne pour chaque pixel à partir de 17² valeurs

FenetreLissage=matrix(1,17,17) # Création de la fenêtre de lissage

# Lissage avec filtre mean ou median (fun)médian

# Ce lissage va éviter une pixelisation

vatikasrtm_mask<- focal(vatikasrtm_mask_brut,fun=mean,w=FenetreLissage,pad=T,padValue=0)

# On impose les contours administratifs comme masque pour éliminer les pixels hors de la carte

vatikasrtm_hd <- mask(vatikasrtm_mask ,mask= adm, inverse=F); #plot(vatikasrtm_mask)

# On récupère les zones potentiellement ombragées

vatikaaspect = terrain(vatikasrtm_hd, opt='aspect')# Récupérer les zones ombragées

# On récupère les pentes

vatikaslope = terrain(vatikasrtm_hd, opt='slope')# récupérer les pentes

# On combine le tout et on implique un niveau d'ensoleillement

vatikahill_temp = hillShade(vatikaslope , vatikaaspect , 10, 310) # Soleil couchant (position à 290, hauteur faible donc < 45

# Affichage

plot(adm,xlim=c(23.05,23.20),ylim=c(36.422,36.53))

plot(vatikahill, col=grey(0:100/100), legend=T,add=T);plot(adm,add=T)

Par défaut, le niveau de gris va de 0 à 1 selon une relation linéaire proportionnelle à la pente et son abri par rapport à la lumière.

Toutefois, on peut changer la courbe de niveau de gris pour s'assurer que les pentes les intenses restent bien sombres tandis que le reste soit bien clair.

Il suffit pour cela de jouer sur la fonction fun_grey_col (fonction maison !).

Fonction pour induire un niveau de gris à courbure exponentielle pour colorier les flans de pente.

fun_grey_col <- function(center=0.5,curve=1,plot=F) {

m <- 100 / curve

distribution <- ((exp(0:100/m)-1)/(exp(curve)-1))

ajustement <- center - mean(distribution)

distribution <- distribution + ajustement

temp <- distribution[distribution>center]-center

temp <- temp/max(temp)*(1-center)+center

distribution[distribution>center] <- temp

temp <- distribution[distribution<center]-min(distribution)

temp <- temp/max(temp)*center

distribution[distribution<center] <- temp

grey_col <- grey(distribution)

if (plot==T) {plot(distribution,pch=16,cex=3,col=grey_col)}

# {plot(c(1:length(grey_col)))}

return(grey_col)

}

On peut définir le niveau qui sera le niveau moyen de niveau de gris (ici 0.65 au lieu de 0.5 par défaut) et une fonction de courbure. Plus curve s'approche de 0 (et doit être toujours >0), plus la relation est linéaire. Plus curve est élevée, plus la courbure est exponentielle.

fun_grey_col(0.65,1,plot=T)

plot(adm,xlim=c(23.05,23.20),ylim=c(36.422,36.53))

plot(vatikahill, col=fun_grey_col(0.65,1,plot=F), legend=T,add=T)

plot(adm,add=T)

On voit ainsi que la carte est plus claire et on peut ainsi adapter le niveau de gris à ses besoins.

Si maintenant on veut superposer le grisage des pentes à une carte déjà colorée, il faudrait introduire un niveau de transparence.

plot(vatikasrtm, legend=F, axes=F, xlab="", ylab="",xlim=c(23.05,23.20),ylim=c(36.422,36.53),box=F)

# box = F ; pour ne pas afficher le cadre

plot(vatikahill, col=mon_gris, legend=F,add=T,alpha=0.5)

# alpha pour ajouter de la transparence

On peut aussi additionner un niveau d'alpha directement aux couleurs produits par la fonction grey ou fun_grey_col().

On peut faire cela avec la fonction alpha.col de la librairie plot3D.

require("plot3D")

mon_gris <- alpha.col (col = fun_grey_col(0.65,1,plot=F), alpha = 0.5)

3- Faire une carte de densité

Il ne faut pas confondre un raster (où chaque point correspond à une valeur z à coloriser sur une carte à 2 dimensions (x et y) avec une carte de densité où il va falloir construire un raster qui permettra de coloriser en fonction du nombre de points qui apparaissent dans une même zone.

Mieux vaut partir d'un exemple... (pas sous ggplot2, car ce serait une vraie boîte noire !)

Sous-titre clefs : Comptabiliser les points sur une aire, les points qui se superposent à un polygone .

3.1 - Exemple 1 - en choroplèthe : compter les points pour chaque sous polygone

0) Voici une carte avec des points que l'on va compter région par région

library(raster)

adm <- getData("GADM",country="FRA",level=1)

plot(adm)

# Niveau 1 - bis - Carte de France avec points

bbox(adm)-> xy1

x <- rnorm(100,mean(xy1[1,]),2)

y <- rnorm(100,mean(xy1[2,]))

plot(adm)

points(x,y,col="red",pch=16,cex=2)

1.a) Méthode 1 - Compter les points pour chaque polygone avec sf

library(raster)

adm <- getData('GADM', country='FRA', level=1)

library(sf)

st_as_sf(adm) -> sff # Conversion au format sf

pts<- st_as_sf(data.frame(x,y), coords = c("x", "y"),crs=st_crs(sff))

inter <- st_intersects(sff, pts) # Points par région

niveau1 <- sapply(X = inter, FUN = length) # Comptage

1.b) Méthode 2 - Compter les points pour chaque polygone avec sp

library(raster)

adm <- getData('GADM', country='FRA', level=1)

install.packages("GISTools")

library("sp") ; library("GISTools")

p<-SpatialPoints(data.frame(x,y))

poly.counts(p, adm) -> res

niveau1 <- setNames(res, adm$NAME_1)

2.1) Affichage avec coloration manuelle

# Coloration

cut(niveau1,breaks=6)->niveau

colfunc <- colorRampPalette(c("white","red"))

colors<- colfunc(6) %>% .[as.numeric(niveau)]

plot(adm,col=colors)

Autre exemple pour aller plus loin...

Cliquer pour voir comment réaliser cet exemple

2.2) Affichage - Autre méthode avec choroLayer de la librairie cartography

library(cartography)

sff <- cbind(sff,niveau1) # reprendre niveau plus haut

plot(st_geometry(sff))

choroLayer(sff,var="niveau1")

3.2 - Exemple 2 - en densité (type "heatmap") : avec lissage KernelDensity

Voici une carte avec plein de points, je voudrais en faire une carte de densité car il y a certainement des zones où les points se superposent !

library(raster)

# Contours administratifs

adm_gr <- getData('GADM', country='GRC', level=1)

# La ligne qui suit permet de limiter la zone traitée à une région,

#ici le Péloponnèse (La carte sera plus rapideà tracer)

adm <- adm_gr[match(toupper(adm_gr$NAME_1[7]),toupper(adm_gr$NAME_1)),]

plot(adm,xlim=c(23.05,23.20),ylim=c(36.422,36.53))

# Simulons les points ! (rnorm sert à simuler des valeurs : pas d'inquiétude !)

x1 <- rnorm(400,mean(c(23.05,23.20)),0.026) ; x2 <- rnorm(200,23.17,0.01) ; y1 <- rnorm(400,mean(c(36.422,36.53)),0.019) ; y2 <- rnorm(200,36.46,0.01) ; x <- c(x1,x2) ; y <- c(y1,y2)

points(x,y,col="red",cex=0.5,pch=16)

Je vais maintenant calculer à partir de ces points des zones de densité

library(MASS) # va générer le kernel (penser à l'installer) (kernel density 2D)

my_kernel <- kde2d(x, y,n=200)

# n permet de lisser l'image. n grand, image lisse - n petit, image moche

A ce stade, non seulement, c'est moche, mais il va falloir convertir le kernel pour le superposer à la carte et le recoloriser.

Convertir le kernel et l'afficher sur la carte

my_raster <- raster(my_kernel) # conversion en raster

plot(adm,xlim=c(23.05,23.20),ylim=c(36.422,36.53))

plot(my_raster,add=T,alpha=0.5,legend=F)

Bon, la recolorisation s'est faite par défaut. Alpha permet de rendre le raster transparent pour voir la carte derrière. Légende permet d'éliminer la légende, ou pas.

Recoloriser le raster et le découper aux contours de la carte

Comme ci-dessus sur la page, je créé mon gradient de couleur

liste_colors <- c("black","blue","yellow","red")

colfunc<-colorRampPalette(liste_colors,interpolate="spline")

colors <- colfunc(1000)

Comme ci-dessus, j'utilise mask pour la découpe

my_raster_decoupe <- mask(my_raster,mask= adm, inverse=F)

Et j'affiche ! Et voilà le travail !

plot(my_raster_decoupe,legend=F,axes=F,box=F,col=colors)

plot(adm,add=T,lwd=4,border="#86091D")

3.3 - Exemple 3 - Lorsque la couleur correspond à un lissage de la valeurs de points définis par une variable

1) Voici une carte où chaque point à sa valeur propre (variable dim simulée ici)

adm <- getData("GADM",country="FRA",level=1)

# Simuler des points sur la carte

bbox(adm)-> xy1

x <- rnorm(100,mean(xy1[1,]),2)

y <- rnorm(100,mean(xy1[2,]))

dim <- round(x+y,0) # Dim serait une valeur associée à chaque point (ici on additionne long et lat, ce qui n'a pas de sens...

# Découper les

cut(dim,breaks=10)->niveau

colfunc <- colorRampPalette(c("white","red","black"))

colors<- colfunc(10) %>% .[as.numeric(niveau)]

plot(adm)

points(x,y,col=colors,pch=16,cex=2)

legend(x="topleft",legend = levels(niveau), cex=0.8,fill=colfunc(10),bty="n")

2) Pour chaque région (polygone), faire la moyennes des valeurs des différents points qui se superposent à ce polygone

Option 1 : Réflexion sp

pts<-st_as_sf(data.frame(x,y,dim), coords = c("x", "y"),crs=st_crs(sff))

as(pts,Class="Spatial")->pts_sp

require(rgeos)

cptg_med = aggregate(pts_sp["dim"], adm, median)

cut(cptg_med$dim,breaks=10)->niveau

Option 2 : Réflexion sf

lmean <- function(vector) {

return(mean(dim[vector]))

}

niveau <- sapply(X = inter, FUN = lmean) # Comptage

3) Coloration par région

# coloration

cut(cptg_med$dim,breaks=10)->niveau

colfunc <- colorRampPalette(c("white","red","black"))

colors<- colfunc(10) %>% .[as.numeric(niveau)]

plot(adm,col=colors)

4) Même affichage avec interpolation (on ne tient plus compte des régions)

xy<- coordinates(pts_sp)

#xy <- data.frame(x,y)

colnames(xy)

resolution <-0.1

library(akima)

a <- interp(x=xy[,1], y=xy[,2], z=dim,

xo=seq(min(xy[,1]),max(xy[,1]),by=resolution),

yo=seq(min(xy[,2]),max(xy[,2]),by=resolution), duplicate="mean")

image(a,col = colfunc(50)) # couleurs modifiables

filled.contour(a, color.palette=colfunc,axes=T)

a <- raster(a)

plot(adm)

plot(a,add=T,alpha=0.5,col=colfunc(50))

Si je veux aller plus loin en réalisant des coupes ou des cartes 3D, je passe au niveau 5 :

Google Sites

Report abuse