Reading ISRIC Phosphorus Retention Data

1 Loading Libraries

library(devtools)
library(here)
load_all(pkg = here())
library(dplyr)
library(tidyr)
library(Hmisc)
library(MASS)
library(ggplot2)
library(ade4)

2 Reading ISRIC 2011 data

2.1 Starting with the base raster file

First we load the raster.

data_dir <- system.file("data", package = "soilP", mustWork = TRUE)
extdata_dir <- system.file("extdata", package = "soilP", mustWork = TRUE)
test_tif <- file.path(extdata_dir,
            "Global_distribution_of_soil_phosphorus_retention_potential",
            "Global_distribution_of_soil_phosphorus_retention_potential.tif")
test_raster <- raster::raster(test_tif)

This raster has a sigle layer and its value corresponds to the soil map unit identifier (MUID) for the Soil Map of the World, FAO74. The same identifier can be found in more recent versions of the FAO soil map and in the harmonized soil database.

2.2 Plotting FAO Soil Map Unit Raster

op <- par()
par(mar = c(0, 0, 0, 0), oma = c(0, 0, 0, 0))
plot(test_raster, axes = FALSE, box = FALSE,
     frame.plot = FALSE,
     legend.width = 2,
     legend.shrink = 0.95,
     axis.args = list(cex.axis = 2),
     legend.args = list(text = "", cex = 1))
par(op)

The MUIDs have the lowest values in Africa and then ascend through Asia, the Americas, Oceania and attain the highest values in Europe. Contrary to what the file name seems to imply these ids do not directly represent the phosphorus retention potential. Batjes used the ids and a series of SQL queries to assign phosphorus retention potential classes to each map unit, which then are easily projected over each pixel in the map.

The phosphorus retention potential class was assigned to each of the FAO74 soil units (Low, Moderate, High and very High). Each map unit is composed of up to 8 soil units in diferent percentages, depending on the composition of the map unit and each soil unit compoment phosphorus retention potential class a class was assigned to each map unit. In total 260 different MU phosphorus retention potential clases were assigned to the 4923 map units. These were further grouped into 16 main phosphorus retention potential classes.

These main classes can be extracted from the map unit raster, and the tables in the mdb files. These tables function as raster attribute tables.

2.3 Extracting Raster Attribute Table from Microsoft Access database (mdb file).

For this I need Hmisc and the ODBC drivers. Mac installation is from brew.

The MS Access database has the tables for mapping the soil unit identifiers to soil properties, i.e. phosphorus retention as soil unit class percentage, or map unit class.

mdb_file <- file.path(extdata_dir,
                  "Global_distribution_of_soil_phosphorus_retention_potential",
                  "ISRIC_Phosphorus_Retention_Potential.mdb")
RATs <- read_ISRIC_RAT(mdb_file)
FAO74    <- RATs$FAO74
mapunit  <- RATs$mapunit
soilunit <- RATs$soilunit

3 Rebuilding the arcgis render appeareance

3.1 Applying Raster Attributes for obtainning P Retention Potential MainClass

We will use a manually curated raster attribute table stored in th soilclass dataframe.

See the manual for details.

?soilclass

Using development documentation for soilclass

First we need to merge soilclass to the Soil Map Unit id number.

data(soilclass)
mu_soilclass <- mapunit %>% 
  dplyr::left_join(soilclass, by=c("mainclass" ="main")) %>%
  dplyr::select(MUID,mainclass,ascending)

Column `mainclass`/`main` joining factors with different levels, coercing to character vector

ISRIC2011 <- read_P_ISRIC(tif = test_tif, soilclass = mu_soilclass,
                     is = "MUID",
                     becomes = "ascending",
                     filename = file.path(extdata_dir,"tif","mu_soilclass.tif"))
# Save ISRIC!!!!!

3.2 Plotting Soil Phosphorus Retention Potential Raster

op <- par()
par(mar = c(0, 0, 0, 0), oma = c(0, 0, 0, 0))
plot(ISRIC2011, axes = FALSE, box = FALSE,
     frame.plot = FALSE,
     legend.width = 2,
     legend.shrink = 0.95,
     axis.args = list(breaks=1:14, at = 0:15, cex.axis = 2),
     legend.args = list(text = "", cex = 1))
par(op)

Now the numbers do represent the phosphorus retention potential!!!

However, the raster plot legend above assumes a continuous scale from 0 to 15, while the data is explicitly a categorical variable (although derived from continuous percentages, see Batjes 2011).

To remediate this I manually assigned integers inascending order to the ascending column in the soilclass dataframe as follows: For the Low map unit main class, higher percentage of Low soil unit P retention corresponded to lower integers. For the Moderate, High, and very High map unit main classes, higher integers were assigned to higher percentage of corresponding soil unit P retention. This means the higher the Low soil unit P retention percentage the lower the integer, and complentarily the higher the Moderate, High, and Very High soil unit P retention percentage the higher the integer. This assignation is naive and consequently both, intuitive and ad hoc.

