---
title: "paleo_enkf"
output:
  pdf_document: default
  html_document: default
---

NB this is not turnkey code for every calculation, sensitivity test and graphic in the paper. A single run should reproduce the main result however and various parameters and settings can be edited for sensitivity tests.

This code is provided under GPL-3.0

https://opensource.org/licenses/GPL-3.0

jdannan@blueskiesresearch.org.uk


```{r setup, include=FALSE}

#clear all variables at the start

rm(list = ls())

#choose era from the below, just comment out the one you don't want

#plio stuff will not be included in this release for the paper

era <- "lgm"
#era <- "plio"

knitr::opts_chunk$set(fig.path=paste('figures/',era,'-',sep=""),fig.align='center', fig.show='hold', eval=TRUE, echo=F, fig.height = 6, fig.width = 7,dev=c( 'pdf')) #saving figs if you knit

#ensure reproducibility as far as practicable


set.seed(42)



```


EnKF code for the paleo state estimation.

1. Read in all the models and assemble into a big array

2. Read in the observations and process into an array with the relevant position and error data

3. Do a localised EnKF using perturbed observation technique primarily because I understand it best and it leads to a simple algorithm

Localisation is based on a Gaspari-Cohn distance function with a single length scale over the globe.


Start by setting various parameters:

```{r set_params}
# collect together the parameters to control the analysis

#number of EOFs used to translate the ensemble mean

N_eofs <- 4

#factor for inflating the ensemble
#no inflation used in final version

ens_inflate <- 1. 

#err_scale = scaling factor for the input obs_sst_err and obs_sat_err currently set to 1, i.e. not used in main calculation

sst_err_scale <- 1.
sat_err_scale <- 1.

#lscale = half the cutoff lengthscale for localisation (Gaspari and Cohn)

lscale <- 2500.


```


Now some subroutines

```{r map_plot}
#code to plot a map:
#takes an input array, assumes it's a global field, plots with a fixed red-blue colour scheme based on IPCC AR4 and
#adds mask for SST as appropriate

library(maps)


map_plot<-function(arr,mask,pts=NA,vals=NA,tit,type="lgm"){
  dims<-dim(arr)
  nlon<-dims[1]
  nlat<-dims[2]
  
  lon_long=(0:nlon)*360/nlon-180 #cell edges
  lat_long=(0:nlat)*180/nlat-90 #cell edges again
  
  lon=(1:nlon)*360/nlon-180-180/nlon #cell middles
  lat=(1:nlat)*180/nlat-90-90/nlat #cell middles again

#set the colour ramps for the red and blues and land  
  #these were hand-built to look reasonable
  
  nstep<-11
  nblu<-4
  nblub<-3
  nred<-4
  
  ColorRamp1 <- rgb( seq(0.,0.4,length=nblu),  # Red
                     seq(0.0,0.5,length=nblu),  # Green
                     seq(.9,1.,length=nblu))  # Blue
  ColorRampb <- rgb( seq(0.3,0.8,length=nblub),  # Red
                     seq(0.7,0.9,length=nblub),  # Green
                     seq(1.,1,length=nblub))  # Blue
  
  ColorRamp1[1]<-rgb(0,0,.4)
  
  
  ColorRamp2 <- rgb( seq(1,1,length=nred),  # Red
                     seq(.8,0,length=nred),  # Green
                     seq(.8,0,length=nred))  # Blue
  
  ColorRampland<-rgb(0.6,0.5,0.4)
  
  
  ColorRamp <- c(ColorRamp1,ColorRampb,ColorRamp2)  
  
    
  ColorLevels<-c(-40,-20,-12,-8,-4,-2,-1,0,1,2,4,15)
  
  
  #a different colorbar for the plio which has some neg but mostly warmer
  
  if(type=="plio"){

  ColorLevels<-c(-12,-8,-4,-2,-1,0,1,2,4,8,12,50)
  
 ColorRamp3 <- rgb( seq(.7,.4,length=2),  # Red
                     seq(.0,.0,length=2),  # Green
                     seq(.7,.4,length=2))  # Blue
     
  ColorRamp <- c(tail(ColorRamp,-2),ColorRamp3) 
  }

  #different colourbar, sorry, colorbar for plotting uncertainty maps which are non-negative everywhere
  
  if(type=="unc"){

  ColorLevels<-c(-12,-2,-1,-0.5,-0.2,0,0.2,0.5,1,2,4,50)
  ColorLevels<-c(0,0.2,0.5,1,2,4,50)
  nstep <- length(ColorLevels)-1
  
 ColorRamp3 <- rgb( seq(.7,.4,length=2),  # Red
                     seq(.0,.0,length=2),  # Green
                     seq(.7,.4,length=2))  # Blue
     
  ColorRamp <- c(ColorRamp2,ColorRamp3) 
  nstep <- length(ColorLevels)-1
  }
  
  
    def.par <- par(no.readonly = TRUE) # save default...
 
 par("mar"=c(2,2,6,2))

  par("mgp"=c(2,1,0))

  
  
  layout(matrix(c(1,2), 2, 1, byrow = TRUE),heights=c(3,1))
  
  
  image(lon_long,lat_long,arr,col=ColorRamp,breaks=ColorLevels,axes=F,main=tit,xlab="",ylab="",cex.main=2)
  
  #a hack for the SST plot to add the mask:

if(length(arr)<3000){
    image(lon_long,lat_long,mask,col=rgb(0.6,0.5,0.4),breaks=c(1,2),add=T,axes=F,main=tit,xlab="",ylab="",cex.main=2)
}
  

  tempusr <- par("usr")
  
  par("usr" = c(-180,180,-90,90))
  map(interior=F,add=T,wrap=T,lwd=3)

  par("usr" = tempusr) #make sure we are set up cleanly
  
  axis(1,(-6:6)*30)#,lwd.ticks=-1,labels=NA)
  axis(2,0:6*30-90)#,lwd.ticks=-1,labels=NA)
  axis(3,(-6:6)*30)#,lwd.ticks=-1,labels=NA)
  axis(4,0:6*30-90)#,lwd.ticks=-1,labels=NA)
  
  if(!is.na(pts[1])){
    
    
    cc<-ColorRamp[rowSums(apply(array(ColorLevels,dim=c(length(ColorLevels),length(vals))),1,'<', vals))]
  
    ptsz <- 0.7
  points(pts[,1]-360,pts[,2],col=cc,pch=19,cex=ptsz)
  points(pts[,1],pts[,2],col=cc,pch=19,cex=ptsz)
  points(pts[,1]-360,pts[,2],col="white",lwd=0.7,cex=ptsz*1.1)
  points(pts[,1],pts[,2],col="white",lwd=0.7,cex=ptsz*1.1)
}
  
  par("mar"=c(4,5,1,5))
  
  par("mgp"=c(2.5,1,0))

  
    
  cb<-matrix(ColorLevels[1:(length(ColorLevels)-1)],length(ColorLevels)-1,1)+1e-8


  lab <- "Temperature anomaly"
if(type=="unc")lab <- "Temperature"
    
  image(seq(0,nstep,by=1),c(0,1),cb,breaks=ColorLevels,col=ColorRamp,xlab=lab,ylab="",axes=F,cex.lab=1.5)
  axis(1,0:nstep,tick=T,labels=format(ColorLevels[-c(1,nstep+1)],digits=2),cex.axis=0.9,at=c(1:(nstep-1)))
  axis(1,0:nstep,labels=NA,lwd.ticks=-1)
  axis(2,0:1,lwd.ticks=-1,labels=NA)
  axis(3,0:nstep,lwd.ticks=-1,labels=NA)
  axis(4,0:1,lwd.ticks=-1,labels=NA)

    par(def.par)

    }


```


