#   ____________________________________________________________________________
#   set input and output folders                                            ####

# set the main folder ----
main_folder <- "/home/lb/my_data/prasia/Data"

# Load the reshuffled shapefile ----

in_shp <- read_vect(file.path(main_folder, 
                              "vector/Ricetlas/riceatlas_asia_reshuffled.shp"))

## Set a region for the analysis ----
Region_name <- "Region_3_-_Central_Luzon"

## "1) Get the polygons of a specific region from the shapefile") ----

in_vect <- dplyr::filter(in_shp, Region == Region_name)
in_vect <- in_vect[1:3] %>%
  sf::st_transform(get_proj4string(
    file.path(main_folder,"orig_mosaic/param_series/decirc/eos_decirc.tif"))) %>% 
  unique()
plot_vect(in_vect, fill_var = "ID")

cat("Create cropped rasters and put them in the \"subsets\" subfolder")
## Create cropped rasters and put them in the "subsets" subfolder
in_rast_folder  <- file.path(main_folder, "orig_mosaic/param_series/")
out_folder      <- file.path(main_folder, "subsets", Region_name)
make_folder(out_folder, type = "dirname", verbose = T)

message("Extracting data on ", Region_name)
pr_extract_subarea(in_rast_folder,
                   in_vect,
                   out_folder  = out_folder)

## 2. Extract the data from the different provinces of the region" ----

in_files <- list.files(
  file.path(out_folder, "param_series/decirc"),
  pattern = "*.RData", 
  full.names = TRUE)
in_files
## [1] "/home/lb/my_data/prasia/Data/subsets/Region_3_-_Central_Luzon/param_series/decirc/eos_decirc.RData"
## [2] "/home/lb/my_data/prasia/Data/subsets/Region_3_-_Central_Luzon/param_series/decirc/pos_decirc.RData"
## [3] "/home/lb/my_data/prasia/Data/subsets/Region_3_-_Central_Luzon/param_series/decirc/sos_decirc.RData"
#### test on one file ----
in_rast <- get(load(in_files[1]))
message("Working on : ", in_files[1])
out <- sprawl::extract_rast(in_rast,
                     in_vect,
                     na.value = 0, 
                     join_geom = FALSE, 
                     id_field = "ID", 
                     verbose = FALSE)
out
## $stats
## # A tibble: 392 x 12
##         ID band_n       date n_pix n_pix_val      avg   med       sd   min
##     <fctr>  <dbl>     <date> <int>     <int>    <dbl> <dbl>    <dbl> <dbl>
##  1 PHL_843      1 2003-01-01 68393        33 1101.455  1100 14.48510  1076
##  2 PHL_843      2 2003-01-01 68393        84 1239.000  1249 24.55189  1169
##  3 PHL_843      3 2003-01-01 68393       885 1340.471  1345 28.24250  1257
##  4 PHL_843      4 2003-01-01 68393       807 1387.587  1385 18.32833  1353
##  5 PHL_843      5 2004-01-01 68393        56 1479.714  1477 19.81892  1433
##  6 PHL_843      6 2004-01-01 68393        70 1600.286  1606 27.68411  1534
##  7 PHL_843      7 2004-01-01 68393       700 1687.131  1686 27.46557  1622
##  8 PHL_843      8 2004-01-01 68393       619 1758.388  1750 23.07280  1718
##  9 PHL_843      9 2005-01-01 68393        42 1845.095  1839 18.26073  1815
## 10 PHL_843     10 2005-01-01 68393        17 1953.176  1948 24.32138  1924
## # ... with 382 more rows, and 3 more variables: max <dbl>, OBJECTI <dbl>,
## #   cty <fctr>
## 
## $alldata
## # A tibble: 638,940 x 8
##         ID band_n       date n_pix_val     N value OBJECTI    cty
##     <fctr>  <int>     <date>     <int> <int> <dbl>   <dbl> <fctr>
##  1 PHL_843      1 2003-01-01        33     1  1124     843    PHL
##  2 PHL_843      1 2003-01-01        33     2  1084     843    PHL
##  3 PHL_843      1 2003-01-01        33     3  1092     843    PHL
##  4 PHL_843      1 2003-01-01        33     4  1084     843    PHL
##  5 PHL_843      1 2003-01-01        33     5  1084     843    PHL
##  6 PHL_843      1 2003-01-01        33     6  1108     843    PHL
##  7 PHL_843      1 2003-01-01        33     7  1108     843    PHL
##  8 PHL_843      1 2003-01-01        33     8  1108     843    PHL
##  9 PHL_843      1 2003-01-01        33     9  1092     843    PHL
## 10 PHL_843      1 2003-01-01        33    10  1100     843    PHL
## # ... with 638,930 more rows
## 3. Save results as an RData file for future use" ----
make_folder(file.path(out_folder, "RData"))
save(out, file = file.path(out_folder, "RData",
                           paste(Region_name, "stats.RData", sep = "_")))