In ONF-AURA/geoDalle: Estimation des réserves profondes en milieu karstique

knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)

Objet

En milieu karstique, la possibilité d'accéder a des réserves profondes est fondamentales pour la survie de certaines essences comme le sapin et le cèdre. Cette possibilité dépend de la présence d'une dalle compacte non fissurée. Dans ce cas de figure, le pendage de la strate géologique est souvent conforme (parallèle au sol) et la pente du terrain régulière. Partant de ces caractéristiques, l'objet est de prédire le risque de présence d'une telle contrainte.

Données en entrée

• Modèle Numérique de Terrain à résolution de 5m ou moins
• Polygone sf de la parcelle étudiée
• Polygones des strates géologiques meubles (moraines, molasses, éboulis...) issus de la carte géologique vectorisée du BRGM
• Points d'observation des pendages issus de la carte géologique vectorisée du BRGM

library(geoDalle)

dir <- system.file(package = "geoDalle")
mnt <- raster(file.path(dir, "data","mnt0.tif"))
zn_prf <- read_sf(file.path(dir,"data","zn_prf.shp"))
geo_meuble <- read_sf(file.path(dir,"data","geo_meuble.shp"))
p.pendage <- read_sf(file.path(dir,"data","p.pendage.shp"))

Traitements

Rasters prédictifs

Les calculs sont réalisés sur une zone de 500 à 1000 m autour de la parcelle étudie. La détection des zones homogènes est réalisée d'après les dérivés du MNT lissé sur un rayon de 20m. Ces dérivés sont la pente, l'aspect (scindé en composantes Nord/Sud et Est/Ouest) et le plan de courbure, calculé selon la méthode "bolstad" du package spatialEco.

pte <- terrain(mnt,"slope")
st <- stack(c(mnt= resample(mnt,pte), pte = pte))

# rasters lissés à 20m
fw <- focalWeight(st$mnt,20,type = "circle")
st$mnts <- focal(st$mnt, w=fw)
st$ptes <- terrain(st$mnts,"slope")
st$asps <- terrain(st$mnts,"aspect")
st$aspsx <- cos(st$asps)
st$aspsy <- sin(st$asps)
st$cur <- curvature(st$mnts,"bolstad")

Filtres des zones ne comportant a priori pas de dalle

Les zones de dalles ont a priori un rayon de courbure faible: les fortes courbures sont éliminées, ainsi que les zones de roches meubles d'après la carte géologique (molasses, moraines, dépôts, éboulis...).

# suppression des pentes irrégulières
# jugées selon curvature

st$ind <- st$ptes
st$ind[abs(st$cur)>.01] <- NA

# suppression des comblements

if(nrow(geo_meuble)>0){
  st$ind <- mask(st$ind, geo_meuble %>% st_transform(st_crs(st$ind)), inverse = T)
}

qm(st$ind)

Identification des zones de pente homogènes

Les surfaces restantes sont découpées en zones de pentes homogènes.

# découpage en zones homogènes

st$clump <- clump(st$ind)

# scission des unités de + de 1000 px hétérogènes

find_clu_a_reduire <- function(st){
  val_clump <- data.frame(clump=values(st$clump), pte=values(st$ptes),azx=values(st$aspsx), azy=values(st$aspsy)) %>%
    na.omit() %>% group_by(clump) %>%
    summarise(nb = n(),
              m_pte = mean(pte),
              cv_pte = sd(pte)/m_pte,
              sd_aspx=sd(azx) * (m_pte>.2),
              sd_aspy=sd(azy) * (m_pte>.2))

  val_clump %>%
    filter(nb>1000 & (cv_pte>0.1 | sd_aspx>0.05 | sd_aspy>0.05)) %>%
    pull(clump)
}

clu_a_reduire <- find_clu_a_reduire(st)

while (length(clu_a_reduire)>0) {

  big_cl <- st$clump
  big_cl[! values(big_cl) %in% clu_a_reduire] <- NA
  st$clump[big_cl %>% boundaries() == 1] <- NA
  st$clump <- clump(st$clump)

  clu_a_reduire <- find_clu_a_reduire(st)
}

v <- values(st) %>% as.data.frame()
v$x <- xFromCell(st$mnt, 1:ncell(st$mnt))
v$y <- yFromCell(st$mnt, 1:ncell(st$mnt))

# suppression des zones de moins de 50 px

clump20 <- v %>% na.omit() %>% count(clump) %>% filter(n>50) %>% pull(clump)
v$clump[! v$clump %in% clump20] <- NA

# suppression des clumps avec moins de 3 valeurs de x ou y

clump_rm <- v %>% group_by(clump) %>%
  summarise(count_x = length(unique(x)), count_y = length(unique(y))) %>%
  filter(count_x < 3 | count_y < 3) %>% pull(clump)

v$clump[v$clump %in% clump_rm] <- NA

v$clump <- v$clump %>% as.factor() %>% as.numeric()

clu_names <- unique(v$clump) %>% na.omit()
names(clu_names) <- clu_names

Zones homogènes identifiées:

plot_cl <- st$mnt
values(plot_cl) <- v$clump

qm(plot_cl)

identification des plans de courbure dominants de la zone d'étude

Le plan de courbure de chaque zone est estimé.

# estimation du plan dans lequel se situe chaque zone