```{r deg_to_grid_etc}

#a bunch of functions relating to grid manipulation

#function to convert from lat lon coordinates to a regular grid
#just for the 5 and 2 deg grids used here, not a general function

deg_to_grid <- function(lonlat,deg){
  
  #deg is the grid size, either 5 or 2
  
  #lon may be -180 to 180 or 0 to 360, because different people use different conventions...
  
  xval <- round((lonlat[,1] - (lonlat[,1] > 180)*360 + 180 + deg/2)/deg)
  yval <- round((lonlat[,2] + 90 + deg/2)/deg)

   return(cbind(xval,yval)) 
  
  
}

#function to convert from 2D array indices to lon and lat in order

grid_to_deg <- function(coords,deg){

if(deg == 5){  
  
lon <- JT_lon[coords[1]]  
lat <- JT_lat[coords[2]]  
}
if(deg == 2){  
  
lon <- lgm_pmip2_sat_lon[coords[1]]  +360*(lgm_pmip2_sat_lon[coords[1]]<0)
lat <- lgm_pmip2_sat_lat[coords[2]]  
}

  return(cbind(lon,lat))
}

#function to convert from array indices to position in a vectorised form

grid_to_vec <- function(coords,deg){
  
  #deg is the grid size, 5 or 2
  
  if(is.null(dim(coords))){
  
  
  pos <- coords[1] + (coords[2]-1) * (360/deg)

  }else{
    
  
  pos <- coords[,1] + (coords[,2]-1) * (360/deg)
}
   return(pos) 
  
  
}

vec_to_grid <- function(vec){
  
  #aim is to convert a point in the model vector to the grid coords
  #SST comes first in a 5 deg grid and SAT 2nd on a 2 deg grid

  
  SST_pts <- (360/5)*(180/5)
  
  if (vec > SST_pts){
    #SAT point
  vec <- vec - SST_pts
  
  xpos <- (vec-1) %% (360/2) + 1
  ypos <- floor((vec - 1) / (360/2)) + 1 
  
    }
  else{
  xpos <- (vec-1) %% (360/5) + 1
  ypos <- floor((vec - 1) / (360/5)) + 1 
    
    
    
    
  }  
  
 return(cbind(xpos,ypos))
  
}


# Calculates the geodesic distance between two points specified by latitude/longitude using the
# Spherical Law of Cosines (slc)
gcd.deg <- function(lon1, lat1, lon2, lat2) {
  R <- 6371 # Earth mean radius [km]

  int <- min(max(sin(lat1*pi/180)*sin(lat2*pi/180) + cos(lat1*pi/180)*cos(lat2*pi/180) * cos(lon2*pi/180-lon1*pi/180),-1),1)
  d <- acos(int) * R
  if(is.na(d)){
    print(c(lon1,lat1,lon2,lat2))
    print(c(int,d))
    stop()
  }
  return(d) # Distance in km
}
# Convert degrees to radians
deg2rad <- function(deg) return(deg*pi/180)

gc.val <- function(d,lscale){
  #d is the distance and l is the length scale
  r <- d/lscale  

  output <- 0*r #vector of appropriate length
  
  for (i in 1:length(output)){
      
    taper <- 0 
  
  rr <- 2.*r[i]
  
  if (r[i] <= 0.5){
    taper <- ((( -0.25*rr +0.5 )*rr +0.625 )*rr -5.0/3.0 )*rr**2+1.
  } else if ((r[i] > 0.5) & (r[i] < 1.)){
    taper <- (((( rr/12.0 -0.5 )*rr +0.625 )*rr +5.0/3.0 )*rr -5.0 )*rr  + 4.0 - 2.0 / (3.0 * rr)
  }
  
output[i] <- taper
}
return(output)
}


```

```{r ens_stats}
#calculate some summary statistics based on ensemble and observations

analysis <- function(ens,obs_list){
  
  #area weights for the grid points for SAT
  sat_wgt<-t(replicate(length(lgm_pmip2_sat_lon),(cos(lgm_pmip2_sat_lat*pi/180))))
  sat_wgt<-sat_wgt/sum(sat_wgt)

  vec_wgt <- as.vector(sat_wgt)

  sat_trop <- 0*sat_wgt
  sat_trop[,31:60]<-1 #just #30S-30N
  sat_trop_wgt <- sat_trop*sat_wgt
  sat_trop_wgt <- sat_trop_wgt/sum(sat_trop_wgt)
  vec_trop_wgt <- as.vector(sat_trop_wgt)
  
  
    #area weights for the grid points for SST
  sst_wgt<-t(replicate(length(lgm_pmip2_sst_lon),(cos(lgm_pmip2_sst_lat*pi/180))))
  
  sst_wgt <- sst_wgt * (1-mask)
  sst_wgt<-sst_wgt/sum(sst_wgt)
  vec_sst_wgt <- as.vector(sst_wgt)

  
  sst_trop <- 0*sst_wgt
  sst_trop[,13:24]<-1 #just #30S-30N
  sst_trop_wgt <- sst_trop*sst_wgt
  sst_trop_wgt <- sst_trop_wgt/sum(sst_trop_wgt)
  vec_sst_trop_wgt <- as.vector(sst_trop_wgt)
  

 #just re-checking that we don't have NAN values etc though they should have been masked out already 
  
ens_mean <- rowMeans(ens)

ens_mask <- !is.na(ens_mean) | (abs(ens_mean) < 1000)

obs_use <- ens_mask[obs_list$vec]

if(length(obs_use)<length(obs_list$vals)) print("warning extra masking required")

gmeansat <- sum(ens_mean[(prod(dimsst)+1):(prod(dimsst)+prod(dimsat))]*vec_wgt)

print(c("global mean SAT",gmeansat))

gmeansst <- sum(ens_mean[1:(prod(dimsst))]*vec_sst_wgt)

print(c("global mean SST",gmeansst))


num_mods <- dim(ens)[2]

ens_sat <- rep(0,num_mods)
ens_sst <- rep(0,num_mods)

ens_trop_sat <- rep(0,num_mods)
ens_trop_sst <- rep(0,num_mods)

for (i in  1:num_mods) {
  
  ens_sat[i] <- sum(ens[(prod(dimsst)+1):(prod(dimsst)+prod(dimsat)),i]*vec_wgt)
  ens_sst[i] <- sum(ens[(1:(prod(dimsst))),i]*vec_sst_wgt)

  ens_trop_sat[i] <- sum(ens[(prod(dimsst)+1):(prod(dimsst)+prod(dimsat)),i]*vec_trop_wgt)
  ens_trop_sst[i] <- sum(ens[(1:(prod(dimsst))),i]*vec_sst_trop_wgt)
  
}

print(c("ensemble members SAT",ens_sat))
print(c("Mean and SD",mean(ens_sat),sd(ens_sat)))  
print(c("ensemble members SST",ens_sst))
print(c("Mean and SD",mean(ens_sst),sd(ens_sst)))
print(c("ensemble members TROP SAT",ens_trop_sat))
print(c("Mean and SD",mean(ens_trop_sat),sd(ens_trop_sat)))  
print(c("ensemble members TROP SST",ens_trop_sst))
print(c("Mean and SD",mean(ens_trop_sst),sd(ens_trop_sst)))


  

hist(ens_sat,main=paste("Rank histogram for GSAT anomalies"), col = "red")

points(ens_sat,rep(-0.5,length(ens_sat)),pch=1,col="red")


#do a Shapiro-Wilks test on the ens_sat

print("now doing Shapiro-Wilk test of normality")

print(shapiro.test(ens_sat))

#now the obs mismatch

anom <- obs_list$vals[obs_use] - ens_mean[obs_list$vec[obs_use]]


print(c("Mean and RMSE error of ens_mean",mean(anom,na.rm=T),sqrt(mean(anom^2,na.rm=T))))


plot(ens_mean[obs_list$vec[obs_use]],obs_list$vals[obs_use],main = "Scatter plot of fitted estimate vs obs",xlab="Fitted",ylab="Observed (proxy estimate)")
abline(a=0,b=1)

print(c("Number of obs used",length(obs_use)))
print(c("Correlation",cor(ens_mean[obs_list$vec[obs_use]],obs_list$vals[obs_use])))

#now the obs mismatch for each ensemble member in turn

an_mean <- rep(0,dim(ens)[2])
an_rmse <- rep(0,dim(ens)[2])


for (i in 1:(dim(ens)[2])){

anom <- obs_list$vals[obs_use] - ens[obs_list$vec[obs_use],i]


an_mean[i] <- mean(anom,na.rm=T)
an_rmse[i] <- sqrt(mean(anom^2,na.rm=T))


}

print("mean bias")
print(an_mean)
print(c(mean(an_mean), sqrt(mean(an_mean^2))))
print("mean rmse")
print(an_rmse)
print(mean(an_rmse))

}


```

Read in models and observations