3.3 Adding Raster Attribute Table to ISRIC2011

levels<- S3 method is internal to raster package and could not be exported to the soilP package, which prevented me from making a working package function out of the following snippet:

# Raster Attribute Table
rat <- levels(ISRIC2011)[[1]] %>% 
       dplyr::inner_join(soilclass, by = c("ID" = "ascending"))
levels(ISRIC2011) <- rat # This assignment method is not exported from raster
save(ISRIC2011, file = file.path(data_dir, "ISRIC2011.RData"))

3.4 Using rasterVis to appropiately label categories and show legend

Remember that we assigned to the to the phosphorus retention potential main classes an ad hoc order so we can interpret a direction of ascending retention potential (see above) in the plot legend.

In order to show the right labels in this ascending scale we use the levelplot function from RasterVis that uses the raster attribute table to transform the numerical value of the raster to discrete categories.

rasterVis::levelplot(ISRIC2011, att = "main",
          col.regions = soilclass$color_hex[1:16],
          maxpixels = ncell(ISRIC2011),  scales = list(draw = FALSE),
          xlab = NULL, ylab = NULL,
          main = "Generalized Phosphorus Retention Potential Map")

4 Multivariate analysis of Soil retention Potential

We should establish a data based order of map unit soil phosphorus retention classes instead of postulating an ad hoc order. This can be done through selecting the first discriminant dimension from a Discriminant Analysis of the percentage of soil unit phosphorus retention classes per mapping unit. Furthermore, we can use that discriminant dimention as a continuous variable instead of the discrete classes for downstream analyses

4.1 Map Unit Soil Composition Matrix

Each soil unit percentage can be used as a map unit descriptor variable for multivariate analysis of the phosphorus retention potential per map unit. This results in sparse matrix with essentially orthogonal columns.

mu_comp <- get_mu_composition(mapunit,FAO74)

4.2 Map unit - Phosphorus Retention Join table

mu_ret <- mapunit %>%
  dplyr::left_join(soilunit, by = c("SOIL1" = "key")) %>%
  dplyr::select(MUID,mainclass,Lo,Mo,Hi,VH,MISC,pH,SAND,SILT,CLAY,CECLAY)

Column `SOIL1`/`key` joining factors with different levels, coercing to character vector

is_soil <- function(x){
  !(x %in% c("GL1","RK1","RK2","WR1"))
}
ret_comp <- subset(cbind(mu_ret,mu_comp),
                    is_soil(mu_ret$mainclass))
soilclass <- soilclass[order(soilclass$ascending),]
ascending_order <- as.character(soilclass$main)[5:16]
ret_comp$mainclass <- factor(ret_comp$mainclass,
                             levels = ascending_order)

5 Discriminant analysis using ade4

5.1 Soil Retention Potential Space is essentially a tetrahedron

with Moderate retention at its center.

Only three degrees of freedom, Lo-VH, Mo, Mo-Hi,

# create dudi object for discriminant analysis
pca <- dudi.pca(ret_comp[,-(1:2)],scannf = FALSE, nf= 4)
dis <- discrimin(pca,factor(ret_comp$mainclass),  scannf = FALSE, nf = 3)
palette <- soilclass$color_hex[5:16]
color <- palette[ret_comp$mainclass]
par(mfrow = c(2, 2))
s.class(dis$li, xax = 2, yax = 1, fac = ret_comp$mainclass, col = palette)
s.class(dis$li, xax = 3, yax = 1, fac = ret_comp$mainclass, col = palette)

s.class(dis$li, xax = 3, yax = 2, fac = ret_comp$mainclass, col = palette)
s.corcircle(dis$va)

par(mfrow = c(1, 1))

6 Discriminant analysis using MASS

fit <- lda(mainclass ~ ., data=ret_comp[,-1])

variables are collinear

plda <-  predict(object = fit, # predictions
                 newdata = ret_comp[,-(1:2)])

7 Linearization of Phosphorus retention potential

From categorical variable to continous “linear” scale. Linear combination of map unit soil content.

ret_scale <- ret_comp[,c("MUID","Lo","Mo","Hi","VH","MISC")]
ret_scale$ret <- ret_comp$mainclass
# Weighted scale, ad hoc weights!!!!!!
ret_scale$weighted <- with(ret_scale,
                           {as.numeric((Lo + 2 * Mo + 3 * Hi + 4 * VH) / 400)})
# Linear Discrimant scale, post hoc weights/coefficients
ret_scale$LD2 <-  plda$x[,2]

plot_P_scales(ret_scale, palette = palette,
              scale_x = "LD2",scale_y = "weighted")

scale256 <- function(x){
  as.integer(round(245 * (x - min(x,na.rm = TRUE))/(max(x,na.rm = TRUE) -
                                                      min(x,na.rm = TRUE))) +5)}
P_retention <- mu_soilclass %>% 
  dplyr::left_join(ret_scale, by = "MUID")
