GUInterp from commandline

# ---------- Esempio di utilizzo di GUInterp da riga di comando ---------- #
#                                                                          #
#     Script di esempio per riprodurre a linea di comando una tipica       #
#     interpolazione impostata e avviata normalmente da GUI.               #
#     Alcune opzioni impostabili da GUI non sono state incluse per non     #
#     appesantire il codice.                                               #
#                                                                          #
#     Autore: Luigi Ranghetti                                              #
# ------------------------------------------------------------------------ #

library(guinterp)
library(sf)
library(dplyr)

### Impostazione dei parametri ----

## Vettoriale dei poligoni
borders_path <- "/path/of/id_geom_0045.shp"
select_uid_which <- "attr" # ID univoco delle geometrie ("record", "attr" o "no")
select_uid_attr <- "id_geom" # campo contenente gli ID (se select_uid_which = "attr")

# Punti da interpolare
inputpts_path <- "/path/of/resA_id_geom_0045.shp"
select_inputvar <- "RESAKG" # Variabile da interpolare


## Filtra dati
filter_buttons <- "manual" # Filtraggio dei punti ("manual", "minimal", "no")
filter_rangey <- c(-Inf, Inf) # Valori consentiti (se filtraggio manuale)
filter_zscorey <- c(-Inf, Inf) # Zscore (se filtraggio manuale)
filter_rbiasy <- c(-Inf, Inf) # Rbias (se filtraggio manuale)
filter_rangeq <- c(2, 98) # Intervallo dei quantili (0-100) (se filtraggio manuale)
# non sono stati impostati i filtri per poligono e per selezione manuale dalla mappa

## Formato di output
interp_res <- 5 # Risoluzione (m)
interp_crs <- st_crs(32632) # Sistema di riferimento
# Le opzioni si riferiscono all'impstazione "Definisci risoluzione e proiezione";
# per scegliere da un raster esistente, oltre a scegliere risoluzione e CRS
# del raster bisognerebbe impostare anche "grid_offset" -
# lo documenterò all'occorrenza.

## Opzioni di interpolazione
interp_maxptdist <- 15 # Distanza massima dai punti
interp_method <- "krige" # Metodo di interpolazione ("krige" o "idw")
# Definizione del variogramma: l'esempio funziona nel caso di "Automatica";
# altre modalità richiedono la definizione manuale del variogramma con metodi
# e funzioni del pacchetto {gstat}.
focal_onoff <- TRUE # Livella i raster interpolati

## Opzioni di processamento
interp_parallel <- TRUE # Ottimizza l'interpolazione per uso desktop (FALSE) o server/workstation (TRUE)
out_path <- "/path/of/raster_out.tif" # Raster di output
# Opzioni avanzate
v_maxdist <- NA # Distanza massima punti-pixel (auto: NA)
v_nmax <- 100 # Numero di punti per pixel (auto: 500)
interp_samplesize <- NA # Usa un campione di punti per l'interpolazione (NA/Inf per usarli tutti)
interp_samplescheme <- "strat_area" # Schema di campionamento nei poligoni ("random", "strat_npts", "strat_area" o "strat_prop")


### Esecuzione ----

# Carica i file vettoriali
borders_sf0 <- st_read(borders_path) %>%
  st_transform(4326)
inputpts_sf <- st_read(inputpts_path) %>%
  st_transform(4326)

# Converti il vettoriale dei bordi nel formato usato da guinterp
borders_sf <- if (select_uid_which == "attr") {
  borders_sf0 %>%
    dplyr::select(select_uid_attr) %>%
    dplyr::rename(id_geom = select_uid_attr)
} else if (select_uid_which == "record") {
  borders_sf0 %>%
    dplyr::transmute(id_geom = seq_len(nrow(borders_sf0)))
} else if (select_uid_which == "no") {
  borders_sf0 %>%
    dplyr::transmute(id_geom = "all points")
} %>%
  dplyr::group_by(id_geom) %>% dplyr::summarise()

# Converti il vettoriale dei punti nel formato usato da guinterp
inputpts_dt <- read_inputpts(
  inputpts_sf,
  borders = borders_sf,
  varname = select_inputvar
)

# Filtra i punti
if (filter_buttons == "manual") { # applica i filtri manuali
  inputpts_dt <- filter_pts(
    inputpts_dt, "rangey", filter_rangey,
    id_fieldname="id_geom",
    inlayer = borders_sf,
    byfield = TRUE, samplesize = NA
  )
  inputpts_dt <- filter_pts(
    inputpts_dt, "zscorey", filter_zscorey/100,
    id_fieldname="id_geom",
    inlayer = borders_sf,
    byfield = TRUE, samplesize = NA
  )
  inputpts_dt <- filter_pts(
    inputpts_dt, "rbiasy", filter_rbiasy/100,
    id_fieldname="id_geom",
    inlayer = borders_sf,
    byfield = TRUE, samplesize = NA
  )
  inputpts_dt <- filter_pts(
    inputpts_dt, "rangeq", filter_rangeq/100,
    id_fieldname="id_geom",
    inlayer = borders_sf,
    byfield = TRUE, samplesize = NA
  )
  inputpts_dt <- filter_pts(
    inputpts_dt, "rangey", filter_rangey,
    id_fieldname="id_geom",
    inlayer = borders_sf,
    byfield = TRUE, samplesize = NA
  )
} else if (filter_buttons == "minimal") { # applica il filtro minimo (2°-98° percentile)
  inputpts_dt <- filter_pts(
    "rangeq", c(.02,.98),
    inlayer = borders_sf,
    byfield = TRUE, samplesize = NA
  )
}

# Interpola
out_path <- guinterp_process(
  inputpts_dt, borders_sf,
  id_fieldname = "id_geom",
  out_path = out_path,
  samplesize = interp_samplesize,
  samplescheme = interp_samplescheme,
  parallel = interp_parallel,
  interp_method = interp_method,
  smooth = focal_onoff,
  interp_res = interp_res,
  out_crs = interp_crs,
  buffer_radius = interp_maxptdist,
  vgm = "auto",
  v_nmax = v_nmax,
  v_maxdist = v_maxdist
)

# Visualizza output
mapview::mapview(raster(out_path))