R/mri.R

In soilmapper/mapsRinteractive: Local Adaptation and Evaluation of Digital Soil Maps

#' @title mri
#'
#' @description Local adaptation and evaluation of digital soil maps in raster
#' format by use of point location soil property data.

#' @author Kristin Piikki & Mats Söderström ,  \email{kristin.piikki@@slu.se}

#' @param x Raster dataset. Required. Must be have a defined coordinate system.
#' If the coordinate system is not cartesian, the data will be projected onto
#' the Web Mercator (epsg: 3857) coordinate system before any analyses/tests.
#'
#' @param y SpatialPolygonsDataFrame. Optional. Delineates the area within
#' which the raster layer shall be locally adapted and evaluted. Must be projected. If not
#' projected onto the same coordinate system as x, it will be reprojected to the
#' coordinate system of x. If not provided, the analyses willbe performed within
#' the intersect of the raster and the sampled area."
#'
#' @param z SpatialPointsDataFrame. Required. Must have at least one column
#' with numerical data and these data must be of the same entity unit as the raster
#' (specify this column by argument: field). Must be projected. If not
#' projected to the same coordinate system as x, it will be reprojected to the
#' coordinate system of x.
#'
#' @param field Character value. Required. Name of the column in
#' y with the data that shall be used to locally adapt and evaluate the raster
#'
#' @param edge Numeric value. Optional. Specifies the width (m) of a
#' buffer zone inside the edge of the polygon that is excluded from the analyses.
#' Allowed values are within the closed range of 0-10000.
#'
#' @param filter Positive integer. Optional. No of cells in the side of a square
#' window for mean filtering of x. Filtering is done before any resampling
#' (see argument: resolution). Allowed values are within the closed range of 1-20.
#'
#' @param resolution Positive numeric value. Optional. The resolution (m) to
#' which the imported raster shall be resampled before the adaptation. Allowed
#' values are within the closed range of 0.1-10000. In addition, a resolution that
#' means more than 1E+8 raster cells is not allowed.
#'
#' @param md Character value. Optional. Variogram model type for the stadardized
#' variograms used for ordinary kriging interpolation of observed data or residuals.
#' Variograms are generated by gstat::vgm. Default is "Sph" (spherical model).
#'
#' @param rg Numeric value. Optional. Range of the stadardized variograms used
#' for ordinary kriging interpolation of observed data or residuals. Variograms
#' are generated by gstat::vgm. If no rg is specified it will be set to half of
#' the square root of the mapping area: y (possibly shrinked by edge).
#'
#' @param ng Numeric value. Optional. Nugget of the stadardized variograms
#' used for ordinary kriging interpolation of observed data or residuals.
#' Variograms are  generated by gstat::vgm. The nugget is expressed as a fraction
#' of the sill. A ng = 0.1 means that the nugget is 10 % of the sill. The sill
#' is by default equal to the variance of the data to be krigied (i.e the point
#' observations or the residuals). Allowed values of ng are within the closed
#' range of 0-1.
#'
#' @return A list with:
#'
#' 1) 'maps'. A raster stack with the original raster map ('map'), the map,
#' created by ordinary kriging of observed data ('ordkrig'), by residual kriging
#' ('reskrig') and by regression kriging ('regkrig').
#'
#' 2) 'area'. SpatialPolygonsDataFrame with the polygon delineating the
#' mapped area.
#'
#' 3) 'pts'. SpatialPointsDataFrame with the point locations used for mapping,
#' i.e points falling within the mapped area, excluding points with NA values in
#' the observed values or the values extacted from the original map. The column names mean:
#'     obs = observed values.
#'     map = original map values.
#'     ordkrig_cv = values from the leave-one-out cross validation of the ordinary kriging.
#'     res = residuals (map - obs)
#'     reskrig_cv = values from the leave-one-out cross validation of the residual kriging.
#'     regpred = predicted values from the linear regression (obs = a*map + b)
#'     regres = residuals (regpred - obs)
#'     regkrig_cv = values from the leave-one-out cross validation of the regression kriging.
#'
#' 4) 'evaluation'. a data.frame with evaluation statistics for the original map
#' and the leave-one-out cross-validation of the other mapping methods.
#'
#' 5) 'feedback' a character vector with logged feedback on inputted and used data.