```{r readin_lgm_models}
#read in the LGM_PMIP2,3,4 models and assemble into a single array

#all of these have been preprocessed onto 2/5 degree for SAT/SST respectively

#needs a netcdf library

library(ncdf4)


#PMIP2

#load in the models
#LGM_PMIP2
#SAT and grid

file.nc <- nc_open("models/LGM_PMIP2/tas.all.nc")

lgm_pmip2_sat <-ncvar_get(file.nc,"tas")
lgm_pmip2_sat_lon <-ncvar_get(file.nc,"lon")
lgm_pmip2_sat_lat <-ncvar_get(file.nc,"lat")

nc_close(file.nc)

#SST and grid

file.nc <- nc_open("models/LGM_PMIP2/tos.all.nc")

lgm_pmip2_sst <-ncvar_get(file.nc,"tos")
lgm_pmip2_sst_lon <-ncvar_get(file.nc,"lon")
lgm_pmip2_sst_lat <-ncvar_get(file.nc,"lat")

nc_close(file.nc)

#LGM_PMIP3
#SAT and grid

file.nc <- nc_open("models/LGM_PMIP3/tas.all.nc")

lgm_pmip3_sat <-ncvar_get(file.nc,"tas")
lgm_pmip3_sat_lon <-ncvar_get(file.nc,"lon")
lgm_pmip3_sst_lat <-ncvar_get(file.nc,"lat")

nc_close(file.nc)

#SST and grid

file.nc <- nc_open("models/LGM_PMIP3/tos.all.nc")

lgm_pmip3_sst <-ncvar_get(file.nc,"tos")
lgm_pmip3_sst_lon <-ncvar_get(file.nc,"lon")
lgm_pmip3_sst_lat <-ncvar_get(file.nc,"lat")

nc_close(file.nc)


#LGM_PMIP4
#SAT and grid

file.nc <- nc_open("models/LGM_PMIP4_JY/tas.all.nc")

lgm_pmip4_sat <-ncvar_get(file.nc,"tas")
lgm_pmip4_sat_lon <-ncvar_get(file.nc,"lon")
lgm_pmip4_sst_lat <-ncvar_get(file.nc,"lat")

nc_close(file.nc)

#SST and grid

file.nc <- nc_open("models/LGM_PMIP4_JY/tos.all.nc")

lgm_pmip4_sst <-ncvar_get(file.nc,"tos")
lgm_pmip4_sst_lon <-ncvar_get(file.nc,"lon")
lgm_pmip4_sst_lat <-ncvar_get(file.nc,"lat")

nc_close(file.nc)


#set a few values

num_mods_pmip2 <- dim(lgm_pmip2_sat)[3]
num_mods_pmip3 <- dim(lgm_pmip3_sat)[3]
num_mods_pmip4 <- dim(lgm_pmip4_sat)[3]
num_mods_lgm <- num_mods_pmip2 + num_mods_pmip3 + num_mods_pmip4

dim2sst <- dim(lgm_pmip2_sst)
dim2sat <- dim(lgm_pmip2_sat)
vec_pmip2 <- rbind(array(lgm_pmip2_sst,dim=c(dim2sst[1]*dim2sst[2],dim2sst[3])),array(lgm_pmip2_sat,dim=c(dim2sat[1]*dim2sat[2],dim2sat[3])))

dim3sst <- dim(lgm_pmip3_sst)
dim3sat <- dim(lgm_pmip3_sat)
vec_pmip3 <- rbind(array(lgm_pmip3_sst,dim=c(dim3sst[1]*dim3sst[2],dim3sst[3])),array(lgm_pmip3_sat,dim=c(dim3sat[1]*dim3sat[2],dim3sat[3])))

dim4sst <- dim(lgm_pmip4_sst)
dim4sat <- dim(lgm_pmip4_sat)
vec_pmip4 <- rbind(array(lgm_pmip4_sst,dim=c(dim4sst[1]*dim4sst[2],dim4sst[3])),array(lgm_pmip4_sat,dim=c(dim4sat[1]*dim4sat[2],dim4sat[3])))


dimsst <- dim2sst[1:2]
dimsat <- dim2sat[1:2]



### this is the full ensemble as an array made up of stacked vectors of SST followed by SAT

lgm_ens <- cbind(vec_pmip2,vec_pmip3,vec_pmip4)


##### remove some models with huge sea masks as documented in the manuscript.

lgm_ens <- lgm_ens[,c(-10,-11,-15)]

num_mods_lgm <- dim(lgm_ens)[2]

lgm_ens_mean <- rowMeans(lgm_ens)



```

``` {r remove_duplicates_lgm, eval=T}
#check pairwise correlations of the models to see how similar they are
#we will do this for SAT only,because the ocean is closely tied anyway.

#also do the s-w test for normality on the global mean SAT anomaly, before and after, and plot a histogram

if(era == "lgm"){


sat_wgt<-t(replicate(length(lgm_pmip2_sat_lon),(cos(lgm_pmip2_sat_lat*pi/180))))
  sat_wgt<-sat_wgt/sum(sat_wgt)

  vec_wgt <- as.vector(sat_wgt)

#keep the original ensemble    
save_ens <- lgm_ens

sat_ens <- lgm_ens[(dim4sst[1]*dim4sst[2]+1):length(lgm_ens[,1]),]  
sat_ens <- save_ens[(dim4sst[1]*dim4sst[2]+1):length(save_ens[,1]),]  
  
names <- read.table("models/finalmodelnames.txt",skip=1)
pairs <- read.table("models/pairs.txt",skip=1)

#do the histogram and normality test

anom <- colSums(vec_wgt*sat_ens)

hist(anom,main=paste("Histogram for SAT anomalies with all usable PMIP2/3/4 models"), col = "red",breaks=-0.5*(0:25)-0.25)

print(cbind(names,anom))

num_p <- (num_mods_lgm-1)*num_mods_lgm/2

dists <- NULL
dist2 <- NULL
cors <- NULL


for (i in 1:(num_mods_lgm-1)){
  for (j in (i+1):(num_mods_lgm)){
    
        dist2 <- c(dist2,sqrt(mean((sat_ens[,i]-sat_ens[,j])^2)))
    dists <- c(dists,sqrt(sum(vec_wgt*(sat_ens[,i]-sat_ens[,j])^2)))
    cors <- c(cors,cor(sat_ens[,i],sat_ens[,j]))
  
    if(T){
      
#use a distance of 3.2 as threshold for suspicious, this is about 10% of the pairs 
    if(sqrt(sum(vec_wgt*(sat_ens[,i]-sat_ens[,j])^2))<3.2){
      print("suspicious")
      print(c(names[i,2],names[j,2],cor(sat_ens[,i],sat_ens[,j]),sqrt(mean((sat_ens[,i]-sat_ens[,j])^2)),sqrt(sum(vec_wgt*(sat_ens[,i]-sat_ens[,j])^2))))

    
        }
    
    }
  }
  
  
}


#this is the final selection of models to omit due to similarity, hard-wired into the code
  
  omitted <- c(1,7,8,12,15,16,19,21,27) #CESM2-1 (model 19) also omitted due to non-normality of s-w test

  #omitted <- c(1,7,8,12,15,16,21,27) #test with CESM2-1 included to show non-normal distribution via s-w test
  
  
  new_names <- names[-omitted,]
new_sat_ens <- sat_ens[,-omitted]

new_num_mods <- num_mods_lgm - length(omitted)


new_dists <- NULL
new_dist2 <- NULL
new_cors <- NULL

print("now checking similarities again with reduced set of")
print(new_num_mods)


for (i in 1:(new_num_mods-1)){
  for (j in (i+1):(new_num_mods)){
    
    new_dist2 <- c(new_dist2,sqrt(mean((new_sat_ens[,i]-new_sat_ens[,j])^2)))
    new_dists <- c(new_dists,sqrt(sum(vec_wgt*(new_sat_ens[,i]-new_sat_ens[,j])^2)))
    new_cors <- c(new_cors,cor(new_sat_ens[,i],new_sat_ens[,j]))
    
    #use 3.2 threshold for suspicion
    if(sqrt(sum(vec_wgt*(new_sat_ens[,i]-new_sat_ens[,j])^2))<3.2){
      print("suspicious")
      print(c(new_names[i,2],new_names[j,2]))
    print(c(cor(new_sat_ens[,i],new_sat_ens[,j]),i,j))
    print(sqrt(mean((new_sat_ens[,i]-new_sat_ens[,j])^2)))
    print(sqrt(sum(vec_wgt*(new_sat_ens[,i]-new_sat_ens[,j])^2)))
    }
    
    
  }
  
  
}


#now fix the ensemble by omitting models

if(T){

lgm_ens <- lgm_ens[,-omitted]
num_mods <- dim(lgm_ens)[2]
num_mods_lgm <- dim(lgm_ens)[2]
names <- names[-omitted,]

}

#re-make the mask in case it's changed

miss <- is.na(rowSums(lgm_ens)) | (abs(rowSums(lgm_ens)) > 1000)

mask_lgm <- array(miss[1:(dim2sst[1]*dim2sst[2])],dim = c(dim2sst[1],dim2sst[2]))


#now do the histogram for this subset of models

sat_ens <- lgm_ens[(dim4sst[1]*dim4sst[2]+1):length(lgm_ens[,1]),]  

anom_new <- colSums(vec_wgt*sat_ens)

hist(anom_new,main=paste("Histogram for SAT anomalies with all retained PMIP2/3/4 models"), col = "red",breaks=-0.5*(0:25)-0.25)


#now make histogram for publication:


p1 <- hist(anom,breaks=-0.5*(0:25)-0.25)
p2 <- hist(anom_new,breaks=-0.5*(0:25)-0.25)

plot(p1,col="red",xlab="Temperature anomaly",ylab="Frequency",main="Histogram of global mean SAT anomalies of PMIP models")
plot(p2,col="blue",add=T)

text(-10,3,"Red: 9 omitted models")
text(-10,2,"Blue: 19 retained models")

}


```


