library(LatticeKrig)

inpath  = "estimate1/synoptic/"  # note trailing slash!
outpath = "outkrig1/"

var = "trend"

system(paste("mkdir",outpath))

station <- read.table("station.data",
                      header=TRUE,
                      na.strings="-9999.00",
                      stringsAsFactors=FALSE)

files <- dir(inpath)

print(date())

pb = txtProgressBar(min=1,max=length(files),style=2)

for (f in files){

  setTxtProgressBar(pb,which(files==f))
  cat(" ",f)

  infile <- paste0(inpath,f)

  synoptic <- read.table(infile,
                         header=TRUE,
                         na.strings="-99.99",
                         stringsAsFactors=FALSE)

  data = merge(station, synoptic, by="id")

  # discard stations with missing data in variable or elevation
  filter = is.na(data$elev) | is.na(data[,var])  
  data = data[!filter,]

  x = as.matrix(data[,c("lon","lat")])
  y = as.matrix(data[,var])
  z = as.matrix(data[,"elev"])

  LK = LatticeKrig(x=x, y=y, Z=z)

  etazed = predict(LK)
  eta = predict(LK, drop.Z=TRUE)
  zed = etazed - eta

  output = cbind(data$id, zed, etazed, eta)
  colnames(output)=c("id","zed","etazed","eta")
  
  write.table(output,
              file=paste0(outpath,f),
              col.names=TRUE,
              row.names=FALSE,
              quote=FALSE)

}

print(" ")
print(date())
