Spatialisation

Code

require(conflicted)
require(knitr)
require(magrittr)
require(mapview)
require(pander)
require(plotly)
require(spatstat)
require(tidyverse)
conflict_prefer("set_names", "purrr")

Code

read_chunk("code/spatial/bornage_irstea.r")
read_chunk("code/spatial/bornage_malo2.r")

Code

library(readODS)
library(sf)

Linking to GEOS 3.9.1, GDAL 3.3.2, PROJ 8.0.1; sf_use_s2() is TRUE

Géolocalisation de la grille d’inventaire

Initialement et lors des nouveaux inventaires, les arbres sont positionnés de manière relative au sein des quadrats de 10$$10 m. Pour pouvoir positionner les arbres à l’échelle des parcelles d’inventaire, il faut disposer de la position des piquets

Une partie des piquets délimitant les quadrats a été positionnée au tachéomètre (Bornage IRSTEA), les piquets manquants, car non accessibles avec l’appareil, l’ont été fin 2024 par triangulation simple, par mesure de coordonnées polaires à partir de piquets déjà positionnés.

Bornage IRSTEA

Import des données

Code

## Positions piquets (format IRSTEA)
bornage <- "SWIO_RUN_MALO1_2_bornage_20150427.ods" %>% 
  str_c("data/bornage/", .) %>% 
  read_ods(sheet = 1)
tmp <- "SWIO_RUN_MALO3_bornage_20150404.ods" %>%
    str_c("data/bornage/", .) %>%
  read_ods(sheet = 1) %>% 
  transform(sub_id = as.character(sub_id))

bornage %<>% bind_rows(tmp)

Apercu des données :

Code

head(bornage)

# A tibble: 6 × 6
  piquet     x     y     z sub_id plot_id
   <dbl> <dbl> <dbl> <dbl> <chr>  <chr>  
1      1  860. 1835.  239. A0     MALO1  
2      2  862. 1845.  240. B0     MALO1  
3      3  864. 1854.  242. C0     MALO1  
4      4  865. 1864.  244. D0     MALO1  
5      5  867. 1874.  247. E0     MALO1  
6      6  868. 1885.  251. F0     MALO1

Calcul des angles

Code

# Angle MALO1-2 depuis bornage
theta <- bornage %>% 
  dplyr::filter(! is.na(sub_id)) %>% 
  with((y[sub_id == "A9"] - y[sub_id == "A0"]) / (x[sub_id == "A9"] - x[sub_id == "A0"])) %>% 
  atan()

# Angle MALO3 depuis bornage
theta2 <- bornage %>% 
  dplyr::filter(! is.na(sub_id)) %>% 
  with((y[sub_id == "201"] - y[sub_id == "11"]) / (x[sub_id == "201"] - x[sub_id == "11"])) %>% 
  atan()

# Angle MALO1-2 depuis GPS
#theta2 <- atan(with(pspdf, (yutm[Name == "A10"] - yutm[Name == "A0"])/(xutm[Name == "A10"] - xutm[Name == "A0"])))

Changement de repère

MALO1-2

Code

bornage %<>% 
  sf::st_as_sf(coords =  c("x", "y"))
born <- bornage %>% 
  dplyr::filter(plot_id %in% c("MALO1", "MALO2")) %>% 
  add_coords(suffix = "utm", spatial = TRUE) %>% 
  rotate_spdf(theta) %>% 
  add_coords() %>%
    transform(x = x - min(x), y = y - min(y))

Warning: some mark values are NA in the point pattern x

Code

#  class(born) <- class(bornage)
# born %<>%
#   add_coords() %>%
#   transform(x = x - min(x), y = y - min(y))

MALO3

Code

tmp <- bornage %>% 
  dplyr::filter(plot_id == "MALO3") %>% 
  rotate_spdf(theta2) %>%
    add_coords() %>%
    transform(x = x - min(x), y = y - min(y))

Warning: some mark values are NA in the point pattern x

Tableau complet

Code

born %<>% bind_rows(tmp)

Cartographie

Fonctions d’appui

Code

# Rotation du semis de points
rotate_spdf <- function(x, theta)
{
    x %<>%
        as.ppp() %>%
        spatstat.geom::rotate(- theta, centre = "bottomleft") %>%  
        st_as_sf() %>% 
    dplyr::filter(! is.na(plot_id)) # bug: rétroconversion depuis ppp, ajout ligne POLYGON ?
    names(x) %<>% str_replace("marks.", "")
    x
}

# Ajouts des coordonnées
add_coords <- function(x, spatial = FALSE, suffix = "")
{
    res <- x %>% 
        sf::st_coordinates() %>%
        data.frame()
  n <- str_c(c("x", "y"), suffix) 
    names(res)[1 : 2] <- n
    res <- bind_cols(x, res)
  if(spatial) res %<>% st_as_sf(coords = n)
    res 
}

Compléments terrain

Chargement

Code

library(readODS)
library(sf)

Import des données

Code

# Import des relevés terrain
dat <- read_ods("data/bornage/20241203_relevés bornage MALO2.ods", sheet = "correc") %>% 
  mutate(x = NA, y = NA)


# Import des position de référence (bornage IRSTEA)
# born <- st_read("data/sig/bornage_MALO_IRSTEA.gpkg")
# names(born) %<>% str_to_lower()

# Conversion azimuth (A) > trigonométrie
dat %<>% 
  mutate(azim = ((360 - azim) + 90) %% 360, # chgt d'origine et sens
         azimrad = azim * (pi / 180)) # radians

Fonction de calcul des coordonnées

Code

coord_rel_pol <- function(df, group = c("cible_sub_id", "cible_dir")) {
  df %<>% 
  mutate(x = x_ref + dis_hor * cos(azimrad),
           y = y_ref + dis_hor * sin(azimrad)) %>% 
  dplyr::select(- ends_with("_ref")) %>% 
  group_by(df[,!! group]) %>% 
  summarise(x = mean(x, na.rm = FALSE),
            dx = diff(range(x, na.rm = FALSE)),
            y = mean(y, na.rm = FALSE),
            dy = diff(range(y, na.rm = FALSE)),
            n = n()) %>%
  ungroup()
  df
}

Tableau des cibles

Nombre de piquets à positionner :

Code

## Nombre de cibles
dat %>% 
  distinct(cible_sub_id, cible_dir) %>% 
  nrow()

[1] 24

Première fusion

Utilisation des coordonnées des piquets déjà positionnés (manip IRSTEA)

Code

# first round : match with coordinates
# from IRSTEA (object 'born')
df <- dat %>%
  inner_join(born, by = c("ref_sub_id" = "sub_id"), suffix = c("", "_ref")) %>% 
  coord_rel_pol()

`summarise()` has grouped output by 'cible_sub_id'. You can override using the
`.groups` argument.

Références manquantes après première fusion :

Seconde fusion

Positionnement à partir des cibles nouvellement positionnées :

Code

# second round : match with estimated coordinates
tmp <- dat %>%
  anti_join(df, by = c("cible_sub_id", "cible_dir")) %>% 
  left_join(df, by = c("ref_sub_id" = "cible_sub_id", "ref_dir" = "cible_dir"), suffix = c("", "_ref")) %>% 
  dplyr::filter(! is.na(x_ref)) %>%
  coord_rel_pol()

`summarise()` has grouped output by 'cible_sub_id'. You can override using the
`.groups` argument.

Code

df %<>% bind_rows(tmp)

Références manquantes après seconde fusion :

[1] 0

Cartographie

Code

df %>% 
  sf::st_as_sf(coords =  c("x", "y")) %>% 
  mapview()