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

• Afficher une région de la carte des altitudes

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)

• 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)

• 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)

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

# 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)

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

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.

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 .

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) 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