library(devtools)
## Loading required package: usethis
library(cwatmRutils)
library(raster)
## Loading required package: sp
## The legacy packages maptools, rgdal, and rgeos, underpinning this package
## will retire shortly. Please refer to R-spatial evolution reports on
## https://r-spatial.org/r/2023/05/15/evolution4.html for details.
## This package is now running under evolution status 0
library(ggplot2)
library(ncdf4)

Irrigation data: https://www.fao.org/aquastat/en/geospatial-information/global-maps-irrigated-areas/latest-version/ Stefan Siebert, Verena Henrich, Karen Frenken and Jacob Burke (2013). Global Map of Irrigation Areas version 5. Rheinische Friedrich-Wilhelms-University, Bonn, Germany / Food and Agriculture Organization of the United Nations, Rome, Italy

Irrigation for year 2005

## FAO Irrigation data
fao_crop = raster("H:/MyDocuments/Validation data/gmia_v5_aei_pct.asc")
plot(fao_crop)

fao_crop_malawi <- crop(fao_crop, extent(31.5  , 37, -18.5 , -7.5 ), snap="out")
#writeRaster(fao_crop_malawi,"H:/MyDocuments/Validation data/fao_crop_percent.tif")
plot(fao_crop_malawi)

fraction_land_cover <- ncdf2raster(pth= "C:/Users/hinton/Documents/GitHub/MalawiLake/cwatm_data/maps/landsurface/fractionLandcover.nc", transpose = TRUE)
                                   
fraction_land_cover <- ncdf2raster(pth= "C:/Users/hinton/Documents/GitHub/MalawiLake/cwatm_data/maps/landsurface/fractionLandcover.nc", transpose = TRUE, time = c(1,50), na.rm= TRUE, varName= c('fracirrNonPaddy', 'fracirrPaddy'))

### 2005 is timestep "37986"
### Extract irrigation estimate for non_paddy irrigation for given time period
non_paddy_irr<- (fraction_land_cover$`21549`$fracirrNonPaddy)*100
non_paddy_irr_later<- (fraction_land_cover$`39447`$fracirrNonPaddy)*100
non_paddy_irr_2005<- (fraction_land_cover$`37986`$fracirrNonPaddy)*100
paddy_irr_2005<- (fraction_land_cover$`37986`$fracirrPaddy)*100 

## All irrigation types summed
all_irr_2005<- paddy_irr_2005+non_paddy_irr_2005

### Compare growth of irrigation over time
frac_irr_times <- ncdf2raster(pth= "C:/Users/hinton/Documents/GitHub/MalawiLake/cwatm_data/maps/landsurface/fractionLandcover.nc", transpose = TRUE, time = c(1,50), fun= mean, na.rm= TRUE, varName= c('fracirrNonPaddy', 'fracirrPaddy'))
frac_irr_times$date<- as.Date(frac_irr_times$time, origin= "1901-01-01")
ggplot(frac_irr_times, aes(date, value, colour=var))+geom_line()

###
# Worldbank estimation of irrigated land
#Irr_percent_worldbank<- read.csv('H:/MyDocuments/Validation data/Irrigation/Irrigation_percent_agricultural_land.csv')
Irr_percent_worldbank_malawi<- c(0.579150579150579,0.529801324503311,0.74222668004012,0.533578656853726)

Compare Fraction land cover irrigation estimate for “all irrigation” for 2005 with FAO 2005 estimate

cellStats(all_irr_2005,mean)
## [1] 0.209814
cellStats(fao_crop_malawi, mean)
## [1] 0.2281341
## Mean for FAO estimate slightly higher. Both report very low % of irrigation (0.2%).
# Calculate root mean square error
# Need to align the 2 rasters 
all_irr_cropped_2005<- crop(all_irr_2005, fao_crop_malawi)
all_irr_cropped_2005_resample <- resample(all_irr_cropped_2005, fao_crop_malawi)
# Define RMSE function
RMSE <- function(x, y) { sqrt(mean(!is.na((x - y)^2))) } 

#Remove 0 values
values_all_irr_2005<- values(all_irr_cropped_2005_resample)
values_all_irr_2005[values_all_irr_2005==0]<- NA

values_fao_crop_malawi<- values(fao_crop_malawi)
values_fao_crop_malawi[values_fao_crop_malawi==0]<- NA

RMSE (values_all_irr_2005,values_fao_crop_malawi)
## [1] 0.3524153
#################################################
#Increase all_irr_cropped_2005_resample to match same mean as FAO (8.7%)
all_irr_cropped_2005_resample_adjust<- all_irr_cropped_2005_resample*((cellStats(fao_crop_malawi, mean))/cellStats(all_irr_2005,mean))

values_all_irr_2005_adjust<- values(all_irr_cropped_2005_resample_adjust)
values_all_irr_2005_adjust[values_all_irr_2005_adjust==0]<- NA

RMSE (values_all_irr_2005_adjust,values_fao_crop_malawi)
## [1] 0.3524153
##################################################
diff_adjust <- overlay(fao_crop_malawi, all_irr_cropped_2005_resample_adjust, fun=function(a,b) return(a==b))
plot(diff_adjust,
     col=c('#FFE4E1','#228B22'),
     legend=FALSE,
     axes=FALSE)
legend("left", legend=c("Agree", "Disagree"),
       col=c("#228B22", "#FFE4E1"), pch = 15, cex=0.8)

#writeRaster(difference_quant_adjusted, "H:/MyDocuments/Validation data/diff_quant_all_irr_cropped_2005_resample_adjust.tif")

#difference_quant_adjusted<- (all_irr_cropped_2005_resample_adjust-fao_crop_malawi)
#writeRaster(difference_quant, "H:/MyDocuments/Validation data/irrigation_diff_quant.tif",overwrite=TRUE)