# --- Load necessary R packages ----
# --- Define repository from which R packages will be downloaded.
options(repos=c(CRAN="http://lib.stat.cmu.edu/R/CRAN/"), error=traceback)
# only need ismev package for computing GEV parameters
# and fields package for plotting the spatial maps
# and rgdal to read in the shapefile for A basin
list.of.packages <- c("geoR", "ncdf4", "fields", "lubridate", "MASS")
for (pack in list.of.packages) {
if (!require(pack, character.only = TRUE)) {
install.packages(pack)
require(pack, character.only = TRUE)
}
}
rm(list.of.packages, pack)
setwd("/Users/andrew/Desktop/workshop-BA/5_metamodel")
PP <- read.table("../1_setup/data/Salado_prcp.dat",header=T)
yr <- year(PP[,1])
mo <- month(PP[,1])
da <- day(PP[,1])
uyr <- unique(yr)
nyr <- length(uyr)
rm(PP)
dd.in <- nc_open("../1_setup/data/2015-11-11-BigPampa-observed-interpolation-1961-2013.nc")
dd.out <- nc_open("../1_setup/data/mike_output.nc")
dd.et <- nc_open("../1_setup/data/mike_output_et.nc")
# setwd("../")
prcp <- ncvar_get(dd.in, "prcp")
tmax <- ncvar_get(dd.in, "tmax")
tmin <- ncvar_get(dd.in, "tmin")
srad <- ncvar_get(dd.in, "srad")
et0 <- ncvar_get(dd.in, "et0")
napa <- ncvar_get(dd.out, "wdepth")
elev <- ncvar_get(dd.out, "elevation")
et <- ncvar_get(dd.et, "wdepth")
lon.GK <- ncvar_get(dd.out, "x_coord")/1000
lat.GK <- ncvar_get(dd.out, "y_coord")/1000
nc_close(dd.in)
nc_close(dd.out)
nc_close(dd.et)
predloc <- read.table("../1_setup/data/grids.txt",header=T)
lon <- matrix(predloc[,3],nrow=dim(napa)[2],ncol=dim(napa)[3],byrow=T)
lat <- matrix(predloc[,4],nrow=dim(napa)[2],ncol=dim(napa)[3],byrow=T)
# number of days in data set
tt <- nrow(napa)
# 12 months per year, nyr years in data set
tt.month <- 12*nyr
yr.month <- rep(uyr,each=12)
mo.month <- rep(1:12,nyr)
nx <- dim(napa)[2]
ny <- dim(napa)[3]
xx.id <- seq(from=0,to=nx,by=2)
yy.id <- seq(from=0,to=ny,by=2)
boundary.mask <- basin.mask <- matrix(data=NA,nrow=dim(napa)[2],ncol=dim(napa)[3])
for(xx in 1:nx){
for(yy in 1:ny){
if(is.na(mean(napa[2:tt,xx,yy])) | (mean(napa[2:tt,xx,yy]) == -1)){
if(!is.na(mean(napa[2:tt,xx,yy])) & mean(napa[2:tt,xx,yy]) == -1){
boundary.mask[xx,yy] = -1
basin.mask[xx,yy] = 1
}else{
next
}
}
basin.mask[xx,yy] = 1
}
}
new.napa <- array(data=NA,dim=dim(napa)[c(1:3)])
for(k in 1:tt){
print(k)
look <- image.smooth(napa[k,,], dx=5, dy=5, theta=5)
new.napa[k,,] <- look$z
new.napa[k,,][is.na(basin.mask)] = NA
new.napa[k,,][boundary.mask==-1] = NA
}
napa.m.s <- array(data=NA,dim=c(tt.month,dim(new.napa)[2:3]))
for(k in 1:tt.month){
print(k)
month.id <- yr==yr.month[k] & mo==mo.month[k]
napa.m.s[k,,] <- apply(new.napa[month.id,,],2:3,mean)
}
prcp.m <- tmax.m <- tmin.m <- srad.m <- et0.m <- et.m <- napa.m <-
array(data=NA,dim=c(tt.month,dim(napa)[2:3]))
for(k in 1:tt.month){
print(k)
month.id <- yr==yr.month[k] & mo==mo.month[k]
prcp.m[k,,] <- apply(prcp[month.id,,],c(2:3),mean)
tmax.m[k,,] <- apply(tmax[month.id,,],c(2:3),mean)
tmin.m[k,,] <- apply(tmin[month.id,,],c(2:3),mean)
srad.m[k,,] <- apply(srad[month.id,,],c(2:3),mean)
et0.m[k,,] <- apply(et0[month.id,,],c(2:3),mean)
et.m[k,,] <- apply(et[month.id,,],c(2:3),mean)
napa.m[k,,] <- apply(napa[month.id,,],c(2:3),mean)
}
girasol.y <- readRDS("../1_setup/data/landUse-girasol.rds")
maize.y <- readRDS("../1_setup/data/landUse-maize.rds")
pasturas.y <- readRDS("../1_setup/data/landUse-pasturas.rds")
soja1.y <- readRDS("../1_setup/data/landUse-soja1.rds")
trigo.y <- readRDS("../1_setup/data/landUse-trigo.rds")
trigosoja2.y <- readRDS("../1_setup/data/landUse-trigosoja2.rds")
girasol.m <- maize.m <- pasturas.m <- soja1.m <- trigo.m <- trigosoja2.m <- array(data=NA,dim=c(tt.month,dim(napa)[2:3]))
for(k in 1:tt.month){
print(k)
index.id <- (yr.month[k] - 1960)
girasol.m[k,,] <- t(girasol.y[index.id,,])
maize.m[k,,] <- t(maize.y[index.id,,])
pasturas.m[k,,] <- t(pasturas.y[index.id,,])
soja1.m[k,,] <- t(soja1.y[index.id,,])
trigo.m[k,,] <- t(trigo.y[index.id,,])
trigosoja2.m[k,,] <- t(trigosoja2.y[index.id,,])
}
save.image("metamodel-setup.RData")
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.