vplan <- map(clu_names,function(i){

  m <- glm(mnt~poly(x,2)+ poly(y,2), data = v %>% na.omit() %>%
             filter(clump==i))
  p <- predict(m,v)
  p-mean(p)
})

vplan <- do.call(rbind, vplan) %>% as.data.frame()

Les zones sont classées et regroupées selon la ressemblance de leur plan de courbure.

# classification des plans

d <- dist(vplan %>% select(sample.int(ncol(vplan),1000)))
hd <- hclust(d, method = "ward.D")
brk_k <- classBreaks(hd$height,2,"std")[2]

k <- sum(hd$height>=brk_k)
cl_hd <- cutree(hd,k)

plot(hd)
rect.hclust(hd,k=k,border="blue")

v$m_cl <- NULL
v <- v %>% left_join(data.frame(clump = as.numeric(as.character(names(cl_hd))), m_cl = cl_hd))

plan_cl <- st$mnt
values(plan_cl) <- v$m_cl

sum_clu <- v %>% na.omit() %>%  group_by(m_cl) %>% summarise(
  nb = round(n()),
  pte_m = mean(pte),pte_sd=sd(pte),
  cur_m = mean(cur)*1000, cur_sd=sd(cur)*1000,
  mnt_m= mean(mnt), mnt_sd=sd(mnt),
  aspx = mean(aspsx), aspy = mean(aspsy)

) %>% arrange(desc(nb))

knitr::kable(sum_clu, digits = 2)

comparaison avec les points d'observation de pendages

Dans le cas de contextes pentus, les points d'observation situés sur un versant opposé à celui de la parcelle étudiée sont éliminés.

# comparaison avec données pendage BRGM -----------------------------------------

# points situés sur le même versant que la parcelle (si pente parcelle >30%)

(val_pt_topo <- extract(st[[c("ptes","aspsx","aspsy")]],
                        p.pendage %>% st_buffer(50),
                        fun = mean, na.rm=TRUE) %>%
   as.data.frame() %>%
   mutate(asp = atan(aspsy/aspsx)/pi*180))
p.pendage$asp <- val_pt_topo$asp

val_topo_prf <- extract(st,zn_prf %>% st_union() %>% as_Spatial(), fun=mean)[1,]

(asp_prf <- atan(val_topo_prf["aspsy"]/val_topo_prf["aspsx"])/pi*180)

if(atan(val_topo_prf["ptes"])>0.3){
  p.pendage.v <- p.pendage %>% filter(abs(asp-asp_prf)<45) %>%
    mutate(PENDAGE = replace(PENDAGE,PENDAGE == 999,NA))
}else{
  p.pendage.v <- p.pendage
}

Les groupes qui ne s'apparentent à aucun des point d'observation restants sont éliminés.

# clump le plus prche des valeurs moyennes des points restants

diff_azx <- sum_clu$aspx - mean(cos(p.pendage.v$AZIMUT/180*pi))
diff_azy <- sum_clu$aspy - mean(sin(p.pendage.v$AZIMUT/180*pi))
sum_clu$diff_az <- abs(diff_azx * diff_azy)^.5/pi*180

moy_pte_pts <- mean(p.pendage.v$PENDAGE, na.rm=TRUE)
sum_clu$diff_pte <- abs(moy_pte_pts - sum_clu$pte_m *180 / pi)

# élimine les clumps très différents

sum_clu_select <- sum_clu

if("diff_az" %in% colnames(sum_clu)){
  sum_clu_select <- sum_clu_select %>%
    filter((diff_az/mean(diff_az))<2)
}

if("diff_pte" %in% colnames(sum_clu)){
  sum_clu_select <- sum_clu_select %>%
    filter((diff_pte/mean(diff_pte))<2)
}

choix <- sum_clu_select$m_cl[1]

Le plan de courbure du groupe restant le plus représenté est choisi.

# création du plan supposé de la strate géologique

vplan_choix <- vplan[which(cl_hd %in% choix),]
clump_choix <- names(cl_hd[cl_hd %in% choix]) %>% as.numeric()

cnt <- v %>% na.omit() %>% count(clump) %>%
  filter(clump %in% clump_choix)

st$plan <- st$mnt
values(st$plan) <- apply(vplan_choix,2, weighted.mean, w = cnt$n)

qm(st$plan)

Création du raster de prédiction

Les réserves profondes sont estimées selon la pente existante entre le plan de courbure et la surface topographique. Si la pente est nulle, l'indice de réserve est égal à 0. Si la pente est > 100%, ou si la courbure locale est supérieures à 0.02, l'indice est de 100.

# pente du terrain par rapport à ce plan

st$h <- st$mnt - st$plan
st$ph <- terrain(st$h,"slope")

rr <-tan(st$ph)*100
rr[rr>100] <- 100
rr[abs(st$cur)>.02] <-100

rrc <- rr %>% crop(zn_prf %>% st_buffer(100) %>% as_Spatial())
stc <- st %>% crop(zn_prf %>% st_buffer(100) %>% as_Spatial())

qm(rrc) + tmap::tm_shape(zn_prf) + tmap::tm_borders(col="black")

Vue 3D:

options(rgl.printRglwidget = TRUE) # si serveur
rgl::rgl.clear()
rasterVis::plot3D(stc$mnt,drape=rrc,col=hcl.colors)

ONF-AURA/geoDalle documentation built on Jan. 14, 2021, 3:17 a.m.