```{r readin_lgm_obs}


library(ncdf4)

#read in the observations from a few data files and assemble into a single set
#need to reorganise values, locations (deg and grid indices), uncertainties all in straightforward format

#nb the SST data here are the anomaly data sets which we have tested but not used in final version. A later routine loads in the absolute data and makes anomalies...

#these are the original MARGOplus and Bartlein et al data that we used in 2012/13
#formatted as 1 row per point with the data arranged as 
# lat lon val sd mod1 mod2 ... mod9 
#with the models being the PMIP2 values, not used here however in same order

odata<-read.table("data/LGM/MARGOplus_values_2.txt")

MARGO_deg <- cbind(odata[,2],odata[,1])

MARGO_grid <- deg_to_grid(MARGO_deg,5)

#reorder the MARGO data according to the grid we are using, to make it consistent.

MARGO_vec <- grid_to_vec(MARGO_grid,5)

#reordering here???
if(T){
ord <- order(MARGO_vec)

MARGO_grid <- MARGO_grid[ord,]
MARGO_vec <- MARGO_vec[ord]
MARGO_deg <- MARGO_deg[ord,]


MARGO_vals<-odata[ord,3]
MARGO_errs<-odata[ord,4] 
}


#mask out the MARGO according to the ensemble

usable  <- which(!mask_lgm[MARGO_grid])

MARGO_use_vals <- MARGO_vals[usable]
MARGO_use_errs <- MARGO_errs[usable]
MARGO_use_deg <- MARGO_deg[usable,]
MARGO_use_grid <- MARGO_grid[usable,]
MARGO_use_vec <- grid_to_vec(MARGO_use_grid,5)

MARGO_list  <- list(vals=MARGO_use_vals,errs=MARGO_use_errs,deg=MARGO_use_deg,grid=MARGO_use_grid,vec=MARGO_use_vec,num=length(MARGO_use_vals))

#this is the new Tierney data set

#netcdf on a 5x5 grid so needs a little processing to amalgamate with the other data set

file.nc <- nc_open("data/LGM/Tierney2020_ProxyData_5x5_deltaSST.nc")
JT_sst <- ncvar_get(file.nc,"deltaSST")
JT_sd <- ncvar_get(file.nc,"std")
JT_lat <- ncvar_get(file.nc,"lat")
JT_lon <- ncvar_get(file.nc,"lon")
nc_close(file.nc)


JT_lon <- JT_lon + 360*(JT_lon<0) #needs a bit of rewrapping to 0-360 range

#find the data points and extract their position from the array

JT_posarr <- which(!is.nan(JT_sst),arr.ind=T)
JT_vals <- JT_sst[JT_posarr]
JT_errs <- JT_sd[JT_posarr]

#extract the model values corresponding to the JT data locations

JT_mod_vals <- array(0,dim=c(length(JT_vals),dim(lgm_ens)[2]))

JT_deg <- cbind(JT_lon[JT_posarr[,1]],JT_lat[JT_posarr[,2]])

JT_grid <- deg_to_grid(JT_deg,5)

usable  <- which(!mask_lgm[JT_grid])

JT_use_vals <- JT_vals[usable]
JT_use_errs <- JT_errs[usable]
JT_use_deg <- JT_deg[usable,]
JT_use_grid <- JT_grid[usable,]
JT_use_vec <- grid_to_vec(JT_use_grid,5)

JT_list  <- list(vals=JT_use_vals,errs=JT_use_errs,deg=JT_use_deg,grid=JT_use_grid,vec=JT_use_vec,num=length(JT_use_vals))


#now to calculate the joint list

#where is data in MARGO but not JT

extra <- which(!(MARGO_list$vec %in% JT_list$vec))

joint_list <- list(vals=c(JT_list$vals,MARGO_list$vals[extra]),errs=c(JT_list$errs,MARGO_list$errs[extra]),deg=rbind(JT_list$deg,MARGO_list$deg[extra,]),grid=rbind(JT_list$grid,MARGO_list$grid[extra,]),vec=c(JT_list$vec,MARGO_list$vec[extra]),num=length(JT_list$vals)+length(MARGO_list$vals[extra]))




#land SAT obs from Quantrec

ldata<-read.table("data/LGM/Quantrecplus_values_2_exclpts_newdata.txt")

QR_vals <- ldata[,3]
QR_errs <- ldata[,4] *0 +2 #not using stated values but replacing with 2 degrees
QR_lat <- -91+2*(1:90)
QR_lon <- 179+2*(1:180) 
QR_lon <- QR_lon - 360*(QR_lon > 360)

QR_deg <- cbind(ldata[,2],ldata[,1])

QR_grid <- deg_to_grid(QR_deg,2)

QR_vec <- grid_to_vec(QR_grid,2) + 360*180/(5*5)


QR_list <- list(vals=QR_vals,errs=QR_errs,deg=QR_deg,grid=QR_grid,vec=QR_vec,num=length(QR_vals))



#print some analysis of the MARGO and JT data to look at differences etc....

#data in the overlap where I'm going to use JT

overlap_J <- which((JT_list$vec %in% MARGO_list$vec)) 
overlap_M <- which((MARGO_list$vec %in% JT_list$vec))

#number of points and check they are the same locations!

print(c(length(overlap_J),length(overlap_M),MARGO_list$vec[overlap_M]-JT_list$vec[overlap_J]))

print("Mean, sd and RMS difference where they overlap")
print(mean(MARGO_list$vals[overlap_M] - JT_list$vals[overlap_J]))
print(sd(MARGO_list$vals[overlap_M] - JT_list$vals[overlap_J]))
print(sqrt(mean((MARGO_list$vals[overlap_M] - JT_list$vals[overlap_J])^2)))

#various numbers of data points
print(c(MARGO_list$num,JT_list$num,joint_list$num,length(extra)))

```


Show the initial ensemble and the observations:

```{r show_lgm_obs, eval = T}

#leaving this in just to get the right random numbers!

if(era=="lgm"){

#MARGO first
#plotting SST model map with the data on top..

  map_plot(array(lgm_ens_mean[1:(dim2sst[1]*dim2sst[2])],dim=c(dim2sst[1],dim2sst[2])),mask_lgm,pts = MARGO_list$deg,vals = MARGO_list$vals, tit=paste("SST plot and MARGO data"))

model_vals <- lgm_ens[MARGO_list$vec,]

obs_perturb <- array(rnorm(prod(dim(model_vals)),sd=1*MARGO_list$errs),dim=dim(model_vals))

ranks <- rowSums((model_vals+obs_perturb*0) < MARGO_list$vals) + 1

print("ranks for MARGO")
print(ranks)
ranks_margo <- ranks

old.par <- par(mar=c(5,4,4,2)+0.1)

hist(ranks,main=paste("Rank histogram for MARGO ocean data\nwith all usable PMIP2/3/4 models"), col = "red",breaks=(0:(1+dim(lgm_ens)[2]))-0.5 + 1)


#JT now

  map_plot(array(lgm_ens_mean[1:(dim2sst[1]*dim2sst[2])],dim=c(dim2sst[1],dim2sst[2])),mask_lgm,pts = JT_list$deg,vals = JT_list$vals, tit=paste("SST plot and JT data"))

model_vals <- lgm_ens[JT_list$vec,]

obs_perturb <- array(rnorm(prod(dim(model_vals)),sd=1*JT_list$errs),dim=dim(model_vals))

ranks <- rowSums((model_vals+obs_perturb*0) < JT_list$vals) + 1
print("ranks for JT")
print(ranks)
ranks_JT <- ranks

old.par <- par(mar=c(5,4,4,2)+0.1)

hist(ranks,main=paste("Rank histogram for JT ocean data\nwith all usable PMIP2/3/4 models"), col = "red",breaks=(0:(1+dim(lgm_ens)[2]))-0.5 + 1)


#joint now

  map_plot(array(lgm_ens_mean[1:(dim2sst[1]*dim2sst[2])],dim=c(dim2sst[1],dim2sst[2])),mask_lgm,pts = joint_list$deg,vals = joint_list$vals, tit=paste("SST plot and combined MARGO/JT data"))

model_vals <- lgm_ens[joint_list$vec,]

obs_perturb <- array(rnorm(prod(dim(model_vals)),sd=1*joint_list$errs),dim=dim(model_vals))

ranks <- rowSums((model_vals+obs_perturb*0) < joint_list$vals) + 1

print("ranks for joint")
print(ranks)
ranks_joint <- ranks

old.par <- par(mar=c(5,4,4,2)+0.1)

hist(ranks,main=paste("Rank histogram for combined MARGO/JT ocean data\nwith all usable PMIP2/3/4 models"), col = "red",breaks=(0:(1+dim(lgm_ens)[2]))-0.5 + 1)


if(F){
#now rescaling the ensemble by a factor of ....2? ens_inflate
  #no don't do this at this point...we rescale the mean first

ens_mean <- rowMeans(lgm_ens)
ens_scale <- ens_inflate*lgm_ens - (ens_inflate-1)*ens_mean

print(c(sd(ens_scale[1234,]),sd(lgm_ens[1234,])*ens_inflate))

#joint data in inflated ensemble

  map_plot(array(lgm_ens_mean[1:(dim2sst[1]*dim2sst[2])],dim=c(dim2sst[1],dim2sst[2])),mask_lgm,pts = joint_list$deg,vals = joint_list$vals, tit=paste("SST plot and combined MARGO/JT data"))

model_vals <- ens_scale[joint_list$vec,]

obs_perturb <- array(rnorm(prod(dim(model_vals)),sd=1*joint_list$errs),dim=dim(model_vals))

ranks <- rowSums((model_vals+obs_perturb) < joint_list$vals) + 1

old.par <- par(mar=c(5,4,4,2)+0.1)

hist(ranks,main=paste("Rank histogram for combined MARGO/JT ocean data\nwith all inflated PMIP2/3/4 ensemble"), col = "red",breaks=(0:(1+dim(lgm_ens)[2]))-0.5 + 1)
}


map_plot(array(lgm_ens_mean[((dim2sst[1]*dim2sst[2])+1):((dim2sst[1]*dim2sst[2])+prod(dimsat))],dim=c(dimsat[1],dimsat[2])), pts = QR_list$deg,vals = QR_list$vals,tit=paste("SAT plot and QR obs"),type=era)




}

```


