R/mri.R

In mapsRinteractive: Local Adaptation and Evaluation of Raster Maps

#' @title mri
#'
#' @description Local adaptation and evaluation of maps of continuous variables
#' in raster format by use of point location data.
#'
#' @param x SpatRaster. Required. Must be have a defined Cartesian coordinate 
#' system. Data must be continuous. If more than one layer, the first layer 
#' will be used.
#'
#' @param y SpatVector of polygons. Optional. Delineates the area within
#' which the raster layer shall be locally adapted and evaluated. If not provided,
#' the analyses will be performed within the intersect of the raster and the
#' sampled area. Must be have a defined Cartesian coordinate system (same as x). 
#'
#' @param z SpatVector of points Required. Must have at least one column
#' with numerical data and these data must be of the same entity and unit as x
#' (specify this column by argument: field). Must be have a defined Cartesian 
#' coordinate system (same as 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 (unit of the 
#' coordinate reference system) 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 standardized
#' 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 standardized 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 standardized 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 percent of the sill. The sill
#' is by default equal to the variance of the data to be kriged (i.e the point
#' observations or the residuals). Allowed values of ng are within the closed
#' range of 0-1.
#'
#' @param check.data Logical value. Default is TRUE. Shall attributes, geometries and projections of
#' the input data (arguments x, y and z) be checked.
#'
#' @return A list with:
#'
#' 1) 'maps'. A raster stack of 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'. SpatVector of the polygon delineating the
#' mapped area.
#'
#' 3) 'pts'. SpatVector of 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 raster maps with continuous variables. A SpatRaster and a SpatVector of point
#' data (same variable and unit as the raster) are required. A SpatVector of polygons
#' can optionally be used to delineate the area for local adaptation and evaluation.
#' 
#' It is a requirement that all spatial objects (x, y and z) have the same projection.
#' The analyses require a Cartesian coordinate reference system. 
#'
#' Four maps are (created and) evaluated: 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 SpatVector of polygons (argument: y) and the
#' buffered point locations. The buffer width is 1.5*(next largest distance) between
#' one point and its nearest neighbor).
#'
#' The mapsRInteractive algorithmns have been
#' described ad by Piikki et al.(2017) and Nijbroek et al. (2018), where more details
#' can be found .
#'
#' 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 'terra'
#' @import 'gstat'
#' @importFrom 'stats' 'complete.cases'
#' @importFrom 'stats' 'var'
#' @importFrom 'utils' 'read.csv'
#' @importFrom 'utils' 'write.table'
#' @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
#' #load package
#' require(terra)
#' 
#' #create a synthetic example raster dataset
#' rr1<-rast(nrow=10, ncol=10, 
#'   vals= sample(1:4, 100, replace=TRUE), 
#'   crs=crs("EPSG:3857")
#'   )
#' rr2<-disagg(rr1, 4, 'bilinear')
#'
#' #create an example SpatVector of points
#' p<-spatSample(x=rr1, size=30, values=TRUE, as.points=TRUE)
#'
#' #do local evaluation and adaptation of the raster data based on the point data
#' m<-mri(x = rr2, z = p, field ="lyr.1")
#'
#' ##check evaluation measures
#' print(m$evaluation)
#' plot(m$maps)

#' @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.data=TRUE){

  #check input data
  if(check.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]]; names(ordkrig)<-'ordkrig'; x<-c(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]]; names(reskrig)<-'reskrig'; x<-c(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]]; names(regkrig)<-'regkrig'; x<-c(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))
  }