In mbask/fertplanspatial: Compute And Spatialize NPK Fertilization Plans

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>")
options(warn = -1)
library(fertplanspatial)
library(ggplot2)
library(raster)
library(sp)

Fertilization plans for nitrogen, phosphorus, and potassium nutrients

Nitrogen, phosphorus and potassium are among most needed nutrient elements by agriculture crops. Their demand depends on a few key factors including future and past crop type, soil texture, soil organic matter amount, content in Calcium, amount of rainfall and many others.

The workflow for computing nutrient needs by crops starts with the analysis of the soil, usually by digging soil samples from the field to be fertilized. Samples are analysed to reveal their physical and chemical features that, along with other general variables, enter the computation of the fertilization components for N, P, and K. Farmers are interested in knowing the average demand among the sampled points. More info on this and on the components of the fertilization plans can be found in the documentation for the fertplan package.

Averaging nutrient demand between soil samples may give a reasonable idea of the crop demand but no clear hint at its spatial variation over the crop area. A few factors may shift nutrient demand over the field area, such as the presence of water tables under the soil surface, or the proximity of a hilly area where nutrients can be moved from, or a previous over- or under-fertilization in localized area, and so on. The ability to spatialize soil demand concentrations from the soil sample points over the whole field may allow farmers to precisely calibrate nutrient spreasing.

We will now be using the soil_spatial builtin dataset. This dataset collects nutrient demands for the same 20 soil sample as the soils dataset in package fertplan and includes samples geographic coordinates as "X", and "Y":

data("soils_spatial")
knitr::kable(soils_spatial)

Spatialisation is carried out by ordinary kriging through the function gstat::krige. Function spatial_nutrient is a wrapper around gstat::vgm and automap::autofitVariogram. Which function is used depends on the model argument, defaulting to "auto" and, therefore, automap::autofitVariogram. This latter function automagically fit a variogram to the data whereas the former function enables to fit a specific variogram, for each specific nutrient for even more preciseness, by passing the appropriate named arguments to spatial_nutrient. Ordinary kriging is performed on a spatial grid built on the soil sampling points bounding box by default. Should a specific grid be needed (ie a larger extent than the one where soils sampled were dug) it should be passed as a sp::SpatialPoints() object. We will simply let automap::autofitVariogram select the optimal variogram and pass it to ordinary kriging on the default bounding box, with a 10 metres spatial resolution:

spatials_l <- spatial_nutrient(soils_spatial, spat_res = 10)

Finally we can plot the spatial distribution of the nutrient fertilization plans:

ggplot(as.data.frame(spatials_l$n), aes(x = X, y = Y)) +
  geom_tile(aes(fill = var1.pred)) +
  geom_point(data = as.data.frame(soils_spatial)) +
  coord_equal() +
  scale_fill_distiller(palette = "Spectral") +
  labs(title = "N fertilization plan", fill = "Concentration\n(kg/ha)") +
  theme_bw()
last_plot() %+% as.data.frame(spatials_l$p) %+% labs(title = "P fertilization plan")
last_plot() %+% as.data.frame(spatials_l$k) %+% labs(title = "K fertilization plan")

A rough interpretation of the plans clearly suggests that there is a need for a light fertilization by phosphorus, that potassium is in great excess and nitrogen is more than enough for the successive crop.

For a greater flexibility the variogram step of the spatialization process can be specifically tailored to use a specific function among [gstat::vgm()] (controlled by argument model set to a specific vgm model) and [automap::autofitVariogram()] (controlled by argument model = auto). Further arguments can be passed to either functions:

spatials_l <- c(
  spatial_nutrient(soils_spatial, model = "auto", spat_res = 10, nutrient = "nitrogen", alpha = seq(0, 359, 15)),
  spatial_nutrient(soils_spatial, model = "Ste",  spat_res = 10, nutrient = "phosphorus"),
  spatial_nutrient(soils_spatial, model = "auto", spat_res = 10, nutrient = "potassium", alpha = seq(0, 359, 15)))
last_plot() %+% as.data.frame(spatials_l$n) %+% labs(title = "N fertilization plan")
last_plot() %+% as.data.frame(spatials_l$p) %+% labs(title = "P fertilization plan")
last_plot() %+% as.data.frame(spatials_l$k) %+% labs(title = "K fertilization plan")

The fertilization plan can be easily converted into a raster stack and saved as GeoTIFF image for further GIS elaboration:

# coerce to SpatialPixelsDataFrame
kriges_spdf_l <- lapply(spatials_l, function(spg) { sp::gridded(spg) <- TRUE; spg })
krige_rs      <- raster::stack(lapply(kriges_spdf_l, raster::raster))
raster::writeRaster(krige_rs, "npk_fert_plans.tif", format = "GTiff")