```{r make_anomalies}

#these are the SST data I actually use....


#make anomaly data sets for use in LGM reconstruction
#using absolute temps with a common baseline
#baseline is ERSST PI (1854-83)
#three data sets are MARGO, TIERNEY and JOINT (priority to TIERNEY)
#all of these are available in absolute temp form and I've put them on 5x5 grids

file.nc <- nc_open("data/LGM/Tierney.absolute.lgm.nc")
Tierney_sst <- ncvar_get(file.nc,"tos")
Tierney_lat <- ncvar_get(file.nc,"latitude")
Tierney_lon <- ncvar_get(file.nc,"longitude")
nc_close(file.nc)

file.nc <- nc_open("data/LGM/Margo_LGM_tos.nc")
Margo_sst <- ncvar_get(file.nc,"tos")
Margo_lat <- ncvar_get(file.nc,"lat")
Margo_lon <- ncvar_get(file.nc,"lon")
nc_close(file.nc)

#rotate MARGO to be on the same grid!

Margo_sst <- rbind(Margo_sst[37:72,],Margo_sst[1:36,])
Margo_lon <- c(Margo_lon[37:72],Margo_lon[1:36])


#masking the Margo and Tierney stuff a priori in order that all comparisons
#are made on the same basis

if(T){

  Margo_sst[which(mask_lgm)] <- NA
  Tierney_sst[which(mask_lgm)] <- NA

  
}


file.nc <- nc_open("data/modern/ERSST.tos.1854_1883.ann.nc")
Ersst_sst <- ncvar_get(file.nc,"tos")
Ersst_lat <- ncvar_get(file.nc,"lat")
Ersst_lon <- ncvar_get(file.nc,"lon")
nc_close(file.nc)

Ersst_lon <- Ersst_lon + 360*(Ersst_lon <0)


#make anomalies now

Margo_sst_anom <- Margo_sst - Ersst_sst
Tierney_sst_anom <- Tierney_sst - Ersst_sst

Joint_sst_anom <- Margo_sst_anom
Joint_sst_anom[which(!is.na(Tierney_sst_anom))] <- Tierney_sst_anom[which(!is.na(Tierney_sst_anom))]


#process a set of obs into the correct format for using

#

#LGM_sst <- Margo_sst_anom
#LGM_sst <- Tierney_sst_anom
LGM_sst <- Joint_sst_anom
LGM_lon <- Margo_lon
LGM_lat <- Margo_lat


LGM_posarr <- which(!is.na(LGM_sst),arr.ind=T)
LGM_vals <- LGM_sst[LGM_posarr]
LGM_errs <- LGM_vals *0 + 2 #simple coding to make 1 degree


#extract the model values corresponding to the JT data locations

LGM_mod_vals <- array(0,dim=c(length(LGM_vals),dim(lgm_ens)[2]))


LGM_deg <- cbind(LGM_lon[LGM_posarr[,1]],LGM_lat[LGM_posarr[,2]])

LGM_grid <- deg_to_grid(LGM_deg,5)

usable  <- which(!mask_lgm[LGM_grid])

LGM_use_vals <- LGM_vals[usable]
LGM_use_errs <- LGM_errs[usable]
LGM_use_deg <- LGM_deg[usable,]
LGM_use_grid <- LGM_grid[usable,]
LGM_use_vec <- grid_to_vec(LGM_use_grid,5)

LGM_list  <- list(vals=LGM_use_vals,errs=LGM_use_errs,deg=LGM_use_deg,grid=LGM_use_grid,vec=LGM_use_vec,num=length(LGM_use_vals))



```


Set up the analysis:


```{r setup_analysis}

#ens = ensemble in column vector form dim = c(mod_size,mod_num)

if(era=="lgm"){

ens <- lgm_ens #* 0.5 #here we can introduce a scaling factor as a severe test of ensemble bias!

num_mods <- num_mods_lgm

}

if(era=="plio"){
  
  ens  <- pliomip_ens
  num_mods  <- num_mods_plio
  
}


ens_mean <- rowMeans(ens)

mask_vec <- is.na(ens_mean) | (abs(ens_mean) > 1000)
mask  <- array(mask_vec[1:(72*36)],dim=c(72,36))
ens_mean[mask_vec] <- 0
ens[mask_vec,] <- 0

ens_del <- ens-ens_mean
mod_size <- dim(ens)[1]
mod_num <- dim(ens)[2]


#save the initial ensemble for future analysis

ens_initial <- ens


if(era=="lgm"){

#obs_sst_list <- joint_list #original MARGO/Tierney blend

obs_sst_list <- LGM_list #new MARGO/Tierney blend based on their abs vals and anomalised to ERSST, hopefully masked to models???

#obs_sst_list <- joint_list MARGO/Tierney blend based on their anomalies but masked to models

#obs_sst_list <- JT_list #Tierney based on their anomalies

#obs_sst_list <- MARGO_list #AMRGO based on their anomalies



obs_sat_list <- QR_list

#obs_sat_list <- list(vals=NA,errs=NA,deg=NA,grid=NA,vec=NA,num=0)


#package oservations into all_obs for the analysis

all_obs_list <- list(vals=c(obs_sst_list$vals,obs_sat_list$vals),errs=c(obs_sst_list$errs,obs_sat_list$errs),deg=rbind(obs_sst_list$deg,obs_sat_list$deg),grid=rbind(obs_sst_list$grid,obs_sat_list$grid),vec=c(obs_sst_list$vec,obs_sat_list$vec),num=obs_sst_list$num+obs_sat_list$num)

}

#new block for the plio obs

if(era=="plio"){

obs_sst_list <- plio_list  

obs_sat_list <- list(vals=NA,errs=NA,deg=NA,grid=NA,vec=NA,num=0)

all_obs_list <- list(vals=c(obs_sst_list$vals,obs_sat_list$vals),errs=c(obs_sst_list$errs,obs_sat_list$errs),deg=rbind(obs_sst_list$deg,obs_sat_list$deg),grid=rbind(obs_sst_list$grid,obs_sat_list$grid),vec=c(obs_sst_list$vec,obs_sat_list$vec),num=obs_sst_list$num+obs_sat_list$num)

#no sat observations 

all_obs_list <- obs_sst_list

}

print("initial ensemble prior to thinning")
analysis(save_ens,all_obs_list)
print("ensemble after thinning")
analysis(ens_initial,all_obs_list)
analysis(ens,all_obs_list)




```






Recentering

