Journées Deep Learning, Télédétection et Paysage (PAYOTE/MODELIA) - 28 juin 2019

Objectifs de l'atelier

Ecohydrologie des fossés agricoles

L'atelier présenté ici à vocation à partager notre expérience de la manipulation de nuages de points en 3D (issues de Lidar, SfM-MVS) pour analyser les caractéristiques de sections de fossés agricoles.

L'atelier se base sur un extrait de jeu de données original, acquis en 2015 lors d'une campagne de suivi par photogrammétrie d'un linéaire de 200 mètres de fossé.

Dans cet atelier, à partir de nuages de points 3D, nous apprendrons à:

  • charger des données 3D
  • visualisation 3D
  • réaliser des projections selon différents plans
  • visualisation 2D
  • définir des sections et les analyser

Packages et données

Packages R

# Packages d'analyse et manipulation de données spatiales
library(sp)
library(raster)
library(maptools)
library(rgdal)
library(rgeos)

# Visualisation de données 3D
library(rgl)

# Manipulation avancée de tableaux de données
library(data.table)

# Algorithme de machine learning pour la classification
library(randomForest)

# Transparence des couleurs
library(scales)

# Cartographie interactive
library(tmap)

Chargement des données spatiales

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Système de projection et offset                                                 |----
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
shiftXYZ=c(726000,6266100,80)
proj=CRS("+init=epsg:2154") # lambert 93/RGF93
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Fichiers spatiaux                                                               |----
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
vZON=readOGR("IN/SHP/","areas")
vTAL=readOGR("IN/SHP/","talweg")
vFOOT=readOGR("IN/SHP/","footprints")
projection(vZON)=projection(vTAL)=projection(vFOOT)=proj

Les données spatiales, sous forme de fichiers shapefiles, renseignent sur la zone d'étude et le découpage en cubes du nuage de point complet.

Visualisation de la zone d'étude

mapBasic=
  tm_scale_bar()+ 
  tm_view(set.zoom.limits=c(5,100))+
  tmap_options(basemaps = c(Imagery = "Esri.WorldImagery",
                            Canvas = "OpenStreetMap"))
tmap_mode("view")
mapBasic+tm_shape(vZON,name="fossé")+tm_fill(col="Situation")+
         tm_shape(vTAL,name="talweg")+tm_lines(col="black",lwd=2)+
         tm_shape(vFOOT,name="empreintes des nuages")+tm_fill(col="red",alpha=0.1)+
         tm_borders(col="red",alpha=0.5)+
         tm_view(set.view=17)

Selection des nuages de points

Considérant que les nuages de points sont stockés dans un dossier à part, l'idée est de ne sélectionner que certains nuages couvrant un polygone donné, ici vAREA:

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Selection des nuages                                                            |----
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Définition de l'aire de selection des nuages
vAREA = SpatialPolygons(
  list(Polygons(
    list(Polygon(
      matrix(
        c(726064,726064,726070,726070,
          6266248.7,6266250.5,6266250.5,6266248.7),ncol=2))),"S1")))
projection(vAREA)=proj
# Selection des noms de fichiers correspondant
vFOOT_sel=vFOOT[c(gIntersects(vAREA,vFOOT,byid=T)),"ID"]
namesCloud=as.character(vFOOT_sel$ID)

Chargement et concaténation des nuages de points

Utilisation du package utils

cloud1=rbindlist(lapply(namesCloud[grep("2014-04-25",namesCloud)],function(j)
                        read.table(paste("IN/CLOUDS/",j,sep=""))))
names(cloud1)=c("x","y","z","R","G","B")
cloud2=rbindlist(lapply(namesCloud[grep("2014-06-11",namesCloud)],function(j)
                        read.table(paste("IN/CLOUDS/",j,sep=""))))
names(cloud2)=c("x","y","z","R","G","B")

Utilisation du package data.table

cloud1=rbindlist(lapply(namesCloud[grep("2014-04-25",namesCloud)],function(j)
                        fread(paste("IN/CLOUDS/",j,sep=""))))
names(cloud1)=c("x","y","z","R","G","B")
cloud2=rbindlist(lapply(namesCloud[grep("2014-06-11",namesCloud)],function(j)
                        fread(paste("IN/CLOUDS/",j,sep=""))))
names(cloud2)=c("x","y","z","R","G","B")

