Journées Deep Learning, Télédétection et Paysage (PAYOTE/MODELIA) - 28 juin 2019
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 à:
# 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)
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # 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.
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)
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)
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")
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)
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.
cloud1_proj=aggregate(z~x+y,data=cloud1,FUN=max) cloud1_proj=merge(cloud1_proj,cloud1,all.y=F)
# 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)
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)
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)
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}
# 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)
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)
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")
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")
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