```{r recentre_prior}
#code to shift the prior mean via a variety of methods...using 4 EOFs for main result


eofdecomp <- prcomp(t(ens),center=F)

print("some details of the EOFs")
print(eofdecomp$sdev^2/sum(eofdecomp$sdev^2))
print(cumsum(eofdecomp$sdev^2)/sum(eofdecomp$sdev^2))

print(tail(eofdecomp$sdev^2,-1)/sum(tail(eofdecomp$sdev^2,-1)))
print(cumsum(tail(eofdecomp$sdev^2,-1))/sum(tail(eofdecomp$sdev^2,-1)))


eof_fit <- lm(all_obs_list$vals ~ eofdecomp$rotation[all_obs_list$vec,(1:N_eofs)]+0)


meanfield <- matrix(eofdecomp$rotation[,(1:N_eofs)],nrow=18792) %*% eof_fit$coefficients

fieldfit <- meanfield[all_obs_list$vec]

#this is just a basic check that the code worked

print(c(mean(all_obs_list$vals - fieldfit),mean(eof_fit$residuals)))
print(c(sd(all_obs_list$vals - fieldfit),sd(eof_fit$residuals)))

#this is a check on how the fitted field differs from the ensemble mean
print(c(mean(ens_mean - meanfield),sd(ens_mean - meanfield)))

#save the old ensemble before replacment with the new fitted mean

ens_mean_original  <- ens_mean

ens_original <- ens



#disable this line to omit the translation for a sensitivity test

ens_mean <- meanfield

#above line to put the translation fit at the mean

#use this line to test putting CESM1-2 at the mean of the ensemble, rather than
#using an EOF fit

#ens_mean <- ens[,12]

#above line to put CESM1-2 at the mean

#now we replace the whole ensemble based on the new mean

for (i in 1:num_mods){
  
ens[,i] <-  ens_del[,i]+ens_mean  
  
}

#use this block after the one above to test putting the translated version of CESM1-2 at the mean of the ensemble

#this block now re-uses member 12 (was CESM but may have been translated) as the new mean

if(F){

ens_mean <- ens[,12]

for (i in 1:num_mods){
  
ens[,i] <-  ens_del[,i]+ens_mean  
  
}

}



#now plot the original, new and difference.

 map_plot(array(ens_mean_original[1:(dim2sst[1]*dim2sst[2])],dim=c(dim2sst[1],dim2sst[2])),mask, pts = obs_sst_list$deg,vals = obs_sst_list$vals, tit=paste("Original SST ensemble mean MARGO/Tierney absolute"),type=era)

  map_plot(array(ens_mean[1:(dim2sst[1]*dim2sst[2])],dim=c(dim2sst[1],dim2sst[2])),mask, pts = obs_sst_list$deg,vals = obs_sst_list$vals, tit=paste("Translated SST ensemble mean with data"),type=era)

   map_plot(array(ens_mean[1:(dim2sst[1]*dim2sst[2])]-ens_mean_original[1:(dim2sst[1]*dim2sst[2])],dim=c(dim2sst[1],dim2sst[2])),mask, tit=paste("Difference in means"),type=era)

   
   
   
 map_plot(array(ens_mean_original[(1+dim2sst[1]*dim2sst[2]):((dim2sst[1]*dim2sst[2])+prod(dimsat))],dim=c(dim2sat[1],dim2sat[2])),mask, pts = obs_sat_list$deg,vals = obs_sat_list$vals,tit=paste("Original SAT ensemble mean"),type=era)

  map_plot(array(ens_mean[(1+dim2sst[1]*dim2sst[2]):((dim2sst[1]*dim2sst[2])+prod(dimsat))],dim=c(dim2sat[1],dim2sat[2])),mask,pts = obs_sat_list$deg,vals = obs_sat_list$vals, tit=paste("Translated SAT ensemble mean with data"),type=era)

   map_plot(array(ens_mean[(1+dim2sst[1]*dim2sst[2]):((dim2sst[1]*dim2sst[2])+prod(dimsat))]-ens_mean_original[(1+dim2sst[1]*dim2sst[2]):((dim2sst[1]*dim2sst[2])+prod(dimsat))],dim=c(dim2sat[1],dim2sat[2])),mask, tit=paste("Difference in means"),type=era)


#look quickly at the original and translated ensemble   
   
print("ensemble prior to translation")
  analysis(ens_original,all_obs_list)
print("ensemble after translation")
  analysis(ens,all_obs_list)


```


```{r show_new_rh}


#do rank histograms of original and translated ensemble, without and then with added inflation

#keep this code to maintain prng reproducibility

model_vals <- ens_original[obs_sst_list$vec,]

obs_perturb <- array(rnorm(prod(dim(model_vals)),sd=1*obs_sst_list$errs),dim=dim(model_vals))

ranks <- rowSums((model_vals+obs_perturb*0) < obs_sst_list$vals) + 1

old.par <- par(mar=c(5,4,4,2)+0.1)

hist(ranks,main=paste("Rank histogram for SST data\nwith original ensemble"), col = "red",breaks=(0:(1+dim(ens)[2]))-0.5 + 1)
 

model_vals <- ens_original[all_obs_list$vec,]

obs_perturb <- array(rnorm(prod(dim(model_vals)),sd=1*all_obs_list$errs),dim=dim(model_vals))

ranks <- rowSums((model_vals+obs_perturb*0) < all_obs_list$vals) + 1

old.par <- par(mar=c(5,4,4,2)+0.1)

hist(ranks,main=paste("Rank histogram for all data\nwith original ensemble"), col = "red",breaks=(0:(1+dim(ens)[2]))-0.5 + 1)
 
print(length(ranks))


model_vals <- ens[obs_sst_list$vec,]

obs_perturb <- array(rnorm(prod(dim(model_vals)),sd=1*obs_sst_list$errs),dim=dim(model_vals))

ranks <- rowSums((model_vals+obs_perturb*0) < obs_sst_list$vals) + 1

old.par <- par(mar=c(5,4,4,2)+0.1)

hist(ranks,main=paste("Rank histogram for SST data\nwith translated ensemble"), col = "red",breaks=(0:(1+dim(ens)[2]))-0.5 + 1)



model_vals <- ens[all_obs_list$vec,]

obs_perturb <- array(rnorm(prod(dim(model_vals)),sd=1*all_obs_list$errs),dim=dim(model_vals))

ranks <- rowSums((model_vals+obs_perturb*0) < all_obs_list$vals) + 1

old.par <- par(mar=c(5,4,4,2)+0.1)

hist(ranks,main=paste("Rank histogram for all data\nwith translated ensemble"), col = "red",breaks=(0:(1+dim(ens)[2]))-0.5 + 1)


ranks <- rowSums((model_vals+obs_perturb*1) < all_obs_list$vals) + 1

old.par <- par(mar=c(5,4,4,2)+0.1)

hist(ranks,main=paste("Rank histogram for all data\nwith translated ensemble and perturbed obs"), col = "red",breaks=(0:(1+dim(ens)[2]))-0.5 + 1)
 
print(length(ranks))


 
#now rescaling the ensemble as a test

ens_mean <- rowMeans(ens)
ens_scale <- ens_inflate*ens - (ens_inflate-1)*ens_mean

print(c(sd(ens_scale[1234,]),sd(ens[1234,])*ens_inflate))


model_vals <- ens_scale[obs_sst_list$vec,]

obs_perturb <- array(rnorm(prod(dim(model_vals)),sd=1*obs_sst_list$errs),dim=dim(model_vals))

ranks <- rowSums((model_vals+obs_perturb*0) < obs_sst_list$vals) + 1

old.par <- par(mar=c(5,4,4,2)+0.1)

hist(ranks,main=paste("Rank histogram for SST data\nwith translated and inflated ensemble"), col = "red",breaks=(0:(1+dim(ens)[2]))-0.5 + 1)
 
ranks <- rowSums((model_vals+obs_perturb) < obs_sst_list$vals) + 1

old.par <- par(mar=c(5,4,4,2)+0.1)

hist(ranks,main=paste("Rank histogram for SST data\nwith translated and inflated ensemble with obs errs"), col = "red",breaks=(0:(1+dim(ens)[2]))-0.5 + 1)
 

model_vals <- ens_scale[all_obs_list$vec,]

obs_perturb <- array(rnorm(prod(dim(model_vals)),sd=1*all_obs_list$errs),dim=dim(model_vals))

ranks <- rowSums((model_vals+obs_perturb*0) < all_obs_list$vals) + 1

old.par <- par(mar=c(5,4,4,2)+0.1)

hist(ranks,main=paste("Rank histogram for all proxy-based data\nin prior ensemble, no obs err"), col = "red",breaks=(0:(1+dim(ens)[2]))-0.5 + 1)

model_vals <- ens_scale[all_obs_list$vec,]

obs_perturb <- array(rnorm(prod(dim(model_vals)),sd=1*all_obs_list$errs),dim=dim(model_vals))

ranks <- rowSums((model_vals+obs_perturb*1) < all_obs_list$vals) + 1

old.par <- par(mar=c(5,4,4,2)+0.1)

hist(ranks,main=paste("Rank histogram for all proxy-based data\nin prior ensemble"), col = "red",breaks=(0:(1+dim(ens)[2]))-0.5 + 1)

 
print(length(ranks))


#set the ensemble to be the inflated one for the rest of the analysis
#NB inflation factor set to 1 in default analysis anyway

ens <- ens_scale

ens_del  <- ens - ens_mean

#check the new numbers 

analysis(ens_initial,all_obs_list)
analysis(ens,all_obs_list)

```



