library(LatticeKrig)

path = "estimate1/synoptic/"  # note trailing slash
var = "trend"

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

files <- dir(path)
f = files[1]

infile <- paste0(path,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)

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


zr = range(c(y, eta, zed, y-eta))

# testing plots
par(mfrow=c(2,3))

quilt.plot(x,y,zlim=zr)
world(add=TRUE)
title("raw data (y)")


quilt.plot(x,eta,zlim=zr)
world(add=TRUE)
title("kriged (eta+zed)")


### Need gridded orography to predict with Z to field

library(ncdf4)

nc = nc_open("../orog.nc")
olat = c(ncvar_get(nc,"lat"))
olon = c(ncvar_get(nc,"lon"))
orog = ncvar_get(nc,"data")

orog = orog[,360:1]
olat = rev(olat)

east = olon < 180
west = olon > 180

orog = rbind(orog[west,],orog[east,])
olon = c(olon[west]-360,olon[east])

xmin = min(station$lon)
xmax = max(station$lon)
ymin = min(station$lat)
ymax = max(station$lat)

regx = olon <= xmax & olon >= xmin
regy = olat <= ymax & olat >= ymin

gridlist = list(x=olon[regx], y=olat[regy])
zo = orog[regx,regy]


image.plot(predictSurface(LK,grid.list=gridlist,Z=zo),zlim=zr)
world(add=TRUE)
title("etazed field")


quilt.plot(x,y-eta,zlim=zr)
world(add=TRUE)
title("residual (y - (eta+zed))")


quilt.plot(x,zed,zlim=zr)
world(add=TRUE)
title("climate (eta)")


image.plot(predictSurface(LK,drop.Z=T),zlim=zr)
world(add=TRUE)
title("eta field")

# dev.copy(png,width=6,height=4,units="in",res=120,file="plot.png")
# dev.off()
