Step 03: Extraction of data from raster bricks

##############################################################
######## Project Name : SATSURE Credit Risk Product
######## File Usage : This script generates CSV files from raster images
######## Date : November 20, 2018
######## Author : Actify Data Labs
##############################################################

# How to run the script
# Rscript feature_csv_generate.R <veg_index> <crop_type> <district>

############################################################################################################
# Inputs to be given to this file : feature_csv_generate.R
#
# vegetation_index: name of the vegetation index – ndvi/ndwi/ndre .. etc
# crop_type: maize/ragi/paddy.. etc
# district: name of district – Srikakulum/Vizianagaram .. etc
# years: years for which the extraction is to be done (in character format)
# for eg: years <- seq(from = 2013, to = 2017, by = 1) %>% as.character
#
# shp_dir: Shapefile Directory( filepath where all the mandal shapefiles reside)
# data_dir: Data directory (filepath where all the data corresponding to a particular district resides)
# years_of_missing_data: Year for which data is missing
###########################################################################################################

# rm(list = ls())

process_start_time <- Sys.time()
args <- commandArgs(TRUE)

source(“~/codes/credit_risk_product/utils.R”)
source(“~/codes/credit_risk_product/settings.R”)

# Importing necessary libraries
import_libraries()

# ###### for testing purposes only #######
# vegetation_index <- “ndvi” %>% toupper
# crop_type <- “maize”
# mandal_names <- mandal_names[7]
# year = “2017”
# #######################################

####### Parameters
(vegetation_index <- toupper(args[1]))
(crop_type <- args[2])
(district <- tools::toTitleCase(args[3]))
(shp_dir <- glue::glue(“~/{district}/SHP/mandal_shapefiles”))
(data_dir <- glue::glue(“~/{district}”))
years <- as.character(2013:2017)
names(years) <- paste0(“X”, years)
# (year_missing_data <- c(2015) %>% as.character())
###################

# # removing years of missing data
# (years <- setdiff(years, year_missing_data))

# Getting vector layer names and shapefiles
(vector_layer_names <- read_shape_files(dir = shp_dir, file_pattern = “geojson$|shp$”))
(mandal_names <- vector_layer_names$mandal_names)

# system.time(lapply(years[length(years)], function(year){

for (year in years){

print(glue::glue(“Extraction started for year {year}..”))
year_start_time <- Sys.time()
# Reading yearwise raster files
(dir_name <- file.path(data_dir, year, toupper(vegetation_index), “mandal_level_bricks”))

# # Making district level raster bricks (Old Approach) ———————-
#
#
# print(glue::glue(“Reading brick for year {year} and vegetation index {vegetation_index}”))
# (brick_file_path <- glue::glue(“{data_dir}/{year}/{vegetation_index}/rast_brick_{year}_{vegetation_index}.grd”))
# (raster_brick <- raster::brick(brick_file_path))

# # Making raster brick
# print(paste0(“Making raster brick for the year “,year, “…….”))
# system.time(raster_brick <- brick(raster_list))

# Reading the crop raster ———————————————

print(glue::glue(“Reading the {crop_type} raster for the year {year}…..”))

(crop_folder <- list.files(glue::glue(“{data_dir}/{year}”), full.names = TRUE) %>% grep(pattern = crop_type, x = ., ignore.case = TRUE, value = TRUE))

(crop_filename <- list.files(path = crop_folder, pattern = “img|tif”, full.names = TRUE) %>% grep(pattern = crop_type, x = ., ignore.case = TRUE, value = TRUE))

(plant_raster <- raster(crop_filename))

# proj4string(plant_raster) <- CRS(proj4string(raster_brick))

# Creating a folder for yearly results
(result_folder_path <- file.path(results_dir, district, year, toupper(vegetation_index), crop_type))

if(!dir.exists(result_folder_path)){
dir.create(result_folder_path, recursive = TRUE)
}

# lapply(tail(mandal_names,-4), function(mandal_name){
for (mandal_name in mandal_names){

print(glue::glue(“Extraction started for year {year} and mandal {mandal_name}”))
mandal_start_time <- Sys.time()

# Reading mandal shapefiles (for cropping the crop raster) ———————————————-

(mandal_name_only <- unlist(strsplit(mandal_name, split = ‘.’, fixed = TRUE))[1])
print(paste0(“Reading mandal shapefiles for “, mandal_name_only))
(mandal_file_name <- file.path(shp_dir, mandal_name))
(mandal_shape_file <- readOGR(mandal_file_name))

# Reading mandal raster bricks ——————————————-

raster_brick <- raster::brick(glue::glue(“{dir_name}/{mandal_name_only}_raster_brick.grd”))

# making the CRS of plant raster and raster_brick the same
proj4string(plant_raster) <- CRS(proj4string(raster_brick))

print(“Projecting the shapefile to same CRS as the vegetation index raster …”)
(mandal_shape_file <- spTransform(mandal_shape_file, CRS(proj4string(raster_brick))))
print(“Projected to same CRS!!!”)

print(paste0(“Extracting vegetation index values for “,mandal_name_only, ” ….”))
system.time(index_values <- extract_index_value(vegetation_index_raster_brick = raster_brick,
mandal_shape_file = mandal_shape_file,
district_plant_raster = plant_raster
))
print(paste0(“Index value extraction for “,mandal_name_only,” complete!”))

(file_name <- file.path(result_folder_path, paste0(mandal_name_only,”_”,crop_type,”_”, toupper(vegetation_index), “.csv”)))
# print(file_name)
# print(gc())
print(paste0(“Writing Pixel Matrix for “, mandal_name_only, “…………”))
data.table::fwrite(index_values,
file = file_name,
sep = “,”,
col.names = TRUE,
row.names = FALSE,
append = FALSE)

print(glue::glue(“Successfully written Pixel Matrix for mandal {mandal_name_only}!!”))

print(glue::glue(“Time taken for year {year} and mandal {mandal_name_only} is {as.numeric(Sys.time()-mandal_start_time)}”))

cat(“\n\n”)

} # end of mandal loop

print(glue::glue(“Time taken for year {year} is {as.numeric(Sys.time()-year_start_time)}”))

} # end of year loop

print(glue::glue(“Process completed for district {district}!!”))

print(glue::glue(“The total time taken for this process is {as.numeric(Sys.time()-process_start_time)}”))