```{r do_enkf, eval = T}

#this is the real code here

#actually perform the enkf on the prior and data


#make the perturbations

ens_pert <- ens * 0

if(obs_sst_list$num > 0){
  #do the sst obs perturbations in the case where we have sst obs

  print("making sst perturbations")
  
  for (i in 1:obs_sst_list$num){

    ens_pert[obs_sst_list$vec[i],] <- rnorm(num_mods) * obs_sst_list$errs[i]
    
  
  }  
  
  
}


if(obs_sat_list$num > 0){
  #do the sat obs perturbations
  print("making sat perturbations")

  for (i in 1:obs_sat_list$num){
    
    ens_pert[obs_sat_list$vec[i],] <- rnorm(num_mods) * obs_sat_list$errs[i]
    
  
  }  
  
  
}


ens_pre_enkf <- ens

ens_new <- ens


#loop through positions


ob_used <- rep(0,prod(dimsst))

for (i in 1:(prod(dimsst)+prod(dimsat))){

  #do the analysis for grid point i based on all data within range of it
  
  #calculate location of point i

grid <-   vec_to_grid(i)
lonlat <- grid_to_deg(grid,2+(i<(prod(dimsst)+1))*3)

  ob <- 0

  oblist <- vector()
  ob_used_pos <- vector()
  ob_used_val <- vector()
  decay_list <- vector()
  
  #do SST obs first
  
  if(T){
  for (j in 1:obs_sst_list$num){
  
    #location of observation j
  
    obs_lonlat <- grid_to_deg(obs_sst_list$grid[j,],5)
    

    #check distance between i and j
    
  dist <- gcd.deg(lonlat[1],lonlat[2],obs_lonlat[1],obs_lonlat[2])

  if (dist < 2*lscale){
    

    ob <- ob + 1
    oblist <- c(oblist,j)
    ob_used_pos <- c(ob_used_pos,grid_to_vec(obs_sst_list$grid[j,],5))
    ob_used_val <- c(ob_used_val,obs_sst_list$vals[j])
    decay_list <- c(decay_list,gc.val(dist,2*lscale))
    
  }
  
  
} #end loop through sst positions
    
   
  }
  
    ob_sst_using <- ob

  if(obs_sat_list$num > 0){
  for (j in 1:(obs_sat_list$num)){
  
    #location of observation j
  
    obs_lonlat <- grid_to_deg(obs_sat_list$grid[j,],2)
    

    #check distance between i and j
    
  dist <- gcd.deg(lonlat[1],lonlat[2],obs_lonlat[1],obs_lonlat[2])
  

  if (dist < 2*lscale){
    
    
    ob <- ob + 1
    oblist <- c(oblist,j)
    ob_used_pos <- c(ob_used_pos,(36*72)+grid_to_vec(obs_sat_list$grid[j,],2))
        ob_used_val <- c(ob_used_val,obs_sat_list$vals[j])

    decay_list <- c(decay_list,gc.val(dist,2*lscale))
    
  }
  
  
} #end loop through positions
  }
  ob_sat_using <- ob - ob_sst_using
  
#print(c("Number of obs to use",ob,ob_sst_using,ob_sat_using))
ob_used[i] <- ob
if(ob == 0){
 # print("warning no update here")
 # stop()
  
  } else{
  
#actually do the enkf here

#create the HA' matrix
#small fix here to account for a single observation breaking the matrix notation...

#HA <- ens_del[ob_used_pos,]
HA <- matrix(ens_del[ob_used_pos,],ncol = num_mods)

HPHt <- HA %*% t(HA) / (num_mods-1)

R <- diag(c((obs_sst_list$errs[oblist[1:ob_sst_using]]*sst_err_scale)^2,obs_sat_list$errs[oblist[(ob_sst_using+1):(ob_sst_using+ob_sat_using)]]*sat_err_scale)^2,nrow=ob)
  

inv <- solve(HPHt + R)


PHt_pos <- ens_del[i,] %*% t(HA) / (num_mods - 1)

K_pos <- PHt_pos %*% inv 
 

innovs <- ob_used_val - ens[ob_used_pos,]  + ens_pert[ob_used_pos,]

#distance-based decay of innovations to generate smooth edges
innovs <- innovs * array(decay_list,dim=c(ob,num_mods)) 

ens_new[i,] <- ens_new[i,] + K_pos %*% innovs

}

}



```

```{r plotout, eval=F}


#for (i in 1:mod_num){
#  map_plot(array(ens[1:(dim2sst[1]*dim2sst[2]),i],dim=c(dim2sst[1],dim2sst[2])),mask,pts = JT_pm,vals = JT_vals, tit=paste("SST plot and JT prior",i))
#  map_plot(array(ens_new[1:(dim2sst[1]*dim2sst[2]),i],dim=c(dim2sst[1],dim2sst[2])),mask,pts = JT_pm,vals = JT_vals, tit=paste("SST plot and JT posterior",i))
#  map_plot(array(ens_new[1:(dim2sst[1]*dim2sst[2]),i]-ens[1:(dim2sst[1]*dim2sst[2]),i],dim=c(dim2sst[1],dim2sst[2])),mask,pts = JT_pm,vals = JT_vals, tit=paste("SST plot and JT diff",i))
#}

for (i in 1:num_mods){
  map_plot(array(ens[1:(dim2sst[1]*dim2sst[2]),i],dim=c(dim2sst[1],dim2sst[2])),mask,pts = obs_sst_list$deg,vals = obs_sst_list$vals, tit=paste("SST plot and obs prior",i),type=era)
  map_plot(array(ens_new[1:(dim2sst[1]*dim2sst[2]),i],dim=c(dim2sst[1],dim2sst[2])),mask,pts = obs_sst_list$deg,vals = obs_sst_list$vals, tit=paste("SST plot and obs posterior",i),type=era)
  map_plot(array(ens_new[1:(dim2sst[1]*dim2sst[2]),i]-ens[1:(dim2sst[1]*dim2sst[2]),i],dim=c(dim2sst[1],dim2sst[2])),mask,pts = obs_sst_list$deg,vals = obs_sst_list$vals, tit=paste("SST plot and obs diff",i),type=era)

  map_plot(array(ens[(1+dim2sst[1]*dim2sst[2]):(dim2sst[1]*dim2sst[2]+dim2sat[1]*dim2sat[2]),i],dim=c(dim2sat[1],dim2sat[2])),mask,pts = obs_sat_list$deg,vals = obs_sat_list$vals, tit=paste("SAT plot and obs prior",i),type=era)
  map_plot(array(ens_new[(1+dim2sst[1]*dim2sst[2]):(dim2sst[1]*dim2sst[2]+dim2sat[1]*dim2sat[2]),i],dim=c(dim2sat[1],dim2sat[2])),mask,pts = obs_sat_list$deg,vals = obs_sat_list$vals, tit=paste("SAT plot and obs posterior",i),type=era)
  map_plot(array(ens_new[(1+dim2sst[1]*dim2sst[2]):(dim2sst[1]*dim2sst[2]+dim2sat[1]*dim2sat[2]),i]-ens[(1+dim2sst[1]*dim2sst[2]):(dim2sst[1]*dim2sst[2]+dim2sat[1]*dim2sat[2]),i],dim=c(dim2sat[1],dim2sat[2])),mask,pts = obs_sat_list$deg,vals = obs_sat_list$vals, tit=paste("SAT plot and obs diff",i),type=era)
  
  
  }

```

