Analyse spatiale avec R

Contents
Prise en main de l’analyse spatiale avec R 1
Objectifs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
Analyse en mode vecteur . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
Chargement des packages et ouvertue des jeux de données . . . . . . . . . . . . . . . . . . . . 2
Sélection attributaire . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
Jointure spatiale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
Opérateurs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
Application sur les données de l’exercice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
Visualisation et export du résultat . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
Analyse en mode raster . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
Chargement des packages et lecture des fichiers . . . . . . . . . . . . . . . . . . . . . . . . . . 4
Découper et reclasser le raster . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
Conversion du vecteur au format raster . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
Rééchantillonnage et décision finale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

Traitement des images satellitaires 6

Extraction des valeurs et profil spectral . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
Seuillage et reclassification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
Classification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

Modélisation spatiale 8
Chaines de Markov pour la projection des état futures . . . . . . . . . . . . . . . . . . . . . . . . . 8
Géostatistique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
Random Forest . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
Optimisation des traitements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
Considérations relatives aux capacités de calcul . . . . . . . . . . . . . . . . . . . . . . . . . . 10

11
<div style=“page-break-after: always; visibility: hidden”>
\pagebreak
</div>

Prise en main de l’analyse spatiale avec R

Objectifs

A travers quelques cas, ce recueil vise à initier aux analyses spatiales et au traitement des données spatiales
avec le langage R.

1
Analyse en mode vecteur

Soit l’exercice relative à l’installation de magasin de glace.

Chargement des packages et ouvertue des jeux de données

La lecture des fichiers vectoriels se base sur l’usage du package sf

# chargement des packages
library(dplyr)
library(tidyr)
library(sf)
# Répertoire du travail
fname <- "/Users/said/Desktop/formations/gis_iav/data_osm"
setwd(fname)

# Lectrue du fichier shapefile

cr <- sf::st_read('communes_maroc2.shp')
# découverte du jeu de données
cr[2, "Nom_Commun"]
cr[4, ]

Sélection attributaire

La syntaxe conventionnelle de l’indexation est utilisée pour l’extraction des éléments d’un objet spatial de
type sf.
# Selection attributaire :
ajdir <- cr[cr$Nom_Commun == "Ajdir",]
ajdir

Jointure spatiale

Si on souhaite extraire les communes qui jouxtent la commune d’Ajdir ou éventuellement qui touchent ladite
commune. L’argument op peut référer à l’un des prédicats spatiaux conventionnels.
# selection spatiale
commune_autour_ajdir<-cr[ajdir, op = st_touches] #

Les opérateurs spatiaux supposent que les géométries sont valides. La ligne de code précédente à cause de
deux problèmes topologiques pose des problèmes.
# Choix de seuls les communes dont la géométrie est valide
cr2<-cr[st_is_valid(cr),]
# Selection spatiale
commune_autour_ajdir<-cr2[ajdir, op = st_touches]
# usage des pipes
commune_autour_ajdir2<-cr2 %>% filter(lengths(st_touches(., ajdir)) > 0)

# il est aussi possible de définir la jointure avec une jointure spatiale explicite
p_aj<-st_join(cr2, ajdir, join=st_touches) # si argument join est omis, st_intersects sera utilisé
plot(st_geometry(ajdir))
plot(st_geometry(commune_autour_ajdir), col='red', add=T)

2
Opérateurs

Différents opérateurs spatiaux peuvent être appliqués au objets vectoriels: projection, croisement, tampon . . .

Application sur les données de l’exercice

# definition du chemin de travail