#' @details The mri function is intended for local adaptation and evaluation
#' of large extent digital soil maps. A raster map and point location soil
#' property are required. A SpatialPolygonsDataFrame can optionally be used to
#' delineate the area for local adaptation and evalutaion.
#'
#' All spatial objects must have defined coordinate systems. If the defined
#' coordinate systems are not the same, the point location data and the polygon
#' data will be transformed to the coordinate system of the raster. If the defined
#' coordinate system of the raster is not a cartesian coordinate system, all
#' spatial datasets will be projected onto the Web Mercator coordinate system
#' (epsg: 3857).
#'
#' The four maps are (created and) evaluated are: the original raster map, a map
#' created solely based on the soil samples data (ordinary kriging using a standardized
#' variogram), two maps based on a combination of the raster data and the point
#' observations (regression kriging and residual kriging, both using standardized
#' variograms).
#'
#' The maps are evaluated by leave-one-out cross validation and a number of
#' evaluation measures are computed: the Nash-Sutcliffe modelling efficiency (E),
#' the mean absolute error (MAE; Janssen & Heuberger, 1995), the coefficient of
#' determination of a linear regression between predicted and measured values (r2).
#'
#' The mapped area is the intersection between the original raster
#' map (argument: x), any provided SpatailPolygonsdataframe (argument: y) and the
#' buffered point locations. The buffer width is 1.5*(next largest distance) between
#' one point and its nearest neighbour).
#'
#' The mapsRInteractive algorithmns have been
#' described ad by Piikki et al.(2017) and Nijbroek et al. (2018) (before it was
#' made available as an R package). More details can be found in
#' these publications. It is also implemented in the open Swedish web application
#' for precision agriculture markdata.se and the Sub-Saharan Africa Soil Data Manager.
#'
#' On error: check that required data are provided (arguments x, y, z and field),
#' check that all spatal datasets (arguments x, y, z) are projected, check that
#' they do overlap and check that the arguments edge, filter and resolution have
#' appropriate values.
#'
#' @import 'raster'
#' @import 'gstat'
#' @import 'rgdal'
#' @import 'sp'
#' @import 'rgeos'
#' @importFrom 'stats' 'complete.cases'
#' @importFrom 'stats' 'var'
#' @importFrom 'utils' 'read.csv'
#' @importFrom 'utils' 'write.table'
#' @importFrom 'grDevices' 'is.raster'
#' @importFrom 'stats' 'dist'

#' @references
#' Nijbroek, R., Piikki, K., Söderström, M., Kempen, B., Turner,
#' K. G., Hengari, S., & Mutua, J. (2018). Soil Organic Carbon Baselines for
#' Land Degradation Neutrality: Map Accuracy and Cost Tradeoffs with Respect to
#' Complexity in Otjozondjupa, Namibia. Sustainability, 10(5), 1610.
#' doi:10.3390/su10051610
#'
#' Piikki, K.,Söderström, M., Stadig, H. 2017. Local adaptation of a national digital
#' soil map for use in precision agriculture. Adv. Anim. Biosci. 8, 430–432.
#'
#' Janssen, P.H.M.; Heuberger, P.S.C.1995. Calibration of process-oriented models.
#' Ecol. Model., 831, 55–66.
#'
#' Nash, J.E.; Sutcliffe, J.V. River flow forecasting through conceptual models
#' part I—A discussion of principles. J. Hydrol. 1970, 103, 282–290.

#' @examples
#' ##prepare raster dataset (the soil map to be adapted)
#' data('CLAYr')
#' CLAYr<-data.frame(CLAYr$POINT_X, CLAYr$POINT_Y, CLAYr$clay_percent) #rearrange columns
#' require(raster) #load required package
#' CLAYr<-rasterFromXYZ(CLAYr) #convert to raster
#' prj<-'+proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs' #projection
#' crs(CLAYr)<-crs(prj) #define projection
#' names(CLAYr)<-'clay_percent' #rename (not necessary)
#' ##prepare example point location data
#' data('CLAYs')
#' CLAYs<-data.frame(CLAYs) #convert to data.frame
#' coordinates(CLAYs)<-~ POINT_X + POINT_Y #convert to SpatialPointsDataFrame
#' crs(CLAYs)<-crs(CLAYr) #define projection
#' ##run local adaptation and evaluation
#' mri.out<-mri(x = CLAYr, z = CLAYs, field ='clay_percent')
#' ##check evaluation measures
#' print(mri.out$evaluation)

#' @export
mri<-function( x = NULL, y = NULL, z = NULL, field = NULL, edge = 0, filter = 1,
               resolution = NULL, md = 'Sph', rg = NULL, ng = 0.1){

  #check input data
  a<-check(x=x, y=y, z=z, field=field, edge = edge, filter=filter, resolution=resolution)
  x<-a[[1]]; y<-a[[2]]; z<-a[[3]]; feedback<-a[[4]]

  #interpolate observed data by ordinary kriging using a standardized variogram
  a<-ordkrige(x=x$map, y=y, z=z, field='obs', rg=rg, check_data=F)
  ordkrig<-a[[1]]; x<-stack(x, ordkrig); names(x)[nlayers(x)]<-'ordkrig'
  z<-a[[2]]
  t<-'Ordinary kriging ready.'
  feedback<-fback(t)

  #adapt map by residual kriging using a standardized variogram
  a<-reskrige(x=x$map, y=y, z=z, field='obs', rg=rg, check_data=F)
  reskrig<-a[[1]]; x<-stack(x, reskrig); names(x)[nlayers(x)]<-'reskrig'
  z<-a[[2]]
  t<-'Residual kriging ready.'
  feedback<-fback(t)

  #adapt map by regression kriging using a standardized variogram
  a<-regkrige(x=x$map, y=y, z=z, field='obs', rg=rg, check_data=F)
  regkrig<-a[[1]]; x<-stack(x, regkrig); names(x)[nlayers(x)]<-'regkrig'
  z<-a[[2]]
  t<-'Regression kriging ready.'
  feedback<-fback(t)

  #compute evaluation statistics
  evaluation<-evaluate(df=as.data.frame(z), observed='obs', predicted=c('map', 'ordkrig_cv', 'reskrig_cv', 'regkrig_cv'))
  t<-'Evaluation ready.'
  feedback<-fback(t)
  #return
  return(list(maps=x, area=y, pts=z, evaluation=evaluation, feedback=feedback))
  }