```{r plot-report}

#now do the means
  map_plot(array(rowMeans(ens[1:(dim2sst[1]*dim2sst[2]),]),dim=c(dim2sst[1],dim2sst[2])),mask,pts = obs_sst_list$deg,vals = obs_sst_list$vals, tit=paste("Prior mean and data, SST"),type=era)
  map_plot(array(rowMeans(ens_new[1:(dim2sst[1]*dim2sst[2]),]),dim=c(dim2sst[1],dim2sst[2])),mask,pts = obs_sst_list$deg,vals = obs_sst_list$vals, tit=paste("Posterior mean and data, SST"),type=era)
  map_plot(array(rowMeans(ens_new[1:(dim2sst[1]*dim2sst[2]),]-ens[1:(dim2sst[1]*dim2sst[2]),]),dim=c(dim2sst[1],dim2sst[2])),mask,pts = obs_sst_list$deg,vals = obs_sst_list$vals, tit=paste("SST Difference"),type=era)

  map_plot(array(rowMeans(ens[(1+dim2sst[1]*dim2sst[2]):(dim2sst[1]*dim2sst[2]+dim2sat[1]*dim2sat[2]),]),dim=c(dim2sat[1],dim2sat[2])),mask,pts = obs_sat_list$deg,vals = obs_sat_list$vals, tit=paste("Prior mean and data, SAT"),type=era)
  map_plot(array(rowMeans(ens_new[(1+dim2sst[1]*dim2sst[2]):(dim2sst[1]*dim2sst[2]+dim2sat[1]*dim2sat[2]),]),dim=c(dim2sat[1],dim2sat[2])),mask,pts = obs_sat_list$deg,vals = obs_sat_list$vals, tit=paste("Posterior mean and data, SAT"),type=era)
  map_plot(array(rowMeans(ens_new[(1+dim2sst[1]*dim2sst[2]):(dim2sst[1]*dim2sst[2]+dim2sat[1]*dim2sat[2]),]-ens[(1+dim2sst[1]*dim2sst[2]):(dim2sst[1]*dim2sst[2]+dim2sat[1]*dim2sat[2]),]),dim=c(dim2sat[1],dim2sat[2])),mask,pts = obs_sat_list$deg,vals = obs_sat_list$vals, tit=paste("SAT difference"),type=era)
  
  
  #now for the uncertainties

essdee <- apply(ens,1,sd)
essdee_new <- apply(ens_new,1,sd)
  
  map_plot(array((essdee[1:(dim2sst[1]*dim2sst[2])]),dim=c(dim2sst[1],dim2sst[2])),mask, tit=paste("Prior uncertainties, SST"),type="unc")
  map_plot(array((essdee_new[1:(dim2sst[1]*dim2sst[2])]),dim=c(dim2sst[1],dim2sst[2])),mask, tit=paste("Posterior uncertainties, SST"),type="unc")
  map_plot(array((essdee_new[1:(dim2sst[1]*dim2sst[2])]-essdee[1:(dim2sst[1]*dim2sst[2])]),dim=c(dim2sst[1],dim2sst[2])),mask, tit=paste("SST uncertainty differences"),type="unc")

  map_plot(array((essdee[(1+dim2sst[1]*dim2sst[2]):(dim2sst[1]*dim2sst[2]+dim2sat[1]*dim2sat[2])]),dim=c(dim2sat[1],dim2sat[2])),mask, tit=paste("Prior uncertainties, SAT"),type="unc")
  map_plot(array((essdee_new[(1+dim2sst[1]*dim2sst[2]):(dim2sst[1]*dim2sst[2]+dim2sat[1]*dim2sat[2])]),dim=c(dim2sat[1],dim2sat[2])),mask, tit=paste("Posterior uncertainties, SAT"),type="unc")
  map_plot(array((essdee_new[(1+dim2sst[1]*dim2sst[2]):(dim2sst[1]*dim2sst[2]+dim2sat[1]*dim2sat[2])]-essdee[(1+dim2sst[1]*dim2sst[2]):(dim2sst[1]*dim2sst[2]+dim2sat[1]*dim2sat[2])]),dim=c(dim2sat[1],dim2sat[2])),mask, tit=paste("SAT uncertainty differences"),type="unc")


model_vals <- ens[all_obs_list$vec,]

#obs_perturb <- array(rnorm(prod(dim(model_vals)),sd=1*all_obs_list$errs),dim=dim(model_vals))

ranks <- rowSums((model_vals+obs_perturb*1) < all_obs_list$vals) + 1

old.par <- par(mar=c(5,4,4,2)+0.1)

hist(ranks,main=paste("Rank histogram for all proxy-based data\nin prior ensemble"), col = "red",breaks=(0:(1+dim(ens)[2]))-0.5 + 1)


model_vals <- ens_new[all_obs_list$vec,]

#obs_perturb <- array(rnorm(prod(dim(model_vals)),sd=1*all_obs_list$errs),dim=dim(model_vals))

ranks <- rowSums((model_vals+obs_perturb*1) < all_obs_list$vals) + 1

old.par <- par(mar=c(5,4,4,2)+0.1)

hist(ranks,main=paste("Rank histogram for all proxy-based data\nin posterior ensemble"), col = "red",breaks=(0:(1+dim(ens)[2]))-0.5 + 1)

 

  
```

```{r stats}
#statistics for the ensemble mean before and after
  
  analysis(ens_initial,all_obs_list)
  analysis(ens,all_obs_list)
  analysis(ens_new,all_obs_list)

```

```{r output_nc}
#make some nc files of the outputs
#just ensemble mean and map of uncertainties for now

dir.create("output")

#SST first

sstdim1=ncdim_def("longitude","degrees",(177.5+(1:72)*5) %% 360)
sstdim2=ncdim_def("latitude","degrees",(-92.5+(1:36)*5))
sstdim3=ncdim_def("model","integer",1:num_mods)

sst <- array(rowMeans(ens_new[1:(dim2sst[1]*dim2sst[2]),]),dim=c(dim2sst[1],dim2sst[2]))
sst[mask] <- NA

vartemp=ncvar_def("SST.mean","Degrees Celsius",list(sstdim1,sstdim2),NA)
nc.ex=nc_create(paste("output/SST.mean.",era,".nc",sep=""),vartemp)
ncvar_put(nc.ex,vartemp,sst)
nc_close(nc.ex)

sstsd <- array(essdee_new[1:(dim2sst[1]*dim2sst[2])],dim=c(dim2sst[1],dim2sst[2]))
sstsd[mask] <- NA

vartemp=ncvar_def("SST.sd","Degrees Celsius",list(sstdim1,sstdim2),NA)
nc.ex=nc_create(paste("output/SST.sd.",era,".nc",sep=""),vartemp)
ncvar_put(nc.ex,vartemp,sstsd)
nc_close(nc.ex)

#SAT

satdim1=ncdim_def("longitude","degrees",(179+(1:180)*2) %% 360)
satdim2=ncdim_def("latitude","degrees",(-91+(1:90)*2))

sat <- array(rowMeans(ens_new[(1+(dim2sst[1]*dim2sst[2])):(dim2sst[1]*dim2sst[2]++dim2sat[1]*dim2sat[2]),]),dim=c(dim2sat[1],dim2sat[2]))

vartemp=ncvar_def("SAT.mean","Degrees Celsius",list(satdim1,satdim2),NA)
nc.ex=nc_create(paste("output/SAT.mean.",era,".nc",sep=""),vartemp)
ncvar_put(nc.ex,vartemp,sat)
nc_close(nc.ex)

satsd <- array(essdee_new[(1+(dim2sst[1]*dim2sst[2])):(dim2sst[1]*dim2sst[2]++dim2sat[1]*dim2sat[2])],dim=c(dim2sat[1],dim2sat[2]))


vartemp=ncvar_def("SAT.sd","Degrees Celsius",list(satdim1,satdim2),NA)
nc.ex=nc_create(paste("output/SAT.sd.",era,".nc",sep=""),vartemp)
ncvar_put(nc.ex,vartemp,satsd)
nc_close(nc.ex)


satall <- array((ens_new[(1+(dim2sst[1]*dim2sst[2])):(dim2sst[1]*dim2sst[2]++dim2sat[1]*dim2sat[2]),]),dim=c(dim2sat[1],dim2sat[2],dim(ens_new)[2]))

vartemp=ncvar_def("SAT.all","Degrees Celsius",list(satdim1,satdim2,sstdim3),NA)
nc.ex=nc_create(paste("output/SAT.all.",era,".nc",sep=""),vartemp)
ncvar_put(nc.ex,vartemp,satall)
nc_close(nc.ex)


sstall <- array((ens_new[1:(dim2sst[1]*dim2sst[2]),]),dim=c(dim2sst[1],dim2sst[2],dim(ens_new)[2]))


for(ii in 1:num_mods){
  ssss <- sstall[,,ii]
  ssss[mask] <- NA
  sstall[,,ii] <- ssss
#sstall[mask,ii] <- NA
}
vartemp=ncvar_def("SST.all","Degrees Celsius",list(sstdim1,sstdim2,sstdim3),NA)
nc.ex=nc_create(paste("output/SST.all.",era,".nc",sep=""),vartemp)
ncvar_put(nc.ex,vartemp,sstall)
nc_close(nc.ex)



```