pdir2<- "/Users/said/Desktop/formations/gis_iav/qgs@iav_lectures/tpr/data/Jane_exercice"
setwd(pdir2)
########################
# Counties
## lecture
counties<-st_read('counties.shp')
plot(st_geometry(counties))
# filtre avec tidyr : critere sur les counties
count_ok<-counties%>% filter(AGE_18_64>25000)
# affichage des counties répondant au critère
plot(st_geometry(count_ok), col='green', add=T)
### zones recreatives
##############@
# critere sur les villes
cities<-st_read('cities.shp')
cit_ok <-cities%>% filter(CRIME_INDE<=0.02, UNIVERSITY==1)
cit_okc<-cit_ok[count_ok, op=st_within]
##############@
## Critère sur es zones récréatives

rec<-st_read('recareas.shp')
st_crs(rec)<-4267
rec<-st_transform(rec,4455) ## repojection
dist<-52800 # 10 miles
rec_buf<-st_buffer(rec, dist)
st_crs(rec_buf)<-4455
###################
# routes interstates
inter<-st_read('interstates.shp')
st_crs(inter)<-4267
inter<-st_transform(inter,4455)
dist<-52800*2 # 20 miles
inter_buf<-st_buffer(inter, dist)
st_crs(inter_buf)<-4455

#cit_ok2<-st_transform(cit_ok, 4455)
cit_ok2<-st_transform(cit_okc, 4455)

# jointure spatiale finale approche 1

cit_ok2 <- cit_ok2%>%st_join(., inter_buf, left=F) %>%
st_join(., rec_buf, left=F)
#length(cit_ok2[[1]])
choix_final<-cit_ok2%>%distinct()
###
#approche 2
cit_ok3<-cit_ok2[inter_buf, op=st_within]

3
cit_ok3<-cit_ok3[rec_buf, op=st_within]
#length(cit_ok2[[3]])
choix_final<-cit_ok3%>%distinct()

Visualisation et export du résultat

plot(choix_final) # pourquoi autant de graphes?

plot(st_geometry(choix_final)) # Quelle est la différence?
dev.off()
########
#plot(st_geometry(counties))
plot(st_geometry(st_transform(count_ok, 4455)))
plot(st_geometry(inter_buf), col='yellow', add=T )
plot(st_geometry(rec_buf), col='green', add=T)
plot(st_geometry(st_transform(counties, 4455)))
plot(st_geometry(choix_final), add=T, col='red')

\newpage

Analyse en mode raster

Chargement des packages et lecture des fichiers

library(sf)
library(raster)

# definir le repertoire de travail

pdir2<- "/Users/said/Desktop/anfcc/iav2020/labs/lab7/raster_tp"
setwd(pdir2)
# Lecture des rasters
r1<-raster("aster1.tif")
r2<-raster("aster2.tif")
r3<-raster("dem_rsk.tif")
# mosaicage
rm<-mosaic(r1, r2, r3, fun=mean)
# ecriture du résultat
writeRaster(rm, 'analyse/r_mosaic.tif', overwrite=T)

Découper et reclasser le raster

# lecture du fichier shapefile relatif à la zone d'étude : ZE

pdir2<- "/Users/said/Desktop/anfcc/iav2020/labs/lab7/raster_tp"
setwd(pdir2)

ze0<-st_read('test_ze.shp')
# conversion du type sf au type sp compatible avec le package raster
ze<-as(ze0, 'Spatial')
#rm<-raster('analyse/r_mosaic.tif')
# altitude

4
alt_ok<-crop(rm, ze)
alt_ok[alt_ok<149]<-0
alt_ok[alt_ok>=149]<-1
writeRaster(alt_ok, 'analyse/r_alt_ok.tif', overwrite=T)
#Pente < 15 degre
### ?probleme
#sl2<-terrain(rm)

# Calcul de la pente
sl<-terrain(rm, opt=c('slope'), unit='degrees')

