Setting the main path directories

path0    = "."
path_in  = paste0(path0,"/data-raw")
path_out = paste0(path0,"/data-output")

Import LUCAS data

We have three shapefiles representing the Lucas data sample points for the years 2009, 2015 and 2018 for the European Union member states. The points corresponding to those of the Walloon region are extracted using a shapefile of the Walloon boundaries.

As mentioned above, the attribute table of the LUCAS 2018 shapefile does not contain all the variables of interest, so a join is made with the excel file.

if (op.import.data) {
# Import and read LUCAS data
f_2009 <- paste0(path_in,"/SoilAttr_LUCAS2009/SoilAttr_LUCAS_2009.shp")
f_2015 <- paste0(path_in,"/LUCAS_SOIL_2015/LUCAS_Topsoil_2015_20200323.shp")
f_2018 <- paste0(path_in,"/LUCAS_SOIL_2018/LUCAS-SOIL-2018 .shp")

lucas_eu_2009 <- st_read(f_2009,stringsAsFactors = F)
lucas_eu_2015 <- st_read(f_2015,stringsAsFactors = F)
lucas_eu_2018 <- st_read(f_2018,stringsAsFactors = F)

# Import and read the limits of Walloon
f_limit <- paste0(path_in,"/wallonie/provinces.shp")
wallonie <- st_read(f_limit, stringsAsFactors = F)

# Define LUCAS Walloon data

#-- Convert coordinates L72 to wgs_84
wallonie_wgs84 <- st_transform(wallonie, st_crs(lucas_eu_2015))


#-- Intersection
lucas_wal_2009 <- st_intersection(lucas_eu_2009,wallonie_wgs84)
lucas_wal_2015 <- st_intersection(lucas_eu_2015,wallonie_wgs84)
lucas_wal_2018 <- st_intersection(lucas_eu_2018,wallonie_wgs84)


# Join
f_lucas_eu_2018 <- paste0(path_in,"/LUCAS_SOIL_2018/LUCAS-SOIL-2018.xls")
## excel_lucas_eu_2018 <- read.xlsx2(f_lucas_eu_2018,sheetIndex="LUCAS-SOIL-2018",
##                             header=TRUE,colClasses=c("character",rep("numeric",9),
##                                                      rep("character",8),rep("numeric",2),
##                                                      "character","numeric",rep("character",5)))

excel_lucas_eu_2018 <- readxl::read_excel(f_lucas_eu_2018,
sheet = "LUCAS-SOIL-2018",
col_types = c("text",rep("numeric",9),
             rep("text",8),rep("numeric",2),
             "text","numeric",rep("text",5))
)


lucas_wal_2018 <- left_join(
    lucas_wal_2018, 
    excel_lucas_eu_2018, 
    by = c("POINTID" = "POINTID")
    )
}

Import REQUASUD data

The soil analysis data of Wallonia are imported and attached to an excel file containing the name of the land cover code according to REQAUSUD. The following code is the one that import the “fuzzyfied” data from REQUASUD. But, for this report, we used the original data set, only internally available. The file that is sourced has the same structure of the “fuzzyfied” one.

if (op.import.data) {
if (op.data.file == "fuzzyfied") {
# Import and read REQUASUD data
f_REQ <- paste0(path_in,"/REQ/extraction_a10_2022_corrige.csv")
REQUASUD.1=read.csv2(f_REQ,header = TRUE, dec=",", stringsAsFactors = F,
                   colClasses = c("numeric",
                                  "numeric",
                                  "character","character",
                                  "numeric",
                                  "character","character",
                                  rep("numeric",6),"character",
                                  rep("numeric",5),"character","character",
                                  "numeric","character","character",
                                  rep("numeric",3),"character","character",
                                  "numeric","character","character"
                                  ),
                    na.strings = "NULL", encoding = "latin1")

## import landcover classes
df.lc.class <- read.csv2("./data-raw/REQ/occupation.csv", encoding = "latin1")

REQUASUD.2 <- REQUASUD.1 %>%
left_join(df.lc.class, by = c("occup.sol" = "code"))

REQUASUD.2$sol_pheau[REQUASUD.2$sol_pheau == 0.5] <- NA
#print(qplot(REQUASUD$sol_pheau,geom = "boxplot"))

#-- Convert df to sf
req <- st_as_sf(REQUASUD.2,coords=c("sol_x","sol_y"), crs=31370)
}
if (op.data.file == "clean") {
    source("../walloonlucas_clean/R/clean_data_import.R")
}
}

Coordinates conservation for future representations

This step is just to keep the coordinates of the points of different projection systems in a dataframe for future use. Three systems were chosen: - the World Geodetic System 1984 (used in GPS), - the Belgian Lambert 72 and - the European ETRS89-extended system.

if (op.import.data) {
# LUCAS

#--2009
#----SCR
lucas_wal_2009_l72 <- st_transform(lucas_wal_2009, 31370)
lucas_wal_2009_3035 <- st_transform(lucas_wal_2009, 3035)

#----conservation
lucas_wal_2009$Long_wgs84 <- st_coordinates(lucas_wal_2009)[,1]
lucas_wal_2009$Lat_wgs84 <- st_coordinates(lucas_wal_2009)[,2]

lucas_wal_2009$X_l72 <- st_coordinates(lucas_wal_2009_l72)[,1]
lucas_wal_2009$Y_l72 <- st_coordinates(lucas_wal_2009_l72)[,2]

lucas_wal_2009$X_3035 <- st_coordinates(lucas_wal_2009_3035)[,1]
lucas_wal_2009$Y_3035 <- st_coordinates(lucas_wal_2009_3035)[,2]

#--2015
#----SCR
lucas_wal_2015_l72 <- st_transform(lucas_wal_2015, 31370)
lucas_wal_2015_3035 <- st_transform(lucas_wal_2015, 3035)

#----conservation
lucas_wal_2015$Long_wgs84 <- st_coordinates(lucas_wal_2015)[,1]
lucas_wal_2015$Lat_wgs84 <- st_coordinates(lucas_wal_2015)[,2]

lucas_wal_2015$X_l72 <- st_coordinates(lucas_wal_2015_l72)[,1]
lucas_wal_2015$Y_l72 <- st_coordinates(lucas_wal_2015_l72)[,2]

lucas_wal_2015$X_3035 <- st_coordinates(lucas_wal_2015_3035)[,1]
lucas_wal_2015$Y_3035 <- st_coordinates(lucas_wal_2015_3035)[,2]

#--2018
#----SCR
lucas_wal_2018_l72 <- st_transform(lucas_wal_2018, 31370)
lucas_wal_2018_3035 <- st_transform(lucas_wal_2018, 3035)

#----conservation
lucas_wal_2018$Long_wgs84 <- st_coordinates(lucas_wal_2018)[,1]
lucas_wal_2018$Lat_wgs84 <- st_coordinates(lucas_wal_2018)[,2]

lucas_wal_2018$X_l72 <- st_coordinates(lucas_wal_2018_l72)[,1]
lucas_wal_2018$Y_l72 <- st_coordinates(lucas_wal_2018_l72)[,2]

lucas_wal_2018$X_3035 <- st_coordinates(lucas_wal_2018_3035)[,1]
lucas_wal_2018$Y_3035 <- st_coordinates(lucas_wal_2018_3035)[,2]


#REQUASUD

#--SCR
req_wgs84 <- st_transform(req, 4326)
req_3035 <- st_transform(req, 3035)


#--conservation
req$Long_wgs84 <- st_coordinates(req_wgs84)[,1]
req$Lat_wgs84 <- st_coordinates(req_wgs84)[,2]

req$X_l72 <- st_coordinates(req)[,1]
req$Y_l72 <- st_coordinates(req)[,2]

req$X_3035 <- st_coordinates(req_3035)[,1]
req$Y_3035 <- st_coordinates(req_3035)[,2]
}
rm(list=c("REQUASUD.1","REQUASUD.2"))

