# This file imports the data from the station observation file
# and imports it to R. Finally, it puts the data onto an equal 
# are grid.
# Needed:
# - the instrumental data from CRU
# - Martin Tingley's ANOVA code (MatLab)

# The data file looks like this
# 0000300 -1650  -6820 4060 LA.PAZ.EL.ALTO       BOLIVIA       1918 1990
# 6190-9999-9999-9999-9999-9999-9999-9999-9999-9999-9999-9999-9999
# 1918   84   84   94   95   85   74   63   84   93   98  112  108
# 1919  111  118   98  109   98   73   72   75  101  101  107  102
# 1920  100  103   95   96   85   67   62   75   94  107  117  111
# ....
#
# This is a fixed width format with two different formatted lines.

stn.frmt <- c(0, 7, 13, 20, 25, 47, 61, 65, 70)
stn.frmt <- diff( stn.frmt)
dta.frmt <- c(4, rep(5, 12))


Station <- list( number="", coords=c(NA,NA,NA), name="", range=c(NA,NA) )

ReadStationData <- function( filename, offset){
  stn.frmt <- c(0, 7, 13, 20, 25, 47, 61, 65, 70)
  stn.frmt <- diff( stn.frmt)
  stn.clss <- c( "character", "numeric", "numeric", "numeric", "character", "character", "numeric", "numeric")

#  Station <- list( number="", coords=c(NA,NA,NA), name="", range=c(NA,NA) )
  print( offset)
  meta <- read.fwf( filename, skip=offset, n=1, header=F, widths=stn.frmt, colClasses=stn.clss)
  print( meta)
  Station <- list( number=meta[1], coords=unlist(meta[2:3])/100, altitude=meta[4], name=meta[5], range=unlist(meta[7:8]) )
  nyears=diff( Station$range) + 1

  data <- as.matrix( read.fwf( filename, skip=offset+2, header=F, widths=dta.frmt, n=nyears))
  data[which( data == -9999)] <- NA
  data <- data/10

  Station$data <- data

  return(Station)
}

FileName <- "tmp.1609291510.clean.dtb"
con <- file( FileName)
FileLength <- length(readLines(con))
close( con)

StationData <- list()
FileIdx <- 0
while( FileIdx < FileLength){
  StationData[[ length( StationData) + 1]] <- ReadStationData( FileName, FileIdx)
  FileIdx <- FileIdx + 2 + diff( StationData[[ length( StationData)]]$range)+1
}

save( list=c( "StationData"), file="CRU_tmp_inst_aslist.RData")

##### BREAK

load("CRU_tmp_inst_aslist.RData")


nlocs <- length( StationData)

# get the maximum record length. start small and expand for every record
rangemax <- c(1960,1990)
for( tloc in seq( nlocs)){
  rangemax <- range( c(rangemax, StationData[[ tloc]]$range))
}

Inst.Locs <- matrix( NA, nrow=nlocs, ncol=3)
Inst.Time <- seq( rangemax[1], rangemax[2]+11/12, by=1/12)
Inst.Data <- matrix( NA, nrow=length( Inst.Time), ncol=nlocs)

# something changed in the way R handles row and column names...
rownames( Inst.Locs) <- rownames( Inst.Locs, do.NULL = F)
for( tloc in seq( nlocs)){
  Inst.Locs[ tloc,] <- unlist(c( StationData[[ tloc]]$coords[2:1], StationData[[ tloc]]$altitude))
  Inst.Idx <- which( findInterval( Inst.Time, StationData[[ tloc]]$range+c(0,1)) == 1)
  Inst.Data[ Inst.Idx, tloc] <- c( t(StationData[[ tloc]]$data[,-1]) )
  rownames( Inst.Locs)[tloc] <- StationData[[ tloc]]$name
}

require( ggplot2)

Inst.df <- data.frame( x= Inst.Locs[,1], y=Inst.Locs[,2], alt=Inst.Locs[,3], nyears=colSums( !is.na( Inst.Data))/12, first=c(apply( is.na(Inst.Data), 2, which.min)) )

