Evaluatuion of soil moisture from CCAMM-CABLE simulation, satellite-based models estimates and statellite observations Skukuza and Malopeni fluxters region case study
Floyd V. Khosa, Mohau J. Mateyisi
08/02/2020
#data description & repository: https://data.mendeley.com/datasets/yvsncwrbws/draft?a=4e0ea427-19f4-480e-a4d0-935487a634e2
The following code reads in a shapefile and deliniate the study region
## Loading required package: sp## rgeos version: 0.5-2, (SVN revision 621)
## GEOS runtime version: 3.6.1-CAPI-1.10.1
## Linking to sp version: 1.3-2
## Polygon checking: TRUE## rgdal: version: 1.4-8, (SVN revision 845)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
## Path to GDAL shared files: C:/Users/MMateyisi/Documents/R/win-library/3.6/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
## Path to PROJ.4 shared files: C:/Users/MMateyisi/Documents/R/win-library/3.6/rgdal/proj
## Linking to sp version: 1.3-2## OGR data source with driver: ESRI Shapefile
## Source: "D:\HESS-2018-546\Study-region\SA_GADMIN", layer: "ZAF_adm1"
## with 10 features
## It has 16 fieldssection 2.4.3
The following codes produces a plot of the elevation for the study region
The paths may be replaced to direct the programme to the downloaded datasets
library(ggplot2)
library(rgeos)
library(ncdf4)
library(raster)
library(rgdal)
library(scales)
Sub.regions <- read.csv("D:\\HESS-2018-546\\Study-region\\selected-subregions.csv") # sub-region data
Topography.25km <- raster("D:\\HESS-2018-546\\Study-region\\Topography-25km.tif")
study.region.shapefile<-"D:\\HESS-2018-546\\Study-region\\SA_GADMIN"
plot_elevation <- function(Sub.regions,Topography.25km, study.region.shapefile){
#======================================================================================
# The functions inputs: Topography data at 25km resolution
#
#
# Ouput: elevation for (Figure 2a) with marked sub-regions
#=====================================================================================
region.bound <- readOGR(dsn=path.expand(study.region.shapefile), layer="ZAF_adm1")
#calculates rectengular boundry box encompassing your region
region.env <- gEnvelope(region.bound)
#extract exterior of your region within the boundrary
region.diff <- gDifference(region.env,region.bound)
#convert your exterior region to a data frame
region.diff.df <- fortify(region.diff)
#convert your region shape file to data frame
OVERLAY <- fortify(region.bound, region="Shape_Area")
# convert raster to data frame
Out.DTF <- as.data.frame(matrix(nrow=0,ncol=3))
names(Out.DTF) <-c("long","lat","z")
test_spdf <- as(Topography.25km, "SpatialPixelsDataFrame")
Topography.25km <- as.data.frame(test_spdf)
names(Topography.25km) <-c("Alt","long","lat")
# show only positive values of altitude
Topography.25km <- Topography.25km[(Topography.25km$Alt >= 0),]
Sub.regions$area <- factor(Sub.regions$area,levels = c("a","b","c","d","e","f"))
p <- ggplot()+
geom_raster(data = Topography.25km, mapping = aes(x=long,y=lat,fill=Alt),interpolate = F) +
geom_point(data = Sub.regions, aes(x=lon,y=lat,shape=area))+
geom_polygon(data = OVERLAY, aes(x=long,y=lat, group = group),colour = "black",fill= NA)+
coord_equal(xlim =c(24,33) ,ylim=c(-29,-22.5))+
geom_polygon(data=region.diff.df ,aes(x=long,y=lat,group=group),colour="white",size=1,fill="White")+
labs(shape="Sub-Region",colour="Alt(m)",fill="Alt(m)",x ="longitude",y="latitude")+
scale_fill_gradientn(colours=c("blue","cyan","white", "green","yellow","red"),
values=rescale(c(-2,0-.Machine$double.eps,0,0+.Machine$double.eps,2)))+
guides(fill = guide_colorbar(barwidth = 0.5, barheight = 10))+
theme(strip.background = element_blank())+
theme(strip.text.x = element_text(size = 16,colour="black"))+#
theme(strip.text.y = element_text(size = 16,colour="black"))+#
theme(strip.background = element_blank())+
scale_shape_manual(values=c(0,1,2,3,4,5),labels=c("OceaSa","HuSuSa-HoSeSa","HoSeSa","HoSeGr"
,"CoSeSa(I)",
"CoSeGr(II)"))+
theme(text=element_text(colour="black",size=12))+
theme(panel.border = element_rect(fill = NA, colour = "black"))#+
#guides(shape=TRUE)
p
}
plot_elevation(Sub.regions,Topography.25km, study.region.shapefile)## OGR data source with driver: ESRI Shapefile
## Source: "D:\HESS-2018-546\Study-region\SA_GADMIN", layer: "ZAF_adm1"
## with 10 features
## It has 16 fields
R Markdown
#Altitude (Alt, m) at the study region at a 25 km resolution, (https://www.ngdc.noaa.gov/mgg/fliers/06mgg01.html)
#The following code plots soil types for the study region
This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.
When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:
library(ggplot2)
library(rgeos)
library(rgdal)
soil.types <- read.csv("D:\\HESS-2018-546\\Study-region\\soil-types.csv")
fluxTower.sites <- read.csv("D:\\HESS-2018-546\\Study-region\\coordinates_validation_sites.csv")
Sub.regions <- read.csv("D:\\HESS-2018-546\\Study-region\\selected-subregions.csv") # sub-region data
study.region.shapefile<-"D:\\HESS-2018-546\\Study-region\\SA_GADMIN"
Plot_soil_types <-function(soil.types,Sub.regions,fluxTower.sites,study.region.shapefile){
#======================================================================================
# The functions inputs: longitute, latitude as well atrributes (1) "soil types", (2) "sub.regions", (3)
# "fluxTower.sites" and (4) a shape file of the study region "SA_GADMIN"
#
# Ouput: Soiltypes (Figure 2b)
#=====================================================================================
region.bound <- readOGR(dsn=path.expand(study.region.shapefile), layer="ZAF_adm1")
region.env <- gEnvelope(region.bound) #calculates rectengular boundry box encompassing your region
region.diff <- gDifference(region.env,region.bound) #extract exterior of your region within the boundrary box
region.diff.df <- fortify(region.diff) #convert your exterior region to a data frame
OVERLAY <- fortify(region.bound, region="Shape_Area") #convert your region shape file to data frame
#introduce categorical variable soil moisture type for plotting
soil.types$ST <- replace(soil.types$ST,soil.types$ST==1,"sand")
soil.types$ST <- replace(soil.types$ST,soil.types$ST==2,"loamy sand")
soil.types$ST <- replace(soil.types$ST,soil.types$ST==3,"sandy loam")
soil.types$ST <- replace(soil.types$ST,soil.types$ST==4,"sandy clay")
soil.types$ST <- replace(soil.types$ST,soil.types$ST==5,"loam")
soil.types$ST <- replace(soil.types$ST,soil.types$ST==6,"silt loam")
soil.types$ST <- replace(soil.types$ST,soil.types$ST==7,"silt")
soil.types$ST <- replace(soil.types$ST,soil.types$ST==8,"silty clay loam")
soil.types$ST <- replace(soil.types$ST,soil.types$ST==9,"clay")
soil.types$ST <- replace(soil.types$ST,soil.types$ST==10,"clay loam")
soil.types$ST <- replace(soil.types$ST,soil.types$ST==11,"sandy clay loam")
soil.types$ST <- replace(soil.types$ST,soil.types$ST==12,"silty clay")
soil.types$ST <- factor(soil.types$ST,levels = c("sand","loamy sand","sandy loam","sandy clay","loam","silt loam",
"silt","silty clay loam","clay","clay loam","sandy clay loam","silty clay"))
#LABELS<-c("OceaSa","HuSuSa-HoSeSa","HoSeSa","HoSeGr","CoSeSa(I)","CoSeGr(II)")
#soil types plotting function
ggplot()+
geom_raster(data =soil.types, mapping = aes(x=long,y=lat,fill=ST),interpolate = F) +
geom_point(data = fluxTower.sites, aes(x = longitude, y = latitude), size = 4, shape = 20)+
geom_point(data = Sub.regions, aes(x=lon,y=lat,shape=area))+
geom_polygon(data = OVERLAY, aes(x=long,y=lat, group = group),colour = "black",fill= NA)+
coord_equal(xlim =c(24,33) ,ylim=c(-29,-22.5))+
#scale_shape_manual(values=c(0,1,2,3,4,5))+
scale_fill_manual(values = c("sand"="red",
"loamy sand"="#ff00ff",
"sandy loam"="#FAEBD7",
"sandy clay"="navy",
"loam"="blue",
"silt loam"="#00a0a0",
"silt"="#00eeee",
"silty clay loam"="#00ee00",
"clay"="gray",
"clay loam"="#707030",
"sandy clay loam"="yellow",
"silty clay"="maroon"),name="Soil Types",guide = guide_legend(reverse = T))+
geom_polygon(data=region.diff.df ,aes(x=long,y=lat,group=group),colour="white",size=1,fill="White")+
labs(shape="sub-region", fill="Soil types",x ="longitude",y="latitude")+
theme(axis.text.x = element_text(colour = "black"))+
theme(axis.text.y = element_text(colour = "black"))+
theme(panel.border = element_rect(fill = NA, colour = "black"))+
scale_shape_manual(values=c(0,1,2,3,4,5),labels=c("OceaSa","HuSuSa-HoSeSa","HoSeSa","HoSeGr"
,"CoSeSa(I)","CoSeGr(II)"))
#guides(shape=FALSE)
}
Plot_soil_types(soil.types,Sub.regions,fluxTower.sites,study.region.shapefile)## OGR data source with driver: ESRI Shapefile
## Source: "D:\HESS-2018-546\Study-region\SA_GADMIN", layer: "ZAF_adm1"
## with 10 features
## It has 16 fields
Including Plots
#Dominant soil types (https://soilgrids.org) per grid cell, at a resolution of 25 km
#The following code plot climate types for the study region
You can also embed plots, for example:
## OGR data source with driver: ESRI Shapefile
## Source: "D:\HESS-2018-546\Study-region\SA_GADMIN", layer: "ZAF_adm1"
## with 10 features
## It has 16 fields
Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.
Section 3.1.2
The following code produces the Taylor and Normalised mean bias (NMB)
library(openair)
Monthly.soil.moisture <- read.csv("D:\\HESS-2018-546\\Taylor & NMB\\Taylor-plot.csv")
get_tylor_plots<-function(LongCData){
#===========================================================================================================
# The funcion plots a Taylor diagram
# Arguments: dataframe
# input: date, in situ observations, modelled data, site and depth
# Output: Taytor diagram
#===========================================================================================================
LongCData$date=as.Date(LongCData$date, format = "%d/%m/%Y")
LongCData$site <- factor(LongCData$site, levels = c("Skukuza","Malopeni"))
TaylorDiagram(LongCData, obs = "in.situ", mod = "value",group = "variable", type = c("depth","site"),
key = TRUE,ylab="Standard deviation (%)",key.title = "Model", key.pos = "right",rms.col="black",
cols = c("#009E73","#0072B2","#D55E00","#E69F00","#F0E422"))
}
get_tylor_plots(Monthly.soil.moisture)
### Normalised mean bias
normlised_meanbias <- read.csv("D:\\HESS-2018-546\\Taylor & NMB\\normalised-mean-bias.csv")
get_normalised_mean_bias <-function(normlised_meanbias){
#=============================================================================================
# The function plots normalised mean bias (NMB)
# Arguments: Dataframe
# Input: Model name, "site", "depth" and "NMB"
# Regurns: normalised mean bias pllot
#============================================================================================
normlised_meanbias$site <- factor(normlised_meanbias$site, levels = c("Skukuza","Malopeni"))
cbpalette <- c("#009E73","#0072B2","#D55E00","#E69F00","#F0E422")
p<-qplot(variable, NMB, data = normlised_meanbias, color = variable, size=0.2)+
geom_hline(yintercept=0, linetype="dashed", color = "black", size=1)+
facet_grid(site~depth,scales = "free_x")+#
labs(y="Normalised mean bias", x = "Models",color="Models")+
theme(legend.position = "right", legend.box = "vertical") +
scale_color_manual(values = cbpalette)+
theme(axis.text.y = element_text(),axis.title.x = element_text(),
axis.title.y = element_text())+
theme(panel.grid.major = element_line(size = 0.5, linetype = 'solid',
colour = "grey"),
panel.grid.minor = element_line(size = 0.5, linetype = 'solid',
colour = "grey"),
panel.background = element_blank())+
theme(panel.border = element_rect(fill = NA, colour = "black"))+
theme(strip.text.x = element_text())+
theme(axis.text.x = element_text(colour = "black"))+
theme(strip.background = element_blank())+
guides(size=FALSE)+
guides(colour = guide_legend(override.aes = list(size=4, stroke=1.5)))+
theme(axis.text.x = element_text(angle = 90, hjust = 1,vjust=0.5))+
#theme(axis.title.x=element_blank(),axis.text.x=element_blank())+
theme(text=element_text())
p
}
get_normalised_mean_bias(normlised_meanbias)
3.1.3
The following code computes and plot crosswavelets
library("WaveletComp")##
## Attaching package: 'WaveletComp'## The following object is masked from 'package:ggplot2':
##
## arrowlibrary(date)
# soil moisture files: SSM-comparison-all-products-skukuza-monthly-imputed.csv,
# RZSM-comparison-all-products-skukuza-monthly-imputed.csv"
SSMmonthly.soil.moisture<- read.csv("D:\\HESS-2018-546\\CrossWavelet\\SSM-comparison-all-products-skukuza-monthly-imputed.csv")
RZSmonthly.soil.moisture<-read.csv("D:\\HESS-2018-546\\CrossWavelet\\RZSM-comparison-all-products-skukuza-monthly-imputed.csv")
model.name<-"CCAM.CABLE"
getWaveltAnalysys<-function(monthly.soil.moisture,model.name){
#*****************************************************************************************
# The function calculates and plots crosswavelets
# Arguments: dataframe and model name
# Inputs: date,insitu soil moisture, modelled soil moisture
# Returns: cross-wavelet plots
# possible model names c("CCAM-CABLE","GLEAM3a")
#****************************************************************************************
monthly.soil.moisture$date=as.POSIXct(monthly.soil.moisture$date, tz = "Africa/Johannesburg", format = "%d/%m/%Y")
#cross wavelet transform
my.wc <- analyze.coherency(monthly.soil.moisture, my.pair = c(model.name,"in.situ.imp"),
loess.span = 0,
dt = 1, dj = 1/100,
lowerPeriod = 1/2,
make.pval = TRUE, n.sim = 10)
op <- par(no.readonly = TRUE, font.axis = 2, font.lab = 2) # manipulating labels
# Prints the cross-wavelet plot
wc.image(my.wc, n.levels = 250,
legend.params = list(lab = "cross-wavelet power levels"),
timelab = "", periodlab = "period (months)",
label.time.axis = TRUE,show.date = TRUE, date.format = "%F %T")
### extracting phase angles
# row.closest.to.11 <- which.min((my.wc$Period - 11)^2) # locate period closest to 11
# angle.series <- my.wc$Angle[row.closest.to.11,] # determine angles
# angle.series[1:6] <- NA # omit 6 values at the beginning and 6 values at the end (outside of coi)
# angle.series[(length(angle.series) - 5):length(angle.series)] <- NA
# lead.time.in.minutes <- 30 * (12 * (angle.series / (2*pi))) # convert to minutes: 2*pi corresponds to 12 hours, because period=12
#
# timelag<-mean(lead.time.in.minutes,na.rm = TRUE)
}
getWaveltAnalysys(SSMmonthly.soil.moisture,model.name)## Starting wavelet transformation and coherency computation...
## ... and simulations...
##
|
| | 0%
|
|======= | 10%
|
|============== | 20%
|
|===================== | 30%
|
|============================ | 40%
|
|=================================== | 50%
|
|========================================== | 60%
|
|================================================= | 70%
|
|======================================================== | 80%
|
|=============================================================== | 90%
|
|======================================================================| 100%
## Class attributes are accessible through following names:
## series loess.span dt dj Wave.xy Angle sWave.xy sAngle Power.xy Power.xy.avg Power.xy.pval Power.xy.avg.pval Coherency Coherence Coherence.avg Coherence.pval Coherence.avg.pval Wave.x Wave.y Phase.x Phase.y Ampl.x Ampl.y Power.x Power.y Power.x.avg Power.y.avg Power.x.pval Power.y.pval Power.x.avg.pval Power.y.avg.pval sPower.x sPower.y Ridge.xy Ridge.co Ridge.x Ridge.y Period Scale nc nr coi.1 coi.2 axis.1 axis.2 date.format date.tz
getWaveltAnalysys(RZSmonthly.soil.moisture,model.name)## Starting wavelet transformation and coherency computation...
## ... and simulations...
##
|
| | 0%
|
|======= | 10%
|
|============== | 20%
|
|===================== | 30%
|
|============================ | 40%
|
|=================================== | 50%
|
|========================================== | 60%
|
|================================================= | 70%
|
|======================================================== | 80%
|
|=============================================================== | 90%
|
|======================================================================| 100%
## Class attributes are accessible through following names:
## series loess.span dt dj Wave.xy Angle sWave.xy sAngle Power.xy Power.xy.avg Power.xy.pval Power.xy.avg.pval Coherency Coherence Coherence.avg Coherence.pval Coherence.avg.pval Wave.x Wave.y Phase.x Phase.y Ampl.x Ampl.y Power.x Power.y Power.x.avg Power.y.avg Power.x.pval Power.y.pval Power.x.avg.pval Power.y.avg.pval sPower.x sPower.y Ridge.xy Ridge.co Ridge.x Ridge.y Period Scale nc nr coi.1 coi.2 axis.1 axis.2 date.format date.tz
getWaveltAnalysys(SSMmonthly.soil.moisture,'GLEAM.v3a')## Starting wavelet transformation and coherency computation...
## ... and simulations...
##
|
| | 0%
|
|======= | 10%
|
|============== | 20%
|
|===================== | 30%
|
|============================ | 40%
|
|=================================== | 50%
|
|========================================== | 60%
|
|================================================= | 70%
|
|======================================================== | 80%
|
|=============================================================== | 90%
|
|======================================================================| 100%
## Class attributes are accessible through following names:
## series loess.span dt dj Wave.xy Angle sWave.xy sAngle Power.xy Power.xy.avg Power.xy.pval Power.xy.avg.pval Coherency Coherence Coherence.avg Coherence.pval Coherence.avg.pval Wave.x Wave.y Phase.x Phase.y Ampl.x Ampl.y Power.x Power.y Power.x.avg Power.y.avg Power.x.pval Power.y.pval Power.x.avg.pval Power.y.avg.pval sPower.x sPower.y Ridge.xy Ridge.co Ridge.x Ridge.y Period Scale nc nr coi.1 coi.2 axis.1 axis.2 date.format date.tz
getWaveltAnalysys(RZSmonthly.soil.moisture,"GLEAM.v3a")## Starting wavelet transformation and coherency computation...
## ... and simulations...
##
|
| | 0%
|
|======= | 10%
|
|============== | 20%
|
|===================== | 30%
|
|============================ | 40%
|
|=================================== | 50%
|
|========================================== | 60%
|
|================================================= | 70%
|
|======================================================== | 80%
|
|=============================================================== | 90%
|
|======================================================================| 100%
## Class attributes are accessible through following names:
## series loess.span dt dj Wave.xy Angle sWave.xy sAngle Power.xy Power.xy.avg Power.xy.pval Power.xy.avg.pval Coherency Coherence Coherence.avg Coherence.pval Coherence.avg.pval Wave.x Wave.y Phase.x Phase.y Ampl.x Ampl.y Power.x Power.y Power.x.avg Power.y.avg Power.x.pval Power.y.pval Power.x.avg.pval Power.y.avg.pval sPower.x sPower.y Ridge.xy Ridge.co Ridge.x Ridge.y Period Scale nc nr coi.1 coi.2 axis.1 axis.2 date.format date.tz
The following code computes the boxplots for soil moisture aspects
library(ggplot2)
library(openair)
library(date)
library(lubridate)##
## Attaching package: 'lubridate'## The following object is masked from 'package:base':
##
## datesoil.moisture.subregional <- na.exclude(read.csv("D:\\HESS-2018-546\\Boxplots\\final-SM-aspects-data.csv")) # daily data
get_subregion_BoxAndWhisker_plots<- function(ddayta=soil.moisture.subregional){
#=========================================================================================================
# The function takes daily soil moisture data and return seasonal boxplots
# Arguments: Dataframe
# Inputs: Model name, date, soil moisture values, slope aspect, season and Sub-region
# returns: boxplots for each subregion
#
#=========================================================================================================
ddayta$angle <- factor(ddayta$angle,levels = c("N","E","S","W"))
colnames(ddayta)[4] <- "date"
ddayta$date=as.Date(as.character(ddayta$date))#, format = "%d/%m/%Y")
subdata_DJF <- selectByDate(ddayta, month = c(12,1,2), day = c(1:30))
subdata_DJF$Season <- "DJF"
subdata_JJA <- selectByDate(ddayta, month = c(6,7,8), day = c(1:30))
subdata_JJA$Season <- "JJA"
subdata_MAM <- selectByDate(ddayta, month = c(3,4,5), day = c(1:30))
subdata_MAM$Season <- "MAM"
subdata_SON <- selectByDate(ddayta, month = c(9,10,11), day = c(1:30))
subdata_SON$Season <- "SON"
subdata <- rbind(subdata_JJA,subdata_DJF,subdata_MAM,subdata_SON)
# drop selected columns in the dataframe
#subdata <- subset(subdata,select=-c(X.x,X.1,X.y,X))
subdata$Season <- factor(subdata$Season,levels = c("DJF","MAM","JJA","SON"))
# sort the data by date
Sdata <- subdata[order(as.Date(subdata$date, format="%d/%m/%Y")),]
#write.csv(Sdata, file = "final-soil-moisture-aspect-data.csv")
#Sdata <- read.csv("final-soil-moisture-aspect-data.csv")
Sdata$angle <- factor(Sdata$angle,levels = c("N","E","S","W"))
Sdata$date=as.Date(as.character(Sdata$date))
Sdata$year <- year(Sdata$date)
subdata$area <-as.character(subdata$area)
#introduce intuitive naming for the sub-regions
subdata$area <- replace(subdata$area,subdata$area=="a","OcSa")
subdata$area <- replace(subdata$area,subdata$area=="b","HuSuSa-HoSeSa")
subdata$area <- replace(subdata$area,subdata$area=="c","HoSeSa")
subdata$area <- replace(subdata$area,subdata$area=="d","HoSeGr")
subdata$area <- replace(subdata$area,subdata$area=="e","CoSeGr(I)")
subdata$area <- replace(subdata$area,subdata$area=="f","CoSeGr(II)")
subdata$area <- factor(subdata$area,levels=c("OcSa","HuSuSa-HoSeSa","HoSeSa",
"HoSeGr","CoSeGr(I)","CoSeGr(II)"))
ggplot(subdata,aes(x=Model,y=value,fill=angle))+
geom_boxplot(aes(shape=angle),shape=21,alpha=1,size=0.4,outlier.size=1)+
facet_grid(area~Season,space = "free_y")+
labs(x="Model",y="Soil moisture (%)")+
theme(legend.position="bottom")+
guides(fill=guide_legend(title="Slope aspect"))+
theme(panel.border = element_rect(fill = NA, colour = "black"))+
theme(axis.text.x = element_text(angle = 90,vjust = 0.4,colour = "black"))+
theme(strip.text.y = element_text(angle = 0,colour = "black",size = rel(1)))+
theme(strip.background = element_blank())
}
get_subregion_BoxAndWhisker_plots(soil.moisture.subregional)
The following codes computes mutual information (MI), it has been commented (#) as it takes long to run, uncomment to run the code.
# setwd("D:\\HESS-2018-546\MI\\data-for-computation")
#
# library(raster)
# library(lubridate)
# library(openair)
#
# RootZone.soil.moisture <- list.files(pattern = "*-b.tif")
# Surface.soil.moisture <- list.files(pattern = "*-a.tif")
#
#
# getMI<-function(files){
#
# #************************************************************************************
# # The function computes mutual information from two models namely CCAM-CABLE and GLEAM
# # Arguments: .tif files
# # Inputs: longitude, latitude, and soil moisture
# # returns: raster object with gridded Mutual information calculated per pixel
# #************************************************************************************
#
# # Reads the CCAM-CABLE files
# cordex_25=brick(files[1])
#
# n<-length(files)
# n<-3
# #loop through the GLEAM files
# for (i in 2:n){
#
#
# sm <- brick(files[i])
#
#
# calc.MI <- function(sm, cordex_25){
# # output template
# cor.map <- raster(sm)
# # combine stacks
# T12 <- stack(cordex_25,sm)
#
# rnlayers=nlayers(T12)
# # the function takes a vector, partitions it in half, then correlates
# # the two sections, returning the correlation coefficient.
#
#
#
# calc.eachpixel.MI <- function(myvec,na.rm=T){
# #**************************************************************************
# # Aguments: vector of stacked two model data representing each pixel
# #
# # Returns: Mutual information (MI) between the two models
# #**************************************************************************
#
# # separate the two model data into independent vectors
# myvecT1<-myvec[1:(length(myvec)/2)]
# myvecT2<-myvec[(length(myvec)/2+1):length(myvec)]
#
#
# dates.list <-as.POSIXct(seq.Date(from=as.Date("2011-01-01"), to=as.Date("2014-12-31"), by="day"))
#
#
# #greate data frame with data and two models columns
# SPDATA_DF.Out <- as.data.frame(matrix(nrow=length(dates.list),ncol = 3))
# colnames(SPDATA_DF.Out) <- c("date","cable","gleam")
# SPDATA_DF.Out$date <- dates.list
# SPDATA_DF.Out$cable <- myvecT1
# SPDATA_DF.Out$gleam <- myvecT2
#
# # aggregate daily data to monthly
# SPDATA_DF.Out <- timeAverage(SPDATA_DF.Out, avg.time = "1 month", data.thresh = 80,start.date=as.Date("2011-01-01"))
#
# if((sum((is.na(SPDATA_DF.Out$cable))/length(myvec))*100==0)&&(sum((is.na(SPDATA_DF.Out$gleam))/length(myvec))*100==0)){
#
# #time series detrending functions
# data.ts = ts(SPDATA_DF.Out$cable, frequency=12, start=c(2001,1), end=c(2014,12))
# data.ts2 = ts(SPDATA_DF.Out$gleam, frequency=12, start=c(2001,1), end=c(2014,12))
# fit = stl(data.ts, s.window=7, t.window=11)
# fit2 = stl(data.ts2, s.window=7, t.window=11)
# residual <- fit$time.series[,3]
# residual2 <- fit2$time.series[,3]
#
# #call the mutual information computing function
# res<- mi.data(residual,residual2,discretization.method = "sturges")
#
# }
#
# else{
#
#
# res <- NA
#
# }
#
# return(res)
# }
# #}
# # apply the function above to each cell and write the MI
# # to the output template.
#
# MI.map <- stackApply(T12, indices = rep(1, rnlayers),fun = calc.eachpixel.MI, na.rm = FALSE)
#
# return(MI.map)
#
# }
#
# MI=calc.MI(sm, cordex_25)
#
# writeRaster(MI, filename=paste(files[i],"-MI.tif",sep=""), format="GTiff", overwrite=TRUE,options="INTERLEAVE=BAND")
#
# }
# }
#
# getMI(RootZone.soil.moisture)
# getMI(Surface.soil.moisture)The following codes plot MI resulting from the above code.
library(ggplot2)
library(rgeos)
library(rgdal)
library(raster)
library(scales)
# Mutual information
setwd("D:\\HESS-2018-546\\MI\\data-for-plotting")
Sub.regions <- read.csv("D:\\HESS-2018-546\\Study-region\\selected-subregions.csv") # sub-region data
files <- list.files(pattern = "*-MI.tif")
get.Data.Frame<-function(files){
#=========================================================================================================
# The function takes in raster files of MI, converted them into a dataframe and plot MI
# Arguments: MI raster files
# Inputs: longitude, latitude and MI
# outputs: data frame containing longitude, latitude and MI
#=========================================================================================================
Out.DTF <- as.data.frame(matrix(nrow=0,ncol=5))
names(Out.DTF) <-c("MI","long","lat","Season","Model")
for (i in 1:length(files)) {
coro <-raster(files[i])
Features.file <- unlist(strsplit(files[i],"-"))
model.name <-paste(Features.file[1],Features.file[2],sep="-")
season.name <-Features.file[3]
test_spdf <- as(coro, "SpatialPixelsDataFrame")
test_df <- as.data.frame(test_spdf)
test_df$Season<-season.name
test_df$Model<-model.name
names(test_df) <-c("MI","long","lat","Season","Model")
Out.DTF<- rbind(Out.DTF,test_df)
}
return(Out.DTF)
}
#write.csv(Out.DTF,file = "MI-data-frame-all-models.csv")
#names(Sub.regions) <- c("X" ,"Alt","long","lat" ,"Region")
plot.Maps <-function(files,Sub.regions){#,files2
#===========================================================================================================
# The function plots mutual information computed from daily data
# Arguments: tiff files and Sub-regions data
# Inputs: longitude, latitude, mutual information (MI), model name, soil depth and sub-regions
# Outputs: MI maps with sub-regions
#===========================================================================================================
Out.DTF <-get.Data.Frame(files)
Out.DTF$Season <-factor(as.character(Out.DTF$Season),levels=unique(as.character(Out.DTF$Season)))
Out.DTF$Season <- factor(Out.DTF$Season, levels = c("SMsurf","RZSM"))
region.bound <- readOGR(dsn=path.expand("D:\\HESS-2018-546\\Study-region\\SA_GADMIN"), layer="ZAF_adm1")
region.env <- gEnvelope(region.bound) #calculates rectengular boundry box encompassing your region
region.diff <- gDifference(region.env,region.bound) #extract exterior of your region within the boundrary box
region.diff.df <- fortify(region.diff) #convert your exterior region to a data frame
OVERLAY <- fortify(region.bound, region="Shape_Area") #convert your region shape file to data frame
p <- ggplot()+
geom_raster(data =Out.DTF, mapping = aes(x=long,y=lat,fill=MI),interpolate = F) +
geom_point(data = Sub.regions, aes(x=lon,y=lat,shape=area))+
coord_equal(xlim =c(24,33) ,ylim=c(-29,-22.5))+
facet_grid(Model~Season)+
geom_polygon(data=region.diff.df ,aes(x=long,y=lat,group=group),colour="white",size=1,fill="White")+
labs(fill="MI",x ="longitude",y="latitude")+
scale_fill_gradientn(colours=c("blue","cyan","white", "green","yellow","red"),
values=rescale(c(-2,0-.Machine$double.eps,0,0+.Machine$double.eps,2)))+
guides(fill = guide_colorbar(barwidth = 0.5, barheight = 10))+
theme(strip.background = element_blank())+
theme(panel.border = element_rect(fill = NA, colour = "black"))+
scale_shape_manual(values=c(0,1,2,3,4,5),labels=c("OcSa","HuSuSa-HoSeSa","HoSeSa","HoSeGr"
,"CoSeGr(I)","CoSeGr(II)"))
p
}
plot.Maps(files,Sub.regions) #,files2## OGR data source with driver: ESRI Shapefile
## Source: "D:\HESS-2018-546\Study-region\SA_GADMIN", layer: "ZAF_adm1"
## with 10 features
## It has 16 fields