R/interp_inputpts.R

In ranghetti/guinterp: A Shiny GUI for Spatial Interpolation

#' @title Interpolate points of a single polygon
#' @description TODO
#' @param indata_sf points of a single field (with optional `"idfield"` attribute)
#' @param field_poly polygon of the field
#' @param vgm.fit variogram
#' @param globgrid grid generated by `make_interp_grid()`
#' @param method one between `"krige"` (default) and `"idw"`
#' @param n_cores 1 for single core (default), more for multicore, NA for multicore autodetection
#' @param focal_types smoothing filters to be applied (vector of types of `focalWeight()`)
#' @param focal_d list of filter parameters (see argument d of `focalWeight()`)
#' @param border border (in CRS unit) to be cut (default: 0)
#' @param outdir directory where rasters are stores (default: temporary directory)
#' @param outname name of the output tif
#' @param samplesize maximum size of the sample of the original data to work with (default: 10000; if NA: all the points)
#' @param nmax argument of `krige()` and `idw()`
#' @param maxdist argument of `krige()` and `idw()`
#' @importFrom raster focal focalWeight raster
#' @importFrom methods as
#' @importFrom stars st_dimensions write_stars
#' @importFrom magrittr "%>%"
#' @importFrom sf st_crop st_union st_buffer st_crs st_transform st_set_crs
#' @author Luigi Ranghetti, phD (2018) \email{ranghetti.l@@irea.cnr.it}
#' @note License: GPL 3.0

interp_inputpts <- function(
  indata_sf, # points of a single field (with optional "idfield" attribute)
  field_poly, # polygon of the field
  vgm.fit, # variogram
  globgrid, # grid generated by make_interp_grid()
  method = "krige", # one between "krige" and "idw"
  n_cores = 1, # 1 for singlecore, more for multicore, NA for multicore autodetection
  focal_types = NULL, # smoothing filters to be applied (vector of types of focalWeight())
  focal_d = NULL, # list of filter parameters (of)see argument d of focalWeight())
  border = 0, # border (in crs unit) to be cut
  outdir = tempdir(), # directory where rasters are stored
  outname = NA, # name of the output tif
  samplesize = 1E4, # maximum size of the sample of the original data to work with (default: 10000; if NA: all the points)
  nmax=1E3, # argument of krige() and idw()
  maxdist=Inf # argument of krige() and idw()
) {

  # Checks

  # check that idfield is unique
  if ("idfield" %in% names(indata_sf)) {
    sel_idfield <- unique(indata_sf$idfield)
    if (length(sel_idfield)>1) {
      stop("Error: idfield attribute must contain only one value.")
    }
  } else {
    sel_idfield <- NA
  }
  # for now, no checks are performed on "idfield" attribute in field_poly
  # (the function assumes that field_poly contains only the geometry of sel_idfield)

  # check field_poly crs
  if (st_crs(field_poly) != st_crs(globgrid)) {
    field_poly <- st_transform(field_poly, st_crs(globgrid))
  }
  # check indata_sf crs
  if (st_crs(indata_sf) != st_crs(globgrid)) {
    indata_sf <- st_transform(indata_sf, st_crs(globgrid))
  }

  if (is.na(samplesize)) {
    samplesize <- nrow(indata_sf)
  }
  if (samplesize < nrow(indata_sf)) {
    indata_sf_sub <- indata_sf[indata_sf$sid <= samplesize,]
  } else {
    indata_sf_sub <- indata_sf
  }

  # crop grid on field_poly extent
  if (!is.null(focal_d)) {
    d_buffer <- max(c(sapply(focal_d, max),unlist(focal_d[focal_types=="Gauss"])*3))
  } else {
    d_buffer <- 3
  }

  # cut on borders ("buffer" value plus half of raster resolution, to cut pixels of border)
  # and create buffered grid (for focals)
  field_poly_buf <- st_buffer(st_union(field_poly), st_dimensions(globgrid)[[1]]$delta*d_buffer)
  fieldgrid_buf <- st_crop(globgrid, field_poly_buf)
  fieldgrid_buf <- fieldgrid_buf[1,]
  fieldgrid <- st_crop(
    fieldgrid_buf,
    st_buffer(field_poly, -(border+st_dimensions(globgrid)[[1]]$delta/2)),
    crop = FALSE
  )

  # interpolate
  interp_raster <- krige_par(selvar ~ 1, indata_sf_sub, fieldgrid, vgm.fit, method=method, n_cores=n_cores, nmax=nmax, maxdist=maxdist)

  # export raster
  if (is.na(outname)) {
    if (is.na(sel_idfield)) {
      outname <- paste0('polygon_',str_pad2(sample(1E4,1),4,"left","0"),'.tif')
    } else {
      outname <- paste0("polygon_",sel_idfield,".tif")
    }
  }

  dir.create(file.path(outdir,"interp"), recursive=FALSE, showWarnings=FALSE)

  write_stars(
    interp_raster["var1.pred",],
    file.path(outdir,"interp",outname),
    NA_value = -32768,
    options = "COMPRESS=DEFLATE"
  )
  message(file.path(outdir,outname))

  # extend interp grid, using average value as NA
  raster_crop_pred <- interp_raster["var1.pred",]
  raster_crop_pred[[1]][is.na(interp_raster$var1.pred)] <- mean(interp_raster$var1.pred, na.rm=TRUE)

  # apply filters
  raster_crop_pred_r <- as(raster_crop_pred, "Raster")
  for (i in seq_along(focal_types)) {
    sel_focal_type <- focal_types[i]
    sel_focal_d <- focal_d[[i]]
    sel_focal_raster_pred <- focal(
      raster_crop_pred_r,
      w=focalWeight(raster_crop_pred_r, sel_focal_d, sel_focal_type),
      na.rm=TRUE
    ) %>% st_as_stars() %>% st_set_crs(st_crs(raster_crop_pred))

    # cut on borders ("buffer" value plus half of raster resolution, to cut pixels of border)
    crop_focal_raster_pred <- st_crop(
      sel_focal_raster_pred,
      st_buffer(field_poly, -(border+st_dimensions(globgrid)[[1]]$delta/2))
    )
    # export raster
    subdir <- paste0("focal_",sel_focal_type,sel_focal_d,"x",sel_focal_d)
    dir.create(file.path(outdir,subdir), recursive=TRUE, showWarnings=FALSE)
    write_stars(
      crop_focal_raster_pred,
      file.path(outdir,subdir,outname),
      NA_value = -32768,
      options = "COMPRESS=DEFLATE",
      overwrite = TRUE
    )
  }
  # } # end of field_id cycle

}