sl_ok<-crop(sl, ze)
sl_ok[sl_ok<15]<-1
sl_ok[sl_ok>=15]<-0
writeRaster(sl_ok, 'analyse/r_sl_ok.tif', overwrite=T)
####
# traitement des fichiers relatifs à la vitesse du vent
#### Vitesse de vent > 5ms-1
r1<-raster("wind_gharb.tif")
r2<-raster("wind_rszz.tif")
rm<-mosaic(r1, r2, fun=mean)
wind_ok<-crop(rm, ze)
wind_ok[wind_ok<5]<-0
wind_ok[wind_ok>=5]<-1
writeRaster(wind_ok, 'analyse/r_wind_ok.tif', overwrite=T)

Conversion du vecteur au format raster

#####
# Critère des routes
pdir2<- "/Users/said/Desktop/anfcc/iav2020/labs/lab7/raster_tp"
setwd(pdir2)
routes<-st_read('routes_primaire_second26191.shp')
df2<-st_buffer(routes, 150)
df2<-st_transform(df2, 4326)
dist2route<-wind_ok
dist2route[dist2route>=0]<-0
dist2route<-rasterize(as(df2, 'Spatial'), dist2route, field=0, background=1)
dist2route<-crop(dist2route, ze)
writeRaster(dist2route, 'analyse/r_dist2route.tif', overwrite=T)
####
# contraintes du type d'utilisation du sol
####
library(dplyr)
lu<-st_read('osm_l_use.shp')
lu<- lu%>%filter(type %in% c('cemetery', 'commercial', 'construction', 'industrial', 'residential'))
#lu<-st_join(lu, ze0)
lu<-lu[ze0, op=st_within]
mamo<-st_read('maamo_26191.shp')
maamo<-st_transform(mamo, 4326)
maamo$osm_id<-rep(NA, 3)
maamo$name<-rep('maamora', 3)

5
maamo$type<-rep('foret', 3)
maamo2<-maamo%>%select(osm_id, name, type, geometry)
#luf<-st_union(st_geometry(maamo),st_geometry( lu))
luf<-rbind(lu, maamo2)
luf2<-wind_ok
luf2[luf2>=0]<-0
luf2<-rasterize(as(luf, 'Spatial'), luf2, field=0, background=1)
writeRaster(luf2, 'analyse/r_luf.tif', overwrite=T)

Rééchantillonnage et décision finale

#######
## Rééchantillonnage du jeu final et extraction de carte d'aptitude
######
sl_ok1 = resample(sl_ok, wind_ok, method = "bilinear")
sl_ok1 = resample(sl_ok, wind_ok, method = "bilinear")
## resample : ngb vs bilinear
alt_ok1 = resample(alt_ok, wind_ok, method = "bilinear")
aptitude<-dist2route*wind_ok*sl_ok1*alt_ok1*luf2
# verifier les valeurs : probleme ??
writeRaster(aptitude, 'analyse/r_suita.tif')
plot(aptitude)

\newpage

Traitement des images satellitaires

