Skip to main content

Cours sur les plans d’échantillonnage statistique dans l’espace pour l’étude des sols

3 Données

3.1 Statistiques globales

Pour la partie s’interessant aux statistiques globales, nous basons nos exercices sur une réalité simulée d’une valeur de carbone dans un parcelle agricole. Ce champ a été simulé en utilisant des techniques géostastitques.

L’objet correspond à un tableau comportant 3 colonnes. les coordonnées avec s1 et s2 les valeurs simulées sim1

Chaque coordonnée correspond au centre d’un pixel

champ <- read.csv('data/simulatedField_1Exp25.csv')

head(champ)
##    s1  s2     sim1
## 1 0.5 0.5 7.184881
## 2 1.5 0.5 7.366762
## 3 2.5 0.5 7.417718
## 4 3.5 0.5 7.305573
## 5 4.5 0.5 7.373427
## 6 5.5 0.5 7.161203

Il est donc facile de transformer le tableau en une grille de type SpatialPixelsDataframe. Pour cela, on passe par une structure de type SpatialPointDataframe dans un premier temps. ce format ancien est nécessaire pour le package spcosa. Il faudra changer de format des raster selon les packages. C’est notament le cas pour tmap qui ne lit que du terra

champSP <- champ
# transformer en un vecteur de points
coordinates(champSP) <- c('s1','s2') 
# transformer en raster
gridded(champSP) <- TRUE

Voici la structure de l’objet de type sp.

# structure de l'objet
str(champSP)
## Formal class 'SpatialPixelsDataFrame' [package "sp"] with 7 slots
##   ..@ data       :'data.frame':  10000 obs. of  1 variable:
##   .. ..$ sim1: num [1:10000] 7.18 7.37 7.42 7.31 7.37 ...
##   ..@ coords.nrs : num(0) 
##   ..@ grid       :Formal class 'GridTopology' [package "sp"] with 3 slots
##   .. .. ..@ cellcentre.offset: Named num [1:2] 0.5 0.5
##   .. .. .. ..- attr(*, "names")= chr [1:2] "s1" "s2"
##   .. .. ..@ cellsize         : Named num [1:2] 1 1
##   .. .. .. ..- attr(*, "names")= chr [1:2] "s1" "s2"
##   .. .. ..@ cells.dim        : Named int [1:2] 100 100
##   .. .. .. ..- attr(*, "names")= chr [1:2] "s1" "s2"
##   ..@ grid.index : int [1:10000] 9901 9902 9903 9904 9905 9906 9907 9908 9909 9910 ...
##   ..@ coords     : num [1:10000, 1:2] 0.5 1.5 2.5 3.5 4.5 5.5 6.5 7.5 8.5 9.5 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : NULL
##   .. .. ..$ : chr [1:2] "s1" "s2"
##   ..@ bbox       : num [1:2, 1:2] 0 0 100 100
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:2] "s1" "s2"
##   .. .. ..$ : chr [1:2] "min" "max"
##   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
##   .. .. ..@ projargs: chr NA

Il est assez simple de faire une carte avec le package tmap mais il faut transformer le raster dans le format du package terra comme ceci

champSP.terra <- rast(champSP)
tm_shape(champSP.terra)  +
  tm_raster() +
tm_layout(legend.outside = TRUE)
Notre réalité: un champ aléatoire simulé avec un variogramme de type exponentiel

Figure 3.1: Notre réalité: un champ aléatoire simulé avec un variogramme de type exponentiel

L’histogramme de la propriété ressemble à ceci:

ggplot(champ, 
       aes(x=sim1)) + 
  geom_histogram(bins = 50)
Histogramme de la population

Figure 3.2: Histogramme de la population

La moyenne vraie de la population vaut ainsi

Moy = mean(champ$sim1)

Moy
## [1] 5.298602

3.2 Cartographie

Pour illuster les différents plans possible pour la cartographie, nous prenons de jeus de données en exemple

3.2.1 Contour d’une parcelle

le contour d’une parcelle est fourni par le package spcosa. Sa forme est complexe et permet de bien illustrer les caractéristiques des plans d’échantillonnage proposés.

shpFarmsum.sf <- read_sf(dsn = system.file("maps", 
                                        package = "spcosa"), 
                      layer = "farmsum")

shpFarmsum <- as(shpFarmsum.sf, "Spatial") # conversion sp

tm_shape(shpFarmsum) + 
  tm_borders()
Contour d'une parcelle agricole

Figure 3.3: Contour d’une parcelle agricole

3.2.2 Covariables rasters en Australie

Un ensemble de 5 rasters sont disponibles pour une parcelle située en Australie dans un zone viticole appelée ‘Hunter vallée’.

  • Altitude
  • l’orientation de la pente
  • pente
  • indice de topographie
  • l’indice de végétation normalisé
#Read data with coordinates and other attributes of fine grid (discretization of study area)
grdHunterValley <- readRDS(file="data/grdHunterValley.RDS")

grdHunterValley2 <- grdHunterValley
coordinates(grdHunterValley2) <- c('Easting','Northing')
gridded(grdHunterValley2) <- TRUE

tm_shape( rast(grdHunterValley2) ) + 
  tm_raster(
    col = c("elevation_m","slope_deg","cti","ndvi")
            )  +
tm_layout(legend.outside = TRUE)
Covariables raster en Australie

Figure 3.4: Covariables raster en Australie