tmap_mode("view")
mapBasic+tm_shape(vZON,name="fossé")+tm_fill(col="Situation")+
         tm_shape(vFOOT,name="empreintes des nuages")+tm_borders(col="red",alpha=0.1)+
         tm_shape(vAREA,name="aire selectionnée")+tm_fill(col="red",alpha=0.5)+
         tm_view(set.view=17)

Visualisation du nuage de points en 3D

Le nuage de points concaténé et correspondant à l'emprise vAREA peut ensuite être visualisé grâce à la fonction points3d:

points3d(cloud1,col=rgb(cloud1[,c("R","G","B")],maxColorValue=255)) # in color
points3d(cloud2,col=alpha("red",0.5))

On peut ainsi observer la superposition des deux nuages de points et les éventuelles différences entre nuages pouvant être expliquée par des facteurs biotiques: terriers d'animaux sauvages et végétation. Ainsi l'étude du relief fin permet de comprendre certaines relations biotiques-abiotiques.

Plan horizontal

Projection des nuages de points

Utilisation du package stats

cloud1_proj=aggregate(z~x+y,data=cloud1,FUN=max)
cloud1_proj=merge(cloud1_proj,cloud1,all.y=F)

Utilisation du package data.table

# Utilisation du package data.table
cloud2_proj=cloud2[,j=list(z=max(z)),by=list(x,y),mult="first"]
cloud2_proj=cloud2[cloud2_proj,mult="first",on=.(x,y,z)]
cloud2_proj[,c("x","y","z"):=list(x+shiftXYZ[1],y+shiftXYZ[2],z+shiftXYZ[3])]
# Transformation des points projetes en un raster
coordinates(cloud2_proj) <- ~x+y
gridded(cloud2_proj) <- TRUE
r2=stack(cloud2_proj);  names(r2)=c("height","R","G","B")
projection(r2)=proj
# Aggregation d'un facteur 10 des données (resolution mm -> cm)
r2_ag=aggregate(r2,fact=10)
# Calcul de la carte d'ombrage pour la visualisation
r2_hill=hillShade(terrain(r2_ag,"slope"),terrain(r2_ag,"aspect"),angle=0)

Classification du raster

Une fois la projection réalisée, on cherche à distinguer la végétation du sol en utilisant un algorithme de machine learning de type randomForest:

rORTHO=r2_ag[[c("R","G","B")]]
plotRGB(rORTHO)

trainGROUND=drawPoly()
trainVEG=drawPoly()
dTRAIN=rbind(data.frame(extract(rORTHO,trainGROUND)[[1]],TYPE="ground"),
      data.frame(extract(rORTHO,trainVEG)[[1]],TYPE="vegetation"))
rfMODEL=randomForest(TYPE~R+G+B,dTRAIN)
r2_classes=predict(r2_ag,rfMODEL)

Visualisation

tmap_mode("view")
mapBasic+
  tm_shape(r2_ag[["height"]],name="MNS")+
  tm_raster(palette=topo.colors(20),style="pretty",n=20,legend.show = F)+
  tm_shape(r2_hill,name="Hillshade")+
  tm_raster(palette=colorRampPalette(c("black","white"))(100),
            alpha=0.5,n=100,legend.show = F)+
  tm_shape(r2_classes,name="Classified map")+tm_raster(palette=c("brown","green"),
                                                       alpha=0.5)+
  tm_view(set.view=23)

Plan vertical

Projection

On projete le nuage de point selon le plan vertical (plan xz):

cloud2_proj_vertical=cloud2[,j=list(y=length(y)),by=list(x,z),mult="first"]
coordinates(cloud2_proj_vertical) <- ~x+z
gridded(cloud2_proj_vertical) <- TRUE
r2_vertical=raster(cloud2_proj_vertical)
plot(r2_vertical,legend=F)

Pour obtenir des sections réalistes et perpendiculaires au talweg du fossé, il faut appliquer une fonction de rotation au nuage:

rotateCloud=function(datas=cloudsZ[[1]],theta=optRot$minimum,plan="xy"){
  datasRot=copy(datas)
  rot<-matrix(c(cos(theta), -sin(theta), sin(theta), cos(theta)), nrow=2, byrow=T)
if(plan=="xy")datasRot[,c("x","y")]=data.table(t(rot%*%t(datasRot[,c("x","y"),with=F])))
if(plan=="yz")datasRot[,c("y","z")]=data.table(t(rot%*%t(datasRot[,c("y","z"),with=F])))
if(plan=="xz")datasRot[,c("x","z")]=data.table(t(rot%*%t(datasRot[,c("x","z"),with=F])))
  datasRot}

