#Rcode for ensembleMOS. Michael Angus Dec 2020.
#First install packages and library for ensembleMOS and dealing with Netcdf
install.packages('ensembleMOS')
library(ensembleMOS)
library(ncdf4)
#now, open the matlab created netcdf file. 
ncpath <- "C:/Users/micha/Documents/Work_HEPPI/EMOS/"
ncname_vars <- "map_obs_ens"
ncname_dates <- "dat_vals"
ncfname_vals <- paste(ncpath, ncname_vars, ".nc", sep="")
ncfname_dates <- paste(ncpath, ncname_dates,".nc",sep="")
dname_obs <- "obs_map"
dname_fcst <- "Forecast_ens_map"
dname_enszeros <- "non_zero_ens"
dname_dates <- "dat_vals"
ncin <- nc_open(ncfname_vals)
ncin2 <- nc_open(ncfname_dates)
#load variables
IMD_obs <- ncvar_get(ncin,dname_obs)
fcst <- ncvar_get(ncin,dname_fcst)
ens_zeros <- ncvar_get(ncin,dname_enszeros)
IMD_dates <- ncvar_get(ncin2,dname_dates)
#IMD_obs[IMD_obs==0] <- NA 
fcst[fcst==0] <- NA
#Turn into required data structure. 
#example to follow. 
#data("ensBMAtest",package = "ensembleBMA")
#ensMemNames <- c("gfs","cmcg","eta","gasp","jma","ngps","tcwb","ukmo")
#obs <- paste("T2","obs",sep=".")
#ens <- paste("T2",ensMemNames,sep=".")
#tempTestData <- ensembleData(forecasts=ensBMAtest[,ens],dates = ensBMAtest[,"vdate"],observations = ensBMAtest[,obs],station = ensBMAtest[,"station"],forecastHour=48,initializationTime = "00",exchangeable =c(ens))
#tempTestFit <- ensembleMOS(tempTestData,trainingDays = 25,model="normal")
#Create data object for each station. 
all_points <- dim(fcst)
ens_names <- paste("emp_ens",1:23,sep="")
exGroups<- c(emp_ens1=1,emp_ens2=1,emp_ens3=1,emp_ens4=1,emp_ens5=1,emp_ens6=1,emp_ens7=1,emp_ens8=1,emp_ens9=1,emp_ens10=1,emp_ens11=1,emp_ens12=1,emp_ens13=1,emp_ens14=1,emp_ens15=1,emp_ens16=1,emp_ens17=1,emp_ens18=1,emp_ens19=1,emp_ens20=1,emp_ens21=1,emp_ens22=1,emp_ens23=1)
#initalise variables
a_mat <- 0
B_mat <- 0
c_mat <- 0
d_mat <- 0
q_mat <- 0
for (stat_choice in 1:all_points[1]) {
stat.data <- data.frame(
  emp_ens1 = fcst[stat_choice,1,],
  emp_ens2 = fcst[stat_choice,2,],
  emp_ens3 = fcst[stat_choice,3,],
  emp_ens4 = fcst[stat_choice,4,],
  emp_ens5 = fcst[stat_choice,5,],
  emp_ens6 = fcst[stat_choice,6,],
  emp_ens7 = fcst[stat_choice,7,],
  emp_ens8 = fcst[stat_choice,8,],
  emp_ens9 = fcst[stat_choice,9,],
  emp_ens10 = fcst[stat_choice,10,],
  emp_ens11 = fcst[stat_choice,11,],
  emp_ens12 = fcst[stat_choice,12,],
  emp_ens13 = fcst[stat_choice,13,],
  emp_ens14 = fcst[stat_choice,14,],
  emp_ens15 = fcst[stat_choice,15,],
  emp_ens16 = fcst[stat_choice,16,],
  emp_ens17 = fcst[stat_choice,17,],
  emp_ens18 = fcst[stat_choice,18,],
  emp_ens19 = fcst[stat_choice,19,],
  emp_ens20 = fcst[stat_choice,20,],
  emp_ens21 = fcst[stat_choice,21,],
  emp_ens22 = fcst[stat_choice,22,],
  emp_ens23 = fcst[stat_choice,23,],
  vdates = IMD_dates,
  IMD_obsvals = IMD_obs[stat_choice,]
)
#Calculate data frame
date_ens=ensembleData(forecasts=stat.data[,ens_names],dates=stat.data[,"vdates"],observations=stat.data[,"IMD_obsvals"],forecastHour = 00,initializationTime = "00",exchangeable=exGroups)
#Gamma model
precip_mod_gamma <-ensembleMOS(date_ens, trainingDays = ens_zeros[stat_choice], control = controlMOScsg0(scoringRule = c("crps")),model = "csg0",exchangeable = exGroups)
a_mat[stat_choice]<-precip_mod_gamma[["a"]]
X<-precip_mod_gamma[["B"]]
B_mat[stat_choice]<-X[1]
c_mat[stat_choice]<-precip_mod_gamma[["c"]]
d_mat[stat_choice]<-precip_mod_gamma[["d"]]
q_mat[stat_choice]<-precip_mod_gamma[["q"]]
print(stat_choice)
}
setwd("C:/Users/micha/Documents/Work_HEPPI/EMOS/coefficents")
all_coef_vals=data.frame(a=a_mat,B=B_mat,c=c_mat,d=d_mat,q=q_mat)
write.csv(all_coef_vals,file="IMD_EMOS_coefs.csv")