Data processing and dataframes

All datafames useful for the comparison are finally produced and saved here.

if (op.import.data) {

# LUCAS
LUCAS2009 <- st_drop_geometry(lucas_wal_2009)
LUCAS2015 <- st_drop_geometry(lucas_wal_2015)
LUCAS2018 <- st_drop_geometry(lucas_wal_2018) %>%
dplyr::select(-"OC (20-30 cm)",
 -"CaCO3 (20-30 cm)",
 -"Ox_Al",
 -"Ox_Fe")

# REQUASUD
REQUASUD <- st_drop_geometry(req)
}

Saving and loading imported and pre-processed data

For lowering the computing time of loading data, after a first run of the previous code, we can save all objects in a .RData file. The next run can start here, loading the objects only by reading this saved file (name of the file is a rmarkdown parameter).

if (op.import.data) {
v.objects <- ls()[!(ls() %in% c("params"))]
save(list = v.objects,
     file = paste0("./data-raw/",op.data.file,".RData")
    )
} else {
load(paste0("./data-raw/",op.data.file,".RData"))
}

LUCAS

Soil properties available in LUCAS and their measurement methods from LUCAS Topsoil Survey

The soils are sampled for a 0 to 20 cm increment depth (0-30 cm for the last version).

The data and methods available from LUCAS datasets are:

lucas.methods <- read.csv2(paste0(path_in,"/lucas-data-methods.csv"))

kable(lucas.methods, caption = "Soil properties methods in LUCAS Survey") %>%
  kable_styling(bootstrap_options = "striped", 
  "bordered", 
  full_width = T)

In the next table, we present the 20 first lines of the main LUCAS dataset (2018).

head(LUCAS2018, n = 20) %>%
  datatable(
  options = list(
  pageLength = 5, autoWidth = TRUE, 
  searchHighlight = FALSE), 
  class = 'cell-border stripe', filter = 'bottom'
  )

REQUASUD

The soil properties available at lucas data level are mostly those measured at Requasud network level. However the units and methods may be different.

The methods defined in the Requasud data set for each soil property are quite variable. They are presented in a table in appendix.

Here we present a list of data available in the Requasud dataset:

• pH in KCl (sol_phkcl) [-]
• pH in H2O (sol_pheau) [-]
• organic carbon content (sol_carbone) [% m.s]
• phosphorus available (sol_phosphore) [mg/100 g TS]
• potassium available (sol_potassium) [mg/100 g TS]
• magnesium available (sol_magnesium) [mg/100 g TS]
• calcium available (sol_calcium) [mg/100 g TS]
• nitrogen content (sol_azote) [mg/g TS]
• cation exchange capacity (sol_cec) [meq/100 g TS]
• particle size distribution (“sol_argile”, “sol_limon” and “sol_sable” are for clay, silt and sand content) [%]
• carbonate content (sol_carbonate) [%]

In the next table, we present the 20 last lines of the REQUASUD dataset. We removed the geodata.

tail(REQUASUD, n = 20) %>%
    dplyr::select(!c("Long_wgs84","Lat_wgs84","X_l72","Y_l72","X_3035","Y_3035")) %>%
  datatable(
  options = list(
  pageLength = 5, autoWidth = TRUE, 
  searchHighlight = FALSE), 
  class = 'cell-border stripe', filter = 'bottom'
  )

Only the soil properties data available in the two datasets (LUCAS and REQUASUD) are considered for the comparison in this report.

Organic carbon content

To fit the LUCAS OC unit, we multiplied by a factor of 10 to convert from % (Requasud) to [g/kg] (Lucas).

REQUASUD$conv_sol_carbone <- REQUASUD$sol_carbone * 10

Phosphorus content

To convert from mg/100g (Requasud) to mg/kg (Lucas), we will multiply the Requasud data by a 10 factor.

REQUASUD$conv_sol_phosphore <- REQUASUD$sol_phosphore * 10

Potassium content

The same convention as in the case of phosphorus is applied here.

REQUASUD$conv_sol_potassium <- REQUASUD$sol_potassium * 10

Number of sites and densities

We provide in the following table the number of sites for both LUCAS and REQUASUD datasets over the study area and their associated densities. The area of Wallonia has been set to 16 845 km².

Note that the comparison of the particles size is done with the data from Lucas (2009 and 2015) and all the data from Requasud. While for the other soil properties, we considered the year 2018.

area <- 16845