ggplot( Inst.df) + borders() + geom_point( aes(x,y,col=alt))

save( list=c( "Inst.Locs", "Inst.Data", "Inst.Time"), file="CRU_tmp_inst.RData")

# need to turn the Inst.Data into anomalies over a common period - 
# otherwise different altitudes and whatnot will mess up the grid 
# cell means!
load("CRU_tmp_inst.RData")
load("CRU_tmp_inst_aslist.RData")

nlocs <- length( StationData)
Inst.Locs <- matrix( NA, nrow=nlocs, ncol=3)
rownames(Inst.Locs) <- rownames( Inst.Locs, do.NULL=F, prefix="")
Inst.Time <- seq( rangemax[1], rangemax[2], by=1)
Inst.Data <- array( NA, c( time=length( Inst.Time), station=nlocs, month=12))

for( tloc in seq( nlocs)){
  Inst.Locs[ tloc,] <- unlist(c( StationData[[ tloc]]$coords[2:1], StationData[[ tloc]]$altitude))
  Inst.Idx <- which( findInterval( Inst.Time, StationData[[ tloc]]$range+c(0,1)) == 1)
  Inst.Data[ Inst.Idx, tloc, ] <- StationData[[ tloc]]$data[,-1] 
  rownames( Inst.Locs)[tloc] <- trimws(StationData[[ tloc]]$name)
}

# construct summer (JJA) data.
Inst.Data.Summer <- matrix(NA, length( Inst.Time), nlocs)
for( tloc in seq( nlocs)){
  Inst.Data.Summer[, tloc] <- rowMeans( Inst.Data[ , tloc, 6:8])
}
# construct annual data
Inst.Data.Annual <- matrix(NA, length( Inst.Time), nlocs)
for( tloc in seq( nlocs)){
  Inst.Data.Annual[, tloc] <- rowMeans( Inst.Data[ , tloc, ])
}

Inst.Locs.Annual <- Inst.Locs[- which(colSums( !is.na( Inst.Data.Annual)) == 0),]
Inst.Data.Annual <- Inst.Data.Annual[ , - which(colSums( !is.na( Inst.Data.Annual)) == 0)]

# When saving, convert NA to NaN to please MatLAB...

write.table( Inst.Data.Annual[-c(1:100,316),], file="CRU_Annual_1801-2015.csv", row.names=F, col.names=F, sep=",", dec=".")
# cut off the last year... 2016 is incomplete.

Inst.Data.Annual <- Inst.Data.Annual[ -c(1:100,316),]

Inst.Locs.Summer <- Inst.Locs[- which(colSums( !is.na( Inst.Data.Summer)) == 0),]
Inst.Data.Summer <- Inst.Data.Summer[ , - which(colSums( !is.na( Inst.Data.Summer)) == 0)]

write.table( Inst.Data.Summer[-c(1:100),], file="CRU_Summer_1801-2016.csv", row.names=F, col.names=F, sep=",", dec=".")

Inst.Data.Summer <- Inst.Data.Summer[ -(1:100),]

# Now, use Tingley's ANOVA to get the means of all stations. This way we 
# can do real anomalies wrt. the instrumental period.

break()

LocMeans <- read.table( "./LocMeans.out")
# chose below!!
Inst.Data.Summer <- t( t(Inst.Data.Summer) - rowMeans( LocMeans) )
Inst.Data.Summer <- Inst.Data.Summer - mean( Inst.Data.Summer, na.rm=T)

Inst.Data.Annual <- t( t(Inst.Data.Annual) - rowMeans( LocMeans) )
Inst.Data.Inst.Data.Annual <- Inst.Data.Annual - mean( Inst.Data.Annual, na.rm=T)

# Now, the temperatures are in anomalies wrt. 1801-2016 and we can hopefully safely average!

# Now, regrid on an equal area grid.