Création de sections de fossé

# Definition de la resolution
resol=0.2
vTAL=crop(vTAL,vAREA)

# Semis de points le long de la ligne de talweg
vSAMP=spsample(vTAL, n=round(sum(SpatialLinesLengths(vTAL))/resol), type="regular")

# Estimation de l'angle tangentiel de chaque point à la ligne de talweg
vBUF=gBuffer(vSAMP,byid=T,width=resol*3)
vSAMP$Angle=sapply(1:length(vBUF),function(i){
  vPTS=crop(vSAMP,vBUF[i,])
  vPTS=vPTS[gDistance(vSAMP[i,],vPTS,byid=T)>0]
  vPTS=vPTS[sort(order(gDistance(vSAMP[i,],vPTS,byid=T))[1:2])]
  ds=apply(coordinates(vPTS),2,diff)
  atan2(ds[1], ds[2])
})

# Simulation de sections perpendiculaires au talweg
vREF=SpatialLines(list(Lines(list(Line(coords=data.frame(x=c(-1,1),y=c(0,0)))),ID=1)))
vSECTION=lapply(1:length(vSAMP),function(x){
  tmp=elide(vREF,shift=coordinates(vSAMP)[x,],
            rotate=vSAMP$Angle[x]*180/pi,
            center=c(0,0))@lines[[1]]
  tmp@ID=as.character(x)
  tmp})
vSECTION=SpatialLines(vSECTION)
vSECTION$ID=1:length(vSECTION)
projection(vSECTION)=proj

# Creation d'un buffer autour de la section
vSECTION_buf=gBuffer(vSECTION,byid=T,width=resol/3)

Visualisation

tmap_mode("view")
mapBasic+
  tm_shape(r2_ag[["height"]],name="MNS")+
  tm_raster(palette=topo.colors(20),style="pretty",n=20,legend.show = F)+
  tm_shape(r2_hill,name="Hillshade")+
  tm_raster(palette=colorRampPalette(c("black","white"))(100),
            alpha=0.5,n=100,legend.show = F)+
  tm_shape(vSECTION_buf,name="section large")+tm_fill(col="red")+
  tm_view(set.view=23)

Travail sur les sections

Selection et decoupe de la section

ids=4
minTalweg=min(cloud1$z)
xyS=c(coordinates(vSAMP)[ids,]-t(shiftXYZ[c(1,2)]))
cloud1_sel=copy(cloud1)
cloud1_sel[,c("x","y","z"):=list(x-xyS[1],y-xyS[2],z-minTalweg)]
cloud1_sel=rotateCloud(datas=cloud1_sel,theta=-vSAMP$Angle[ids],plan="xy")
cloud1_sel=cloud1_sel[y>=-resol/2 & y<resol/2,]

cloud2_sel=copy(cloud2)
cloud2_sel[,c("x","y","z"):=list(x-xyS[1],y-xyS[2],z-minTalweg)]
cloud2_sel=rotateCloud(datas=cloud2_sel,theta=-vSAMP$Angle[ids],plan="xy")
cloud2_sel=cloud2_sel[y>=-resol/2 & y<resol/2,]
points3d(cloud1_sel,col="red")
points3d(cloud2_sel,col="green")

Calcul des paramètres écohydrauliques

dXZ=cloud1_sel[,j=list(z=max(z)),by=list(x=round(x,2))]
dXZ=dXZ[dXZ$z<=0.7,]
dXZ=dXZ[order(x),]
areaWet1=SpatialPolygons(list(Polygons(list(Polygon(dXZ)),ID="total")))
plot(cloud1_sel[,c("x","z")],asp=1,pch=".")
points(cloud2_sel[,c("x","z")],col=alpha("green",0.2),pch=".")
plot(areaWet1,add=T,col=alpha("blue",0.5),border="transparent")

Si vous êtes intéressés

Vinatier, F., Bailly, J.-S., & Belaud, G. (2017). From 3D grassy vegetation point cloud to hydraulic resistance: Application to close-range estimation of Manning coefficients for intermittent open channels. Ecohydrology, 10(8). https://doi.org/10.1002/eco.1885

Vinatier, F., Dollinger, J., Rudi, G., Feurer, D., Belaud, G., & Bailly, J.-S. (2018). The Use of Photogrammetry to Construct Time Series of Vegetation Permeability to Water and Seed Transport in Agricultural Waterways. Remote Sensing, 10(12), 2050. https://doi.org/10.3390/RS10122050