##
library(raster)
fname <- "/Users/said/Desktop/formations/gis_iav/GIS-DATA/sat"
setwd(fname)
download.file(url = 'https://raw.githubusercontent.com/GeoScripting-WUR/VectorRaster/gh-pages/data/lands
unzip('data/landsat8.zip', exdir = "data")

landsatPath <- list.files(path = "data/", pattern = glob2rx('LC8*.grd'), full.names = TRUE)

wagLandsat <- brick(landsatPath)
# plotRGB ne supporte pas les valeurs <0
wagLandsat[wagLandsat < 0] <- NA
names(wagLandsat) <- c("band1","band2","band3","band4","band5","band6","band7")

plotRGB(wagLandsat, 5, 4, 3)
nf <- layout(matrix(c(1,0,2), 1, 3, byrow = TRUE), width = c(1,0.2,1), respect = TRUE)
plotRGB(wagLandsat, r = 3, g = 2, b = 1, axes = TRUE, stretch = "lin", main = "CC1")
plotRGB(wagLandsat, r = 4, g = 3, b = 2, axes = TRUE, stretch = "lin", main = "FC")
#### calculs
cellStats(wagLandsat[[2]], stat=max)
cellStats(wagLandsat[[2]], stat=mean)
maxValue(wagLandsat[[2]])
summary(wagLandsat)
### brick avec un nombre réduit de bandes

6
gewata <- brick(wagLandsat[[2]], wagLandsat[[3]], wagLandsat[[4]])

hist(gewata, maxpixels=1000)
par(mfrow = c(1, 1))

pairs(gewata)
par(mfrow = c(1, 1))
# Calcul du ndvi
ndvi <- overlay(wagLandsat[[4]], wagLandsat[[3]], fun=function(x,y){(x-y)/(x+y)})
plot(ndvi)

Extraction des valeurs et profil spectral

# Points d'apprentissage : type de culture/utilisation des terres ....

samp <- readRDS('samples.rds')

df <- extract(wagLandsat,samp)

# affichage de la valeur de reflectance

head(df)

df <- round(df)
# creation df vide dans lequel stocker les valeurs moyennes des RS
ms <- matrix(NA, nrow = length(unique(samp$id)), ncol = nlayers(ss))

for (i in unique(samp$id)){
x <- df[samp$id==i,]
ms[i,] <- colMeans(x)
}
# noms des col/lig
rownames(ms) <- unique(samp$class)
colnames(ms) <- names(ss)

# vecteur des couleurs

mycolor <- c('cyan', 'darkgreen', 'yellow', 'burlywood', 'darkred', 'darkgray', 'blue', 'lightgreen')

# Plot
plot(1, ylim=c(300, 4000), xlim = c(1,10), xlab = "Bandes", ylab = "Reflectance", xaxt='n')

# X
axis(1, at=1:10, lab=colnames(ms))
# Ajouter les valeurs spectrales
for (i in 1:nrow(ms)){
lines(ms[i,], type = "o", lwd = 3, lty = 1, col = mycolor[i])
}
title(main="Profil Spectral ...", font.main = 2)# Titre
legend("topleft", rownames(ms), cex=0.8, col=mycolor, lty = 1, lwd =3, bty = "n")# legende

7
Seuillage et reclassification

veg <- calc(ndvi, function(x){x[x < 0.4] <- NA; return(x)})

plot(veg, main = 'Couver végétal')
vegc <- reclassify(veg, c(-Inf,0.2,0, 0.2,0.3,1, 0.3,0.4,2, 0.4,0.5,3, 0.5, Inf, 4))
plot(vegc,col = rev(terrain.colors(4)), main = 'Classe de couvert (seuils NDVI) ')

Classification

nr <- getValues(ndvi)
nr.km <- kmeans(na.omit(nr), centers = 10, iter.max = 500, nstart = 3, algorithm="Lloyd")
knr <- ndvi
knr[] <- nr.km$cluster
plot(knr, main = 'Classification non supervisée')

###############
# avec le package RStoolbox
library(RStoolbox)
cns<-unsuperClass(ndvi, nSamples=200, nClasses=6, nStarts=5)
ggR(cns$map, geom_raster=T, forceCat=T) # affichage

Classification non supervisée

Classification supervisée Exemples de classification supervisée

https://gfc.ucdavis.edu/events/arusha2016/_static/labs/day4/day4_lab1_remote-sensing.pdf
http://amsantac.co/blog/en/2015/11/28/classification-r.html
https://geoscripting-wur.github.io/AdvancedRasterAnalysis/#Advanced_Raster_Analysis
https://rpubs.com/daniballari/imagenes_eh
##
#library(rgdal)
#terrain<-readOGR (".", "terrain")
#clas_sup<-superClass(image, trainData=terrain, responseCol="classe", model=mlc)

\newpage

Modélisation spatiale

Chaines de Markov pour la projection des état futures

Supposons que vous disposez des occupations des sols à différentes pas temporels et vous voulez calculer la
probabilité de transition d’une utilisation des terres à une autre. Ces transitions vous permettront de projeter
l’état futur de vos occupations des sols.

8
# carte occupation des terres au niveau du répertoire actif

nlcd_2001 = raster('raster/nlcd_2001.tif')
nlcd_2006 = raster('raster/nlcd_2006.tif')
nlcd_2011 = raster('raster/nlcd_2011.tif')

# Lecture des limites de la ZE

# arguments readOGR pour ESRI Shapefile:
# data source name (dsn) = repertoire
# layer = shapefile filename without shp extension
cr = readOGR(dsn='vecteur', layer='communes')
# extraction des valeurs des pixel (ligne) selon les années (colonne)
# choisir une commune, extraire son polygone et son emprise :bbox
com = 'Ajdir'
poly = cr[cr@data$Nom_Commun == com,]
# poly<-as(ajdir, "Spatial") ### Ajdir : de l'exercie precedent
bbox = bbox(poly)

# mask par la commune

lc01 = mask(crop(nlcd_2001, poly), poly)
lc06 = mask(crop(nlcd_2006, poly), poly)
lc11 = mask(crop(nlcd_2011, poly), poly)

v = data.frame(
lc01 = values(lc01),
lc06 = values(lc06),
lc11 = values(lc11))

# matrice from 2001 to 2006

m = table(v[,c('lc01','lc06')])

kable(m, format='pandoc', caption='pixels ayant changé de type de couvert entre 2001 (lignes) et 2006 (c

# Matrice transition
P = as.matrix(m / rowSums(m))
write.csv(P, 'P_lc01to06.csv')

# Table (kable) de transition à utiliser pour les prédictions

kable(P, format='pandoc', caption='Matrice de probabilité de transition (_P_) extraite des CC entre 2001

Calcul des matrices de transition

# prediction de l'état en 2011 sur la base de la matrice de transition 2001-> 2006

projections = table(v$lc06) %*% P

Prédictions

Géostatistique

Les packages conseillés

9
 • spatial
• kriging
• automap
# voir exemple automap sur slides

Random Forest

\newpage
library(randomForest)
# df : dataframe; 1:5 colonnes variables explicatives, classe : variable de réponse
modelRF <- randomForest(x=df[ ,c(1:5)], y=df$classe,
importance = TRUE)

Optimisation des traitements

Considérations relatives aux capacités de calcul

La contrainte des ressources peut imposer soit une parallélisation des calculs et éventuellement un traitement
des rasters sur le disque et non plus en mémoire. L’exemple suivent illustre ces deux cas de figure.
# Paralléisation et calcul sur le disque d
library(gdalUtils)
library(parallel)
library(raster)
###
# Accès aux fichiers
###
mnt<-raster('/Users/said/Desktop/als/drone_azrou/Azrou/MNT/Azrou02_dtm.tif')
mns<-raster('/Users/said/Desktop/als/drone_azrou/Azrou/MNS/Azrou02_dsm.tif')

#Définition des paramètres

beginCluster(n = 10) # nombre de coeur du processeur à utiliser
rasterOptions(todisk = TRUE) # traitement des raster sur le disk
# sans parallélisation
mnt.new = resample(mnt, mns, "bilinear")
# avec parallélisation
mns.new = clusterR (mns, resample, args=list( mnt, "bilinear"))
writeRaster(clusterR (mns, resample, args=list( mnt, "bilinear")),
'/Users/said/Desktop/als/drone_azrou/Azrou/MNT/Azrou02_mns_resamp.tif',
datatype='INT2U')

# parallélisation avec une fonction personnalisée

cat('calcul')
chmf <- function(x, y) (x-y)
chm <- clusterR(mnt, calc, args=list(fun=chmf, y=mns.new))
writeRaster(chm, '/Users/said/Desktop/als/drone_azrou/Azrou/MNT/Azrou02_chm.tif', datatype='INT2U')
rasterOptions(todisk = F)
endCluster()

10
11