require( R2G2)
data(grid5000)
Rad2DD = function(lon, lat){
  lonDD = 360 * round((lon / (2 * pi)), 5) - 180 + 0.01
  latDD = 180 * round(lat / pi, 5) - 90 + 0.01
  cbind(lonDD, latDD)
}
grid5000DD = Rad2DD(grid5000[,2], grid5000[, 3])

# Distance function
EarthDistances <- function(xypoint, grid){
  RR<-6378.137; #radius of the earth in km
  Dist <- RR*2*asin( sqrt(sin((xypoint[2] - grid[,2])*pi/180/2)^2 + cos(xypoint[2]*pi/180)*cos(grid[,2]*pi/180)*sin(abs(xypoint[1] - grid[,1])*pi/180/2)^2) )
  return( Dist)
}

GridIdx <- rep(NA, nrow( Inst.Locs.Annual) )
# The next step is really slow.
for( tloc in seq( nrow( Inst.Locs.Annual) )){
  Distances <- EarthDistances(Inst.Locs.Annual[ tloc,1:2], grid5000DD)
  GridIdx[ tloc] <- which.min( Distances)
}

length( unique(GridIdx))
# 1177 locations on the grid that actually have data inside them. Not too much...
save( GridIdx, file="CRU-Grid5000DD-Indices.RData")

OBJ_rMEAN <- function( x){
  # the stuff commented out below adjusts to the n -> oo variance...
  N.Locs <- ncol( x)
  N.Time <- nrow( x)
  H.Patterns <- mgcv::uniquecombs( !is.na(x))
  H.Years <- attr( H.Patterns, "index")
  CorMat <- cor( x, use="pair")

  y <- rowMeans( x, na.rm=T)
  rbar <- rep( NA, nrow( H.Patterns))
  neff <- rep( NA, nrow( H.Patterns))
  r_bar <- y*NA
  n_eff <- y*NA
  for( H.Idx in seq( nrow( H.Patterns)) ){
    Idx <- which( H.Years == H.Idx)
    thispat <- H.Patterns[ H.Idx,]
    nmeas   <- sum( H.Patterns[ H.Idx,] )
    if( nmeas > 1){
      rbar[ H.Idx] <- mean( (CorMat[ thispat==1,thispat==1 ])[lower.tri( diag(nmeas ))])
      neff[ H.Idx] <- nmeas/(1 + (nmeas-1)*rbar[ H.Idx])
      y[ Idx] <- y[ Idx]*sqrt(neff[ H.Idx]) #* sqrt( rbar[ H.Idx])
      r_bar[ Idx] <- rbar[ H.Idx]
      n_eff[ Idx] <- neff[ H.Idx]
    }
  }
  return( y)# cbind(y, r_bar, n_eff))
}


G500.Data <- matrix(NA, nrow=nrow( Inst.Data.Summer), ncol=nrow( grid5000DD))
for( tloc in unique( GridIdx)){
  tIdx <- which( GridIdx == tloc)
  if( length( tIdx) == 1){
    G500.Data[, tloc] <- Inst.Data.Summer[,tIdx]
  } else if ( length( tIdx) > 1){
    G500.Data[, tloc] <- OBJ_rMEAN( Inst.Data.Summer[, which( GridIdx == tloc)])
  }
}

G500.df <- data.frame( x= grid5000DD[,1], y=grid5000DD[,2], nyears=colSums( !is.na( G500.Data)), first=c(apply( is.na(G500.Data), 2, which.min)) )

ggplot( G500.df) + borders() + geom_point( aes(x,y,col=nyears))

# I will actually do the reconstruction with all grid cells first.
ArcIdx <- which( grid5000DD[,2] >= 60)
Inst.DistMat <- EarthDistances( grid5000DD[ ArcIdx,])
plot( Inst.DistMat, cor( G500.Data[,ArcIdx], use="pair") )

Inst.Data <- G500.Data[,ArcIdx]
Inst.Locs <- grid5000DD[ ArcIdx,]

Inst.Time <- 1801:2016

save( list=c("Inst.Data","Inst.Locs","Inst.Time"), file="CRU_InstData_G5000.RData")