P_retention$weighted <- scale256(P_retention$weighted)
P_retention$weighted[is.na(P_retention$weighted)] <- 0
P_retention$LD2 <- scale256(P_retention$LD2) 
P_retention$LD2[!is_soil(P_retention$mainclass)] <- 0
P_weighted <- read_P_ISRIC(tif = test_tif, soilclass = P_retention,
                     is = "MUID",
                     becomes = "weighted",
                     filename = file.path(extdata_dir,"tif","P_weighted.tif"))
P_LD2 <- read_P_ISRIC(tif = test_tif, soilclass = P_retention,
                     is = "MUID",
                     becomes = "LD2",
                     filename = file.path(extdata_dir,"tif","P_LD2.tif"))

7.1 Weighted P retention Potential

op <- par()
par(mar = c(0, 0, 0, 0), oma = c(0, 0, 0, 0))
plot(P_weighted, axes = FALSE, box = FALSE,
     frame.plot = FALSE,
     legend.width = 2,
     legend.shrink = 0.95,
     axis.args = list(cex.axis = 2),
     legend.args = list(text = "", cex = 1))
par(op)

7.2 Linear Discriminant P retention Potential

op <- par()
par(mar = c(0, 0, 0, 0), oma = c(0, 0, 0, 0))
plot(P_LD2, axes = FALSE, box = FALSE,
     frame.plot = FALSE,
     legend.width = 2,
     legend.shrink = 0.95,
     axis.args = list(cex.axis = 2),
     legend.args = list(text = "", cex = 1))
par(op)

8 Soil Map Phosphorus Retention Potential as Percentages of Map Units

P_retention[!is_soil(P_retention$mainclass),c("Lo","Mo","Hi","VH")] <- 0
P_retention[!is_soil(P_retention$mainclass),"MISC"] <- 100
ISRIC2011_Lo <- read_P_ISRIC(tif = test_tif, soilclass = P_retention,
                     is = "MUID",
                     becomes = "Lo",
                     filename = file.path(extdata_dir,"tif","ISRIC2011_Lo.tif"))
ISRIC2011_Mo <- read_P_ISRIC(tif = test_tif, soilclass = P_retention,
                     is = "MUID",
                     becomes = "Mo",
                     filename = file.path(extdata_dir,"tif","ISRIC2011_Mo.tif"))
ISRIC2011_Hi <- read_P_ISRIC(tif = test_tif, soilclass = P_retention,
                     is = "MUID",
                     becomes = "Hi",
                     filename = file.path(extdata_dir,"tif","ISRIC2011_Hi.tif"))
ISRIC2011_VH <- read_P_ISRIC(tif = test_tif, soilclass = P_retention,
                     is = "MUID",
                     becomes = "Hi",
                     filename = file.path(extdata_dir,"tif","ISRIC2011_VH.tif"))
ISRIC2011_MISC <- read_P_ISRIC(tif = test_tif, soilclass = P_retention,
                     is = "MUID",
                     becomes = "MISC",
                     filename = file.path(extdata_dir,"tif","ISRIC2011_MISC.tif"))

8.1 Lo P retention Potential

op <- par()
par(mar = c(0, 0, 0, 0), oma = c(0, 0, 0, 0))
plot(ISRIC2011_Lo, axes = FALSE, box = FALSE,
     frame.plot = FALSE,
     legend.width = 2,
     legend.shrink = 0.95,
     axis.args = list(cex.axis = 2),
     legend.args = list(text = "", cex = 1))
par(op)

8.2 Mo P retention Potential

op <- par()
par(mar = c(0, 0, 0, 0), oma = c(0, 0, 0, 0))
plot(ISRIC2011_Mo, axes = FALSE, box = FALSE,
     frame.plot = FALSE,
     legend.width = 2,
     legend.shrink = 0.95,
     axis.args = list(cex.axis = 2),
     legend.args = list(text = "", cex = 1))
par(op)

8.3 Hi P retention Potential

op <- par()
par(mar = c(0, 0, 0, 0), oma = c(0, 0, 0, 0))
plot(ISRIC2011_Hi, axes = FALSE, box = FALSE,
     frame.plot = FALSE,
     legend.width = 2,
     legend.shrink = 0.95,
     axis.args = list(cex.axis = 2),
     legend.args = list(text = "", cex = 1))
par(op)

8.4 VH P retention Potential

op <- par()
par(mar = c(0, 0, 0, 0), oma = c(0, 0, 0, 0))
plot(ISRIC2011_VH, axes = FALSE, box = FALSE,
     frame.plot = FALSE,
     legend.width = 2,
     legend.shrink = 0.95,
     axis.args = list(cex.axis = 2),
     legend.args = list(text = "", cex = 1))
par(op)

8.5 MISC P retention Potential

op <- par()
par(mar = c(0, 0, 0, 0), oma = c(0, 0, 0, 0))
plot(ISRIC2011_MISC, axes = FALSE, box = FALSE,
     frame.plot = FALSE,
     legend.width = 2,
     legend.shrink = 0.95,
     axis.args = list(cex.axis = 2),
     legend.args = list(text = "", cex = 1))
par(op)