# Title : Groundwater Sustainability Policy for the United States -GSP4
# By: Iman Haqiqi

# Read required libraries
library(raster)

# set the working directory 
setwd("c:/h/")

# ---- load the input files ----
# This is the irrigated area from USGS MIrAD aggregated to 5 arc-min
mirad_5min = raster("USGS_MIrAD_5min_2010.tif")

# This is the groundwater extraction to recharge rates from Reitz et al (2017) 
rate_5min  = raster("Reitz_rate_5min_2010.tif")

# This is the grid IDs of the SIMPLE-G-US model
grid_5min  = raster("SIMPLEG_GRID_5min.tif")


# excluding the grid cells outside MIrAD
X2R <- (mirad_5min/mirad_5min)*rate_5min

# Resample to ensure consistency
X2R <- resample(X2R, grid_5min)

# calculate the sustainability shock and save as CSV
st = stack(grid_5min, X2R)

st$shock = 100*(1- st$layer)/st$layer

st$shock[st$shock > 0] = 0

df <- as.data.frame(rasterToPoints(st))

# Write CSV 
write.csv(df, file = "GSP4_GroundwaterSustainabilityPolicy.csv")

writeRaster(st$shock,
            "GSP4_GroundwaterSustainabilityPolicy_pct.tif",
            format="GTiff",
            overwrite=T)

writeRaster(st$shock,
            "GSP4_GroundwaterSustainabilityPolicy_pct.nc",
            format="CDF",
            overwrite=T,
            varname="WAT4_QWATGRD",
            longname="GSP4 Groundwater Sustainability Policy in percentage")