In lbusett/phenoriceR: Routines for analysis of PhenoRice outputs

knitr::opts_chunk$set(echo = TRUE)
library(sprawl)
library(tibble)
library(raster)
library(phenoriceR)
library(dplyr)
library(sf)
library(dplyr)

#   ____________________________________________________________________________
#   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"))

#'  # Suppose you want to extract data for Region: "Region_3_-_Central_Luzon" :
#'  # --> extract it from the full shapefile

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)
# this is needed to remove duplicate polygons and keep pnly useful columns
# of the vector
in_vect <- unique(in_vect[c(1:4, 19)]) 
# reproject to the CRS of the rasters
in_vect <- reproj_vect(in_vect,get_proj4string(
    file.path(main_folder,"orig_mosaic/param_series/decirc/eos_decirc.tif")))

# in_vect contains one polygon for each sub_region of "Region_3_-_Central_Luzon"
plot_vect(in_vect, fill_var = "ID_name")

message("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)
out_folder

# you have all you need to crop the mosaics: folder of the mosaic rasters, a 
# vector where to cut and an output folder. Cropped rasters will be available 
# in out_folder

message("Extracting data on ", Region_name)
pr_extract_subarea(in_mosaics_folder = in_rast_folder,
                   in_mask           = 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

#### test on one file: in this case it is eos_decirc.RData ----
in_rast <- get(load(in_files[1]))
message("Working on : ", in_files[3])
out <- sprawl::extract_rast(in_rast,
                     in_vect,
                     na.value = 0, 
                     join_geom = FALSE, 
                     id_field = "ID_name", 
                     verbose = FALSE)
out


## 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 = "_")))