ggranga

GUInterp from commandline

Oct 27th, 2020 (edited)
1,041
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. # ---------- Esempio di utilizzo di GUInterp da riga di comando ---------- #
  2. #                                                                          #
  3. #     Script di esempio per riprodurre a linea di comando una tipica       #
  4. #     interpolazione impostata e avviata normalmente da GUI.               #
  5. #     Alcune opzioni impostabili da GUI non sono state incluse per non     #
  6. #     appesantire il codice.                                               #
  7. #                                                                          #
  8. #     Autore: Luigi Ranghetti                                              #
  9. # ------------------------------------------------------------------------ #
  10.  
  11. library(guinterp)
  12. library(sf)
  13. library(dplyr)
  14.  
  15. ### Impostazione dei parametri ----
  16.  
  17. ## Vettoriale dei poligoni
  18. borders_path <- "/path/of/id_geom_0045.shp"
  19. select_uid_which <- "attr" # ID univoco delle geometrie ("record", "attr" o "no")
  20. select_uid_attr <- "id_geom" # campo contenente gli ID (se select_uid_which = "attr")
  21.  
  22. # Punti da interpolare
  23. inputpts_path <- "/path/of/resA_id_geom_0045.shp"
  24. select_inputvar <- "RESAKG" # Variabile da interpolare
  25.  
  26.  
  27. ## Filtra dati
  28. filter_buttons <- "manual" # Filtraggio dei punti ("manual", "minimal", "no")
  29. filter_rangey <- c(-Inf, Inf) # Valori consentiti (se filtraggio manuale)
  30. filter_zscorey <- c(-Inf, Inf) # Zscore (se filtraggio manuale)
  31. filter_rbiasy <- c(-Inf, Inf) # Rbias (se filtraggio manuale)
  32. filter_rangeq <- c(2, 98) # Intervallo dei quantili (0-100) (se filtraggio manuale)
  33. # non sono stati impostati i filtri per poligono e per selezione manuale dalla mappa
  34.  
  35. ## Formato di output
  36. interp_res <- 5 # Risoluzione (m)
  37. interp_crs <- st_crs(32632) # Sistema di riferimento
  38. # Le opzioni si riferiscono all'impstazione "Definisci risoluzione e proiezione";
  39. # per scegliere da un raster esistente, oltre a scegliere risoluzione e CRS
  40. # del raster bisognerebbe impostare anche "grid_offset" -
  41. # lo documenterò all'occorrenza.
  42.  
  43. ## Opzioni di interpolazione
  44. interp_maxptdist <- 15 # Distanza massima dai punti
  45. interp_method <- "krige" # Metodo di interpolazione ("krige" o "idw")
  46. # Definizione del variogramma: l'esempio funziona nel caso di "Automatica";
  47. # altre modalità richiedono la definizione manuale del variogramma con metodi
  48. # e funzioni del pacchetto {gstat}.
  49. focal_onoff <- TRUE # Livella i raster interpolati
  50.  
  51. ## Opzioni di processamento
  52. interp_parallel <- TRUE # Ottimizza l'interpolazione per uso desktop (FALSE) o server/workstation (TRUE)
  53. out_path <- "/path/of/raster_out.tif" # Raster di output
  54. # Opzioni avanzate
  55. v_maxdist <- NA # Distanza massima punti-pixel (auto: NA)
  56. v_nmax <- 100 # Numero di punti per pixel (auto: 500)
  57. interp_samplesize <- NA # Usa un campione di punti per l'interpolazione (NA/Inf per usarli tutti)
  58. interp_samplescheme <- "strat_area" # Schema di campionamento nei poligoni ("random", "strat_npts", "strat_area" o "strat_prop")
  59.  
  60.  
  61. ### Esecuzione ----
  62.  
  63. # Carica i file vettoriali
  64. borders_sf0 <- st_read(borders_path) %>%
  65.   st_transform(4326)
  66. inputpts_sf <- st_read(inputpts_path) %>%
  67.   st_transform(4326)
  68.  
  69. # Converti il vettoriale dei bordi nel formato usato da guinterp
  70. borders_sf <- if (select_uid_which == "attr") {
  71.   borders_sf0 %>%
  72.     dplyr::select(select_uid_attr) %>%
  73.     dplyr::rename(id_geom = select_uid_attr)
  74. } else if (select_uid_which == "record") {
  75.   borders_sf0 %>%
  76.     dplyr::transmute(id_geom = seq_len(nrow(borders_sf0)))
  77. } else if (select_uid_which == "no") {
  78.   borders_sf0 %>%
  79.     dplyr::transmute(id_geom = "all points")
  80. } %>%
  81.   dplyr::group_by(id_geom) %>% dplyr::summarise()
  82.  
  83. # Converti il vettoriale dei punti nel formato usato da guinterp
  84. inputpts_dt <- read_inputpts(
  85.   inputpts_sf,
  86.   borders = borders_sf,
  87.   varname = select_inputvar
  88. )
  89.  
  90. # Filtra i punti
  91. if (filter_buttons == "manual") { # applica i filtri manuali
  92.   inputpts_dt <- filter_pts(
  93.     inputpts_dt, "rangey", filter_rangey,
  94.     id_fieldname="id_geom",
  95.     inlayer = borders_sf,
  96.     byfield = TRUE, samplesize = NA
  97.   )
  98.   inputpts_dt <- filter_pts(
  99.     inputpts_dt, "zscorey", filter_zscorey/100,
  100.     id_fieldname="id_geom",
  101.     inlayer = borders_sf,
  102.     byfield = TRUE, samplesize = NA
  103.   )
  104.   inputpts_dt <- filter_pts(
  105.     inputpts_dt, "rbiasy", filter_rbiasy/100,
  106.     id_fieldname="id_geom",
  107.     inlayer = borders_sf,
  108.     byfield = TRUE, samplesize = NA
  109.   )
  110.   inputpts_dt <- filter_pts(
  111.     inputpts_dt, "rangeq", filter_rangeq/100,
  112.     id_fieldname="id_geom",
  113.     inlayer = borders_sf,
  114.     byfield = TRUE, samplesize = NA
  115.   )
  116.   inputpts_dt <- filter_pts(
  117.     inputpts_dt, "rangey", filter_rangey,
  118.     id_fieldname="id_geom",
  119.     inlayer = borders_sf,
  120.     byfield = TRUE, samplesize = NA
  121.   )
  122. } else if (filter_buttons == "minimal") { # applica il filtro minimo (2°-98° percentile)
  123.   inputpts_dt <- filter_pts(
  124.     "rangeq", c(.02,.98),
  125.     inlayer = borders_sf,
  126.     byfield = TRUE, samplesize = NA
  127.   )
  128. }
  129.  
  130. # Interpola
  131. out_path <- guinterp_process(
  132.   inputpts_dt, borders_sf,
  133.   id_fieldname = "id_geom",
  134.   out_path = out_path,
  135.   samplesize = interp_samplesize,
  136.   samplescheme = interp_samplescheme,
  137.   parallel = interp_parallel,
  138.   interp_method = interp_method,
  139.   smooth = focal_onoff,
  140.   interp_res = interp_res,
  141.   out_crs = interp_crs,
  142.   buffer_radius = interp_maxptdist,
  143.   vgm = "auto",
  144.   v_nmax = v_nmax,
  145.   v_maxdist = v_maxdist
  146. )
  147.  
  148. # Visualizza output
  149. mapview::mapview(raster(out_path))
  150.  
RAW Paste Data

Adblocker detected! Please consider disabling it...

We've detected AdBlock Plus or some other adblocking software preventing Pastebin.com from fully loading.

We don't have any obnoxious sound, or popup ads, we actively block these annoying types of ads!

Please add Pastebin.com to your ad blocker whitelist or disable your adblocking software.

×