# density and number of site table
density_table <- cbind.data.frame(Dataset = c("LUCAS.2015 (particule size)","REQUASUD (particule size)",
                                              "LUCAS.2018","REQUASUD.2018"),
  Nb.site = c(nrow(data.granuloLUCAS),nrow(data.granuloREQUASUD),
                                   nrow(data.chemicalLUCAS),nrow(data.chemicalREQUASUD)
                                   ),
                 Density = c(ceiling((area)/nrow(data.granuloLUCAS)),
                             ceiling((area)/nrow(data.granuloREQUASUD)),
                                     ceiling((area)/nrow(data.chemicalLUCAS)),
                             ceiling((area)/nrow(data.chemicalREQUASUD))
                             )
  )%>% 
  #PRENDRE LE NOMBRE DE POINT PAR KM2 A CHANGER
  kbl(caption = "Number of sites of LUCAS and REQUASUD datasets and their density
      (1 obs for X km²) over study area") %>%
  kable_styling(bootstrap_options = "striped", full_width = F)

density_table

Spatial distribution representation

f.wl.spatial(wallonie,
             data.granuloREQUASUD,
             data.granuloLUCAS,
             "Long_wgs84",
             "Lat_wgs84" ,
             coord.crs = 4326,
             data.used = "Data_Particle.size"
)

f.wl.spatial(wallonie,
             data.chemicalREQUASUD,
             data.chemicalLUCAS,
             "Long_wgs84",
             "Lat_wgs84" ,
             coord.crs = 4326,
             data.used = "Data_Chemical.property"
)

Interactive maps (LUCAS only)

m_lucas <- mapview(lucas.sf.l72, col.regions = viridis(3)[1], cex = 2, lwd = 0.5)

• Individual maps of LUCAS and SIMS spatial points maps.

m_lucas

Sampling point density

We can calculate sampling point density in France for LUCAS in a radius of 10km.

## densité points
w <- as.owin(wallonie) ## creating a window
w.km <- rescale.owin(w, s = 1000) ## change the units in kilometers

## lucas point density
s.lucas  <- as.ppp(lucas_wal_2018_l72) ## points in projected system
marks(s.lucas) <- NULL
s.lucas <- rescale.ppp(s.lucas, s = 1000) # convert units in km
Window(s.lucas) <- w.km ## making the same window

density.lucas <- density.ppp(s.lucas, sigma=10) # calculating point density for LUCAS with 16 km bandwidth
plot(density.lucas, main="Point density of LUCAS samples", las=1)
contour(density.lucas, add=TRUE)

## Requasud point density
s.req  <- as.ppp(req) ## points in projected system
marks(s.req) <- NULL
s.req <- rescale.ppp(s.req, s = 1000) # convert units in km
Window(s.req) <- w.km ## making the same window

density.req <- density.ppp(s.req, sigma=10) # calculating point density for LUCAS with 16 km bandwidth
plot(density.req, main="Point density of REQ samples", las=1)
contour(density.req, add=TRUE)

Main observation, and consequences for the whole following report

We observe that the spatial distribution of the data are not the same in the two data sets.

For the particule size analyses (PSA), the LUCAS points are equally distributed in the whole Walloon region while the REQUASUD points are concentrated around the province of Brabant wallon, whose main laboratory (La Hulpe) is the one in charge with this analysis in the REQUASUD network. This region is also know as the main one for intensive arable farming, in the loamy belt of Belgium.

For the chemical properties, we face also with two contrasted situations : the LUCAS datasets hold a limited number of sampling points, but equally distributed in the Walloon region. On the other side, the REQUASUD datasets hold a huge number of point, but with a high concentration of points in the intesive arable cropping systems of four provinces : East of Hainaut, North of Namur, Brabant wallon and west of Liège. We see also some spots in the grassland farming systems in Luxembourg province and East of the Liège province.

This contrasted situation find its origin in the two contrasting “goals” of these data sets : on the one hand, LUCAS aims at monitor the soil of Europe in general, and in Wallonia for this report ; on the other hand, REQUASUD is a service proposed to farmers, and is done in the framework on some administrative demands or on a voluntary basis.

This contrasted situation has logical consequence on the data analyses that will be done here after. It will be an important source of variablity between the datasets in following comparisons.

Full datasets

LUCAS2018$LandCover <- as.factor(LUCAS2018$LC0_Desc)
req$LandCover <- as.factor(req$LC0_Desc)

# Extracting the sites per land cover classes 
    myTabLC <- dplyr::left_join(
    as.data.frame(table(LUCAS2018$LandCover)),
    as.data.frame(table(req$LandCover)),
    by=c("Var1"="Var1")
    )

colnames(myTabLC) <- c("LandCover", "nb.LUCAS2018", "nb.REQUASUD.full")

# Computing percentage
myTabLC$nb.LUCAS2018.perc <- round((myTabLC$nb.LUCAS2018/nrow(LUCAS2018))*100,2)
myTabLC$nb.REQUASUD.2018.perc <- round((myTabLC$nb.REQUASUD.full/nrow(req))*100,2)

# Creating the table
SitesLU <- myTabLC %>%
  kbl(caption = "Number of sites by Land Cover Classes") %>%
  kable_styling(bootstrap_options = "striped",
                full_width = F)

SitesLU

Filtered datasets

As REQUASUD is a service mainly proposed to Walloon farmers, two main land covers are dominants : cropland and grassland. Bareland is an uncommon land cover and is considered, in this report, as a farming plot (crop or grassland) that is temporary not covered. For trying to have as comparable as possible datasets, we filtered the data on these three land covers.

For chemical data at least, the REQUASUD dataset is also filtered based on the year (only 2018) and on the sampling depth (15 cm and 25 cm), close to the LUCAS sampling protocol.

This is why the numbers of sampling points in the second table are lower for the REQUASUD column.

data.chemicalLUCAS$LandCover <- as.factor(data.chemicalLUCAS$LC0_Desc)
data.chemicalREQUASUD$LandCover <- as.factor(data.chemicalREQUASUD$LC0_Desc)

# Extracting the sites per land cover classes
    myTabLC <- dplyr::left_join(
    as.data.frame(table(data.chemicalLUCAS$LandCover)),
    as.data.frame(table(data.chemicalREQUASUD$LandCover)),by=c("Var1"="Var1")
    )

colnames(myTabLC) <- c("LandCover", "nb.LUCAS2018", "nb.REQUASUD2018)")

# Computing percentage
myTabLC$nb.LUCAS2018.per <- ceiling((myTabLC$nb.LUCAS2018/nrow(data.chemicalLUCAS))*100)
myTabLC$nb.REQUASUD.2018.per <- ceiling((myTabLC$nb.REQUASUD2018/nrow(data.chemicalREQUASUD))*100)

# Creating the table
SitesLU <- myTabLC %>%
  kbl(caption = "Number of sites by Land Cover Classes") %>%
  kable_styling(bootstrap_options = "striped",
                full_width = F)

SitesLU

We provide also the graphical representation.

# Vertical format table for graphical representation using ggplot2
TableCover  <- 
  cbind.data.frame(
    LandCover = c(as.factor(myTabLC$LandCover),
                  as.factor(myTabLC$LandCover)),
    Dataset = as.factor(c(rep("REQUASUD",3),
                          rep("LUCAS",3))),
    pc.site =c(myTabLC[,5],
               myTabLC[,4])
    )

# Graphical representation
SitesLU.graph <- ggplot(data=TableCover,
                        aes(x=LandCover, 
                            y=pc.site,
                            fill=Dataset 
                            )
                        ) +
   geom_bar(stat="identity",
            position=position_dodge()
            ) +
   scale_fill_viridis(direction=-1,discrete = TRUE) +
  ylab("Percentage of sites") + 
  theme_minimal()

SitesLU.graph

Repartition of LUCAS and REQUASUD sites according to the land cover classification

Observation

The two major soil units of Wallonia are Orthic Luvisols (Lo, approx. 70.5%), mainly in cropland, and Dystric Cambisols (Bd, approx. 19.5%), mainly in grassland regions.

Comparison

Graphical representation of the difference between distribution of dominant soil classes in Beglium/Wallonia (i.e. Reference frequency) and the soil classes frequencies in LUCAS and REQUASUD data sets.

dom.soil <- dom.soil.freq %>%
  filter(! Dominant_soil %in% c("555","444","333","222", "111","")) %>%
  filter(Frequency > 0.1) %>% ## filter classes with frequency is above 0.1 % in the FAO 
  rowwise() %>%
  mutate(lucas = freq.lucas - Frequency,
         requasud = freq.req - Frequency) %>% ## difference between soil classes frequencies
  pivot_longer(., cols = c("lucas","requasud"), values_to = "Diff.frequency", names_to = "Program") %>%
  unique() %>%
  arrange(Program, Diff.frequency)
dom.soil$Dominant_soil <- factor(dom.soil$Dominant_soil, levels = unique(dom.soil$Dominant_soil))


dom.soil %>%
  group_by(Program) %>%
  dplyr::summarise(median = median(Diff.frequency, na.rm = T),
            min = min(Diff.frequency, na.rm = T),
            max = max(Diff.frequency, na.rm = T)) %>%
  kbl(caption = "Statistics of differences between global frequencies of dominant soil classes and datasets", digits = 3) %>%
  kable_styling(bootstrap_options = "striped", full_width = F)

dom.soil %>%
  ggplot(., aes(x = Dominant_soil, y = Diff.frequency, color = Program)) +
  geom_point(stat = "identity", size = 3) +
  geom_segment( aes(x=Dominant_soil, xend=Dominant_soil, y=0, yend=Diff.frequency), color="grey") +
  coord_flip() +
  geom_vline(xintercept = 0) +
  theme_bw() +
  theme(legend.position = "none")+
  facet_wrap(~ Program) +
  labs(title = "Difference between frequencies of dominant soil classes",
       y = "Reference frequency - Monitoring frequency",
       x = "Dominant soil classes")

High negative of positive values implies either an under- or over-representation of dominant soil classes.

Observations

Statistically, the Orthic luvisol (cropland) is overrepresented, and the Dystric cambisol is underrepresented in the two datasets, but more clearly in the REQUASUD data set. It is linked to the goal of the requasud service, more popular in the cropland area of Wallonia.

The reasons of the choice

As the Lucas 2018 particle size data is not available, we will use the particle size data from other available years for comparison.

For the Requasud data, we noted that the particle size data started to be measured recently (see figure below). Since particle size evolves slowly on the human time scale, we used all available data at the Requasud level.

knitr::include_graphics("./data-raw/evolution_granulo_anim.gif")

Animated map showing the clay content mesurement from 2008 to 2021

Central tendency

We estimate central tendency of LUCAS and REQUASUD dataset

f.wl.stat.table(summary_granuloLUCAS,summary_granuloREQUASUD,'central')

Measures of spread

In this section, we estimate the range, the quantiles (25% and 75%), the standard deviation, the skewness and kurtosis. These indicators will provide information to qualify the distribution of a soil property either spread or not.

f.wl.stat.table(summary_granuloLUCAS,summary_granuloREQUASUD,'spread')

Comparison graphical representations

df.granulo.awc.lucas <- f.soil.awc.ptf(data.granuloLUCAS)
g.tt.lucas <- f.soil.tritext(df.granulo.awc.lucas)


TT.plot(
  class.sys = "BE.TT",
  ## tri.data = rand.text,
  tri.data = na.omit(df.granulo.awc.lucas),
  ## z.name = "OC",
  css.names = c("Clay", "Silt", "Sand"),
  main = "Soil texture triangle (BE) & data lucas",
  tri.sum.tst = FALSE
)

  geo <- TT.geo.get()


  ## Avec densité de probabilité
  kde.res <- TT.kde2d(
    geo = geo,
    tri.data = na.omit(df.granulo.awc.lucas),
    css.names = c("Clay", "Silt", "Sand"),
    tri.sum.tst = FALSE
  ) 

  TT.plot(
    class.sys = "BE.TT",
    ## tri.data = rand.text,
    tri.data = na.omit(df.granulo.awc.lucas),
    ## geo = geo,
    grid.show = FALSE,
    col = "gray",
    ## z.name = "OC",
    css.names = c("Clay", "Silt", "Sand"),
    main = "Soil texture triangle (BE) & data lucas",
    tri.sum.tst = FALSE
  )
  TT.contour(
    x = kde.res,
    geo = geo,
    main = "Probability density estimate of the texture data",
    lwd = 1,
    add = TRUE, # <<-- important
    col = "red"
  ) 

g.tt.lucas %>%
    ggtern::ggsave("./fig/tt_lucas.png", . ,
           width=9, height=7, units="in"
           )

print(g.tt.lucas)

df.granulo.awc.requasud <- f.soil.awc.ptf(data.granuloREQUASUD)
g.tt.requasud <- f.soil.tritext(df.granulo.awc.requasud,
                 colored.by = "AWC")


  ## Test the data (?? text.tol not working in TT.plot !!!
  ##
  ## TT.data.test(df.granulo.awc.requasud,
  ##         css.names = c("Clay", "Silt", "Sand"),
  ##         text.tol = 0.01
  ##         )

TT.plot(
  class.sys = "BE.TT",
  ## tri.data = rand.text,
  tri.data = df.granulo.awc.requasud,
  ## z.name = "OC",
  css.names = c("Clay", "Silt", "Sand"),
  main = "Soil texture triangle (BE) & data requasud",
  tri.sum.tst = FALSE
)

  kde.res <- TT.kde2d(
    geo = geo,
    tri.data = df.granulo.awc.requasud,
    css.names = c("Clay", "Silt", "Sand"),
    tri.sum.tst = FALSE
  ) 

  TT.plot(
    class.sys = "BE.TT",
    ## tri.data = rand.text,
    tri.data = df.granulo.awc.requasud,
    ## geo = geo,
    grid.show = FALSE,
    col = "gray",
    ## z.name = "OC",
    main = "Soil texture triangle (BE) & data requasud",
    css.names = c("Clay", "Silt", "Sand"),
    tri.sum.tst = FALSE
  )
  TT.contour(
    x = kde.res,
    geo = geo,
    main = "Probability density estimate of the texture data",
    lwd = 1,
    add = TRUE, # <<-- important
    col = "red"
  ) 

  g.tt.requasud %>%
    ggtern::ggsave("./fig/tt_requasud.png", . ,
           width=9, height=7, units="in"
           )

print(g.tt.requasud)

  requasud.texture.class <- df.granulo.awc.requasud
  requasud.texture.class$texture.class <- TT.points.in.classes( 
    tri.data    = df.granulo.awc.requasud, 
    class.sys   = "BE.TT", 
    PiC.type    = "t",
    css.names = c("Clay", "Silt", "Sand"),
    tri.sum.tst = FALSE
  )
  requasud.texture.class$dataset <- "requasud"

  lucas.texture.class <- na.omit(df.granulo.awc.lucas)
  lucas.texture.class$texture.class <- TT.points.in.classes( 
    tri.data    = na.omit(df.granulo.awc.lucas), 
    class.sys   = "BE.TT", 
    PiC.type    = "t",
    css.names = c("Clay", "Silt", "Sand"),
    tri.sum.tst = FALSE
  )
  lucas.texture.class$dataset <- "lucas"

  ## Graphical representation
  df.texture.classes <- rbind(requasud.texture.class[c("texture.class","dataset")],lucas.texture.class[c("texture.class","dataset")]) %>%
    group_by(texture.class,dataset) %>%
    summarise(n = n(),.groups="keep") %>%
    ungroup(texture.class) %>%
    mutate(freq = n / sum(n))
  SitesTT.graph <- ggplot(data = df.texture.classes,
              aes(x = texture.class,
                  y = freq,
                  fill = dataset
                  )
              ) +
    geom_bar(stat = "identity",
         position = position_dodge()
         ) +
    viridis::scale_color_viridis() +
    ylab("Percentage of sites for each sol texture classes") + 
    theme_minimal() + 
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
  
  print(SitesTT.graph)

  TT.plot(
    class.sys = "USDA.TT",
    main = "Soil texture triangle (USDA)"
  )

  requasud.texture.class$texture.class <- TT.points.in.classes( 
    tri.data    = df.granulo.awc.requasud, 
    class.sys   = "USDA.TT", 
    PiC.type    = "t",
    css.names = c("Clay", "Silt", "Sand"),
    tri.sum.tst = FALSE
  )
  requasud.texture.class$dataset <- "requasud"

  lucas.texture.class$texture.class <- TT.points.in.classes( 
    tri.data    = na.omit(df.granulo.awc.lucas), 
    class.sys   = "USDA.TT", 
    PiC.type    = "t",
    css.names = c("Clay", "Silt", "Sand"),
    tri.sum.tst = FALSE
  )
  lucas.texture.class$dataset <- "lucas"

  ## Graphical representation
  df.texture.classes <- rbind(requasud.texture.class[c("texture.class","dataset")],lucas.texture.class[c("texture.class","dataset")]) %>%
    group_by(texture.class,dataset) %>%
    summarise(n = n(),.groups="keep") %>%
    ungroup(texture.class) %>%
    mutate(freq = n / sum(n))
  SitesTT.graph <- ggplot(data = df.texture.classes,
              aes(x = texture.class,
                  y = freq,
                  fill = dataset
                  )
              ) +
    geom_bar(stat = "identity",
         position = position_dodge()
         ) +
    viridis::scale_color_viridis() +
    ylab("Percentage of sites for each sol texture classes") + 
    theme_minimal() + 
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

  SitesTT.graph

data.granuloREQUASUD %>% 
  dplyr::select(Clay, Silt, Sand) %>%
  pivot_longer(cols = Clay:Sand,
               names_to = "property",
               values_to = "value") %>%
  mutate(Dataset = "REQUASUD") %>%
  bind_rows(
    data.granuloLUCAS %>% 
      dplyr::select(Clay, Silt, Sand) %>% 
      pivot_longer(cols = Clay:Sand,
                   names_to = "property",
                   values_to = "value") %>%
      mutate(Dataset = "LUCAS") 
  ) %>%
  mutate(property = factor(property,levels = c("Clay","Silt","Sand"))) %>%
  ggplot(
    aes(x = Dataset ,
        y = value,
        fill = Dataset)) + 
  geom_boxplot() +
  labs(title="Plot of particles size by type of soil monitoring network",x="Soil monitoring", y = "Soil porperty") +
  scale_fill_viridis(direction=-1,discrete = TRUE) +
  theme_classic() +
  facet_wrap( ~ property ,
              scales = "free_y"
  )

f.wl.distr(data.granuloREQUASUD,
           data.granuloLUCAS,
           wl.param = 'Clay',
           wl.bin = 4,
           wl.xlab = "Clay content [%]"
)

f.wl.distr(data.granuloREQUASUD,
           data.granuloLUCAS,
           wl.param = 'Silt',
           wl.bin = 6,
           wl.xlab = "Silt content [%]"
)

f.wl.distr(data.granuloREQUASUD,
           data.granuloLUCAS,
           wl.param = 'Sand',
           wl.bin = 6,
           wl.xlab = "Sand content [%]"
)

Statistical inferences

We provide as well the KS test that evaluates if LUCAS and REQUASUD represent the same population. The KS is provided with a p-value.

KS.df <- cbind.data.frame(soil.property = c("Clay","Silt","Sand"),
                 KS.test = round(c(
                   as.numeric(ks.test(x=data.granuloREQUASUD$Clay,
                                      y = data.granuloLUCAS$Clay)$statistic),
                   
                   as.numeric(ks.test(x =data.granuloREQUASUD$Silt,
                                     y=data.granuloLUCAS$Silt)["statistic"]),
                   
                   as.numeric(ks.test(x = data.granuloREQUASUD$Sand,
                                      y=data.granuloLUCAS$Sand)["statistic"])),3),
                 
                 p.value = round(c(
                   as.numeric(ks.test(x =data.granuloREQUASUD$Clay,
                                      y = data.granuloLUCAS$Clay)["p.value"]),
                   
                   as.numeric(ks.test(x = data.granuloREQUASUD$Silt,
                                      y = data.granuloLUCAS$Silt)["p.value"]),
                   
                   as.numeric(ks.test(x = data.granuloREQUASUD$Sand,
                                     y = data.granuloLUCAS$Sand)["p.value"])),10))

KS.df %>% 
  kbl(caption = "Kolmogorov-Smirnov (KS) test values and p-value") %>%
  kable_styling(bootstrap_options = "striped", full_width = F)

Comment

The two datasets are different. The clay content shows the highest differences and is higher in the REQUASUD samples.

The reasons of the choice

Soil chemical properties are dependent on the land use. In order to better understand and analyze the chemical property data made available, especially those of REQUASUD, we performed a 5Km buffer around each LUCAS point. These have made it possible to select REAQUSUD points for a finer observation of those according to land use.

Comparison of the OC between requasud and lucas.

  meanCreq_per_lucas <- req.sf.l72 %>%
    st_join(prem_buffer_5k, left = FALSE) %>%
    group_by(POINTID) %>%
    summarise(meanOCreq = mean(OC),
              sdOCreq = sd(OC),
              n = length(POINTID)
              ) %>%
    left_join(select(as_tibble(lucas.sf.l72),c(POINTID,OC,LC0_Desc))) %>%
    mutate(err = abs(1 - (OC / meanOCreq)))

  property <- 'OC'
  my.limit <- +Inf
  my.color <- "LC0_Desc"
  g <- meanCreq_per_lucas %>%
    dplyr::filter(OC < my.limit, meanOCreq < my.limit) %>%
    ggplot(aes(x = OC, y = meanOCreq)) +
    geom_abline(intercept = 0, slope = 1, linewidth = 1, color = "darkgrey") +
    geom_pointrange(aes(
      ymin = (meanOCreq - sdOCreq),
      ymax = (meanOCreq + sdOCreq),
      color = get(my.color)),
      size = 0.1) +
    geom_point(
      aes(color = .data[[my.color]]),
      size = 2
    ) +
    coord_cartesian(
      xlim = c(0, min(my.limit,100)),
      ylim = c(0, min(my.limit,100))
    ) +
    labs(title = "Comparaison of continuous data",
         subtitle = "Buffered Requasud data around Lucas",
         x = "OC (LUCAS)",
         y = "OC (Requasud) mean +/- sd"
         ) +
    scale_color_manual(
      values = c(viridis::viridis(5)[5:1]),
      name = "Land Cover",
      drop = FALSE
    ) +
    theme_classic()
  print(g)

Map of the relative differences between datasets for OC based on land use.

  map.data <- lucas.sf.l72 %>%
    left_join(select(as_tibble(meanCreq_per_lucas),
                     c(POINTID,err,LC0_Desc))
              )

  map.diffrel.lucas <- map.data %>%
    ggplot() +
    geom_sf(data = wallonie, fill = NA) +
    geom_sf(aes(color = LC0_Desc, size = err)) + 
    scale_color_manual(
      values = c(viridis::viridis(5)[5:1]),
      name = "Land Cover",
      drop = FALSE
    ) +
    theme_minimal()
  print(map.diffrel.lucas)

Based on this observation, a filter is applied to the data to take only croplands and grasslands at the level of the two data sets.

To do this we have selected pH(H2O), organic carbon, Phosphorus, Nitrogen and Potassium contents from the LUCAS and REQUASUD datasets after the filter, to compute summary statistics.

#Summary statistic
summary_statLUCAS <- f.summary(data.chemicalLUCAS,param.select =c("pH_H2O","OC","P","N","K"))

summary_statREQUASUD <- f.summary(data.chemicalREQUASUD,param.select =c("pH_H2O","OC","P","N","K"))

Central tendency

We estimate central tendency of LUCAS and REQUASUD dataset

f.wl.stat.table(summary_statLUCAS,summary_statREQUASUD,'central')

Measures of spread

We do it in the same way for these soil properties.

f.wl.stat.table(summary_statLUCAS,summary_statREQUASUD,'spread')

Comparison graphical representations

data.chemicalREQUASUD %>% 
  dplyr::select(pH_H2O,OC, N, P, K) %>%
  pivot_longer(cols = pH_H2O:K,
               names_to = "property",
               values_to = "value") %>%
  mutate(Dataset = "REQUASUD") %>%
  bind_rows(
    data.chemicalLUCAS %>% 
      dplyr::select(pH_H2O,OC, N, P, K)%>% 
      pivot_longer(cols = pH_H2O:K,
                   names_to = "property",
                   values_to = "value") %>%
      mutate(Dataset = "LUCAS")
  ) %>%
  mutate(property = factor(property,levels = c("pH_H2O","OC","P","N","K"))) %>%
  ggplot(aes(x = Dataset ,
             y = value,
             fill = Dataset)
  ) + 
  
  geom_boxplot()+
  labs(title="Plot of soil property by type of soil monitoring network",x="Soil monitoring", y = "Soil porperty")+
  scale_fill_viridis(direction=-1,discrete = TRUE) +
  theme_classic()+
  facet_wrap( ~ property ,
              scales = "free_y"
  )

f.wl.distr(data.chemicalREQUASUD,
           data.chemicalLUCAS,
           wl.param = 'pH_H2O',
           wl.bin = 0.2,
           wl.xlab = "pH_in_H2O [-]"
)

f.wl.distr(data.chemicalREQUASUD,
           data.chemicalLUCAS,
           wl.param = 'OC',
           wl.bin = 10,
           wl.xlab = "Organique carbone [g/Kg]"
)

f.wl.distr(data.chemicalREQUASUD,
           data.chemicalLUCAS,
           wl.param = 'N',
           wl.bin = 0.5,
           wl.xlab = "Nitrogen [g/Kg]"
)

f.wl.distr(data.chemicalREQUASUD,
           data.chemicalLUCAS,
           wl.param = 'K',
           wl.bin = 50,
           wl.xlab = "Potassium [mg/Kg]"
)

f.wl.distr(data.chemicalREQUASUD,
           data.chemicalLUCAS,
           wl.param = 'P',
           wl.bin = 50,
           wl.xlab = "Phosphorus [mg/Kg]"
)

KS.df <- cbind.data.frame(soil.property = c("pH_H2O","OC", "P","N","K"),
                 KS.test = round(c(
                   as.numeric(ks.test(x = data.chemicalREQUASUD$pH_H2O,
                       y=data.chemicalLUCAS$pH_H2O)["statistic"]),
                   
                   as.numeric(ks.test(x = data.chemicalREQUASUD$OC,
                                      y= data.chemicalLUCAS$OC)["statistic"]),
                   
                   as.numeric(ks.test(x = data.chemicalREQUASUD$N,
                                    y=data.chemicalLUCAS$N)["statistic"]),
                                    
                   as.numeric(ks.test(x = data.chemicalREQUASUD$P,
                                    y= data.chemicalLUCAS$P)["statistic"]),
                                    
                   as.numeric(ks.test(x = data.chemicalREQUASUD$K,
                                    y= data.chemicalLUCAS$K)["statistic"])),3),
                 
                 p.value = round(c(
                   
                   as.numeric(ks.test(x = data.chemicalREQUASUD$pH_H2O,
                                    y = data.chemicalLUCAS$pH_H2O)["p.value"]),
                   
                   as.numeric(ks.test(x = data.chemicalREQUASUD$OC,
                                      y = data.chemicalLUCAS$OC)["p.value"]),
                   
                   as.numeric(ks.test(x = data.chemicalREQUASUD$N,
                                    y = data.chemicalLUCAS$N)["p.value"]),
                                    
                   as.numeric(ks.test(x = data.chemicalREQUASUD$P,
                                    y = data.chemicalLUCAS$P)["p.value"]),
                                      
                   as.numeric(ks.test(x = data.chemicalREQUASUD$K,
                                    y = data.chemicalLUCAS$K)["p.value"])),10))

KS.df %>% 
  kbl(caption = "Kolmogorov-Smirnov (KS) test values and p-value") %>%
  kable_styling(bootstrap_options = "striped", full_width = F)

Spatial distribution of clay values

f.wl.property (wallonie,
               data.granuloREQUASUD,
               data.granuloLUCAS,
               property = 'Clay',
               "Long_wgs84",
               "Lat_wgs84",
               coord.crs = 4326,
               title = 'Clay [%]'
) 

Spatial distribution of silt values

f.wl.property (wallonie,
               data.granuloREQUASUD,
               data.granuloLUCAS,
               property = 'Silt',
               "Long_wgs84",
               "Lat_wgs84",
               coord.crs = 4326,
               title = 'Silt [%]'
) 

Spatial distribution of sand values

f.wl.property (wallonie,
               data.granuloREQUASUD,
               data.granuloLUCAS,
               property = 'Sand',
               "Long_wgs84",
               "Lat_wgs84",
               coord.crs = 4326,
               title = 'Sand [%]'
)  

Spatial distribution of pH values

f.wl.property (wallonie,
               data.chemicalREQUASUD,
               data.chemicalLUCAS,
               property = 'pH_H2O',
               "Long_wgs84",
               "Lat_wgs84",
               coord.crs = 4326,
               title = 'pH in water [-]'
)  

Spatial distribution of OC values

f.wl.property (wallonie,
               data.chemicalREQUASUD,
               data.chemicalLUCAS,
               property = 'OC',
               "Long_wgs84",
               "Lat_wgs84",
               coord.crs = 4326,
               title = 'Organique\ncarbone [g/Kg]'
)

Spatial distribution of P values

f.wl.property (wallonie,
               data.chemicalREQUASUD,
               data.chemicalLUCAS,
               property = 'P',
               "Long_wgs84",
               "Lat_wgs84",
               coord.crs = 4326,
               title = 'Phophorus [mg/Kg]'
)

Spatial distribution of N values

f.wl.property (wallonie,
               data.chemicalREQUASUD,
               data.chemicalLUCAS,
               property = 'N',
               "Long_wgs84",
               "Lat_wgs84",
               coord.crs = 4326,
               title = 'Nitrogen [g/Kg]'
)

Spatial distribution of K values

f.wl.property (wallonie,
               data.chemicalREQUASUD,
               data.chemicalLUCAS,
               property = 'K',
               "Long_wgs84",
               "Lat_wgs84",
               coord.crs = 4326,
               title = 'Potasium [mg/Kg]'
)

‘./R/wl_personal_functions.R’

cat(readLines('./R/wl_personal_functions.R'), sep = '\n')

## Filter REQUASUD data
f.requasud <- function(f.year = 'all',
                       f.depth = op.depth){
  
  if (f.year == 'all') {
    df <- data.REQUASUD %>% 
      filter(Depth %in% f.depth)
  } else {
    df <- data.REQUASUD %>% 
      filter(year == f.year & Depth %in% f.depth)
    
  }
  
  return(df)
}

get_png <- function(filename) {
    grid::rasterGrob(png::readPNG(filename), interpolate = TRUE,
                     x = unit(1,"npc"),
                     y = unit(-0.15,"npc"),
                     ## width = 0.33,
                     height = 0.12,
                     vjust = 0,
                     hjust = 1
                     )
}
logos.data <- get_png("./img/montage_logos_graph.png")
text.data <- grobTree(textGrob('Data source:',
                               x=unit(0.66, "npc"), y=unit(-0.075,"npc"),
                               hjust=1))

## Graph distribution (histogram)
## ------------------
f.wl.distr <- function(wl.data.req,
                       wl.data.lucas,
                       wl.param,
                       wl.bin,
                       wl.xlab,
                       wl.title = 'Histogram comparing LUCAS and REQUASUD'){
    g <- ggplot() +
        geom_histogram(data = wl.data.req,
                       aes(y = -(..density..), x = get(wl.param), fill = "Requasud"),
                       col = "black", binwidth = wl.bin) +
        geom_histogram(data = wl.data.lucas,
                       aes(y = (..density..), x = get(wl.param), fill = "Lucas"),
                       col = "black", binwidth = wl.bin) +
        labs(fill ="Dataset", x = wl.xlab, y = " ",title=wl.title) +
        scale_fill_viridis(direction=-1,discrete = TRUE) +
        theme_classic() + 
        theme(legend.position = "right",
              plot.margin = unit(c(0.02, 0.02, 0.12, 0.02), "npc")
              ) +
        coord_flip(clip = "off") +
        annotation_custom(text.data) +
        annotation_custom(logos.data) ## , ymin = 0.05, ymax = 0.1, xmin = -7, xmax = -12)
    return(g)
}

### Statistic
## ------------------
f.summary <- function(data.df,
                      param.select){
  summary_stat <- data.df %>% 
    dplyr::select(param.select) %>%
    describe(ranges = T, quant = c(.25,.75))%>%
    as.data.frame()
  return(summary_stat)
  }

### Table of summary
## ------------------
f.wl.stat.table <- function(data.stat.lucas,
                            data.stat.req,
                            show.stat = 'central' ){
  if (show.stat == 'central') {
    stat.table <- cbind.data.frame(data.stat.lucas %>% 
                              dplyr::select(Mean.LUCAS = mean),
                            data.stat.req %>% 
                              dplyr::select(Mean.REQUASUD = mean),
                            data.stat.lucas %>% 
                              dplyr::select(Median.LUCAS = median),
                            data.stat.req %>% 
                              dplyr::select(Median.REQUASUD = median)) %>%
    round(digits = 1) %>%
    kbl(caption ="Central tendency for LUCAS and REQUASUD datasets" ,row.names = T) %>%
    kable_styling(bootstrap_options = "striped", full_width = F)
  }
  else {
    stat.table <- cbind.data.frame(data.stat.lucas %>% 
                                     dplyr::select(Range.LUCAS = range),
                                   data.stat.req %>% 
                                     dplyr::select(Range.REQUASUD = range),
                                   data.stat.lucas %>% 
                                     dplyr::select(Quantile25.LUCAS = Q0.25),
                                   data.stat.req %>% 
                                     dplyr::select(Quantile25.REQUASUD = Q0.25),
                                   data.stat.lucas %>% 
                                     dplyr::select(Quantile75.LUCAS = Q0.75),
                                   data.stat.req %>% 
                                     dplyr::select(Quantile75.REQUASUD = Q0.75),
                                   data.stat.lucas %>% 
                                     dplyr::select(SD.LUCAS = sd),
                                   data.stat.req %>% 
                                     dplyr::select(SD.REQUASUD = sd),
                                   data.stat.lucas %>% 
                                     dplyr::select(Skewness.LUCAS = skew),
                                   data.stat.req %>% 
                                     dplyr::select(Skewness.REQUASUD = skew),
                                   data.stat.lucas %>% 
                                     dplyr::select(Kurtosis.LUCAS = kurtosis),
                                   data.stat.req %>% 
                                     dplyr::select(Kurtosis.REQUASUD = kurtosis)) %>%
      round(digits = 1) %>%
      kbl(caption = "Spread indicators for LUCAS and REQUASUD datasets",row.names = T) %>%
      kable_styling(bootstrap_options = "striped", full_width = F)
    
  }
  return(stat.table)
}


## Spatial distribution
## ------------------


logos.data.map <- grid::rasterGrob(png::readPNG("./img/montage_logos_graph.png"),
                                   interpolate = TRUE,
                     x = unit(0.01,"npc"),
                     y = unit(0.01,"npc"),
                     ## width = 0.33,
                     height = 0.12,
                     vjust = 0,
                     hjust = 0
                     )
text.data.map <- grobTree(textGrob('Data source:',
                                   x=unit(0.01, "npc"),
                                   y=unit(0.15,"npc"),
                               hjust=0,vjust=0))

f.wl.spatial <- function (region,
                          data.spt.req,
                          data.spt.lucas,
                          x.coord,
                          y.coord,
                          coord.crs,
                          data.used
                          ){
    region = st_transform(region,crs = coord.crs)
    
    plot.spatial.dis <- (ggplot(region) +
                         geom_sf(fill='white') +
                         geom_point(data = data.spt.req,
                                    aes(x = .data[[x.coord]],
                                        y = .data[[y.coord]],
                                        fill = "Requasud"),
                                    size = 2,
                                    shape = 23,alpha = 0.2) +
                         geom_point(data = data.spt.lucas,
                                    aes(x = .data[[x.coord]],
                                        y = .data[[y.coord]],
                                        fill = "Lucas"),
                                    size = 2,
                                    shape = 23) +
                         labs(fill = data.used) +
                         scale_fill_viridis(direction=-1,discrete = TRUE) +
                         annotation_scale(location = "bl",
                                          width_hint = 0.1,
                                          pad_x = unit(0.01, "npc"),
                                          pad_y = unit(0.22, "npc")) +
                         annotation_north_arrow(location = "bl", which_north = "true",
                                                pad_x = unit(0.01, "npc"),
                                                pad_y = unit(0.25, "npc"),
                                                style = north_arrow_fancy_orienteering) +
                         coord_sf(clip = "off") +
                         theme_bw() + 
                         theme(legend.position = "right"
                               ) +
                         annotation_custom(text.data.map) +
                         annotation_custom(logos.data.map)
    )
    return(plot.spatial.dis)
    
}


## Spatial distribution by property
## ------------------

f.wl.property <- function (region,
                           data.spt.req,
                           data.spt.lucas,
                           property,
                           x.coord,
                           y.coord,
                           coord.crs,
                           title
){
  region <- st_transform(region,crs = coord.crs)
  
  req.sf <- st_as_sf(data.spt.req,
                     coords=c(x.coord,y.coord),
                     crs= coord.crs)
  
  lucas.sf <- st_as_sf(data.spt.lucas,
                     coords=c(x.coord,y.coord),
                     crs= coord.crs)
  
  map.req <- tm_shape(region) + 
    tm_fill(col = "white") + 
    tm_polygons () +
    tm_shape(req.sf) + 
    tm_dots(property,palette = viridis(14,direction = -1),size = 0.3,colorNA = NULL) +
    tm_legend(legend.outside = TRUE,
              title = paste0("REQUASUD / ", title))
  
  
  map.lucas <- tm_shape(region) + 
    tm_fill(col = "white") + 
    tm_polygons () +
    tm_shape(lucas.sf) + 
    tm_dots(property,palette = viridis(14,direction = -1),size = 0.3,colorNA = NULL) +
    tm_legend(legend.outside = TRUE,
              title = paste0("LUCAS / ", title))
  
  maplucasVSrequausud <- tmap_arrange(map.lucas , map.req)
  
  return(maplucasVSrequausud)
}

‘./R/wl_personal_functions_fred.R’

cat(readLines('./R/wl_personal_functions_fred.R'), sep = '\n')

f.map.raster.fact <- function(raster, variable, title, unit) {
    map <- ggplot() +
        geom_spatraster(data = raster, aes(fill = .data[[variable]])) +
        geom_sf(data = wallonie, fill = NA,
                colour = "black", linewidth = 1) +
        scale_fill_viridis_d(na.value = NA) +
        labs(fill = unit) +
        ggtitle(title) +
        theme_bw()
    return(map)
}

f.map.raster <- function(raster, variable, title, unit) {
    map <- ggplot() +
        geom_spatraster(data = raster, aes(fill = .data[[variable]])) +
        geom_sf(data = wallonie, fill = NA,
                colour = "black", linewidth = 1) +
        scale_fill_viridis_c(na.value = NA) +
        labs(fill = unit) +
        ggtitle(title) +
        theme_bw()
    return(map)
}

  f.soil.awc.ptf <- function(data.granulo, text.tol = 0.01) {
    ## Estimation of available water capacity
    df.granulo.AWC <- data.granulo %>%
      mutate(sum.granulo = ((Clay+Silt+Sand)/100)) %>% 
      ## CHECK IF NEEDED TODO !!!!
      ## mutate(Clay = Clay/10) %>%
      ## mutate(Silt = Silt/10) %>%
      ## mutate(Sand = Sand/10) %>%
      ## Pedotransfer function from Roman Dobarco et al., 2019
      mutate(FC = 0.270 + ((2.45*10^-3)*Clay) + ( (-1.35*10^-3)*Sand) ) %>% 
      mutate(PWP = 0.08 + ((4.01*10^-3)*Clay) + ( (-2.93*10^-4)*Sand)) %>%
      ## simplified estimation of available water capacity
      mutate(AWC = (FC - PWP)*100) %>%
      filter(!is.na(sum.granulo),
             !(sum.granulo < (1-text.tol)),
             !(sum.granulo > (1+text.tol))
             )

    return(df.granulo.AWC)
  }

  f.soil.tritext <- function(data.granulo,
                             colored.by = "AWC"
                             ) {
      TT <- (ggtern(
          data = data.granulo, ggtern::aes(
                                           x = Sand,
                                           y = Clay,
                                           z = Silt
                                       )) +
          geom_point(
              aes(color = .data[["AWC"]]), ## no variable available BUG
              size=2, alpha=0.5
          ) +
          ## scale_color_continuous(
          ##   low="tomato",
          ##   high="dodgerblue"
          ## ) +
          scale_color_viridis_c(option = "plasma", trans = "sqrt") +
          labs(
              title="Soil Texture Triangle",
              subtitle="Difference between field capacity and wilting point",
              color="Available water capacity (vol %)"
          ) +
          ggtern::theme_bw()
      )
  }