1
|
#setup
|
2
|
mo=9
|
3
|
day=2
|
4
|
year=2010
|
5
|
datelabel=format(ISOdate(year,mo,day),"%b %d, %Y")
|
6
|
#
|
7
|
# Step 1 - 10 year monthly averages
|
8
|
#
|
9
|
library(raster)
|
10
|
old<-getwd()
|
11
|
setwd("c:/data/benoit/data_Oregon_stations_Brian_04242012")
|
12
|
l=list.files(pattern="mean_month.*rescaled.tif")
|
13
|
molst<-stack(l)
|
14
|
setwd(old)
|
15
|
molst=molst-273.16 #K->C
|
16
|
idx <- seq(as.Date('2010-01-15'), as.Date('2010-12-15'), 'month')
|
17
|
molst <- setZ(molst, idx)
|
18
|
layerNames(molst) <- month.abb
|
19
|
themolst<-raster(molst,mo)
|
20
|
plot(themolst)
|
21
|
#
|
22
|
# Step 2 - Weather station means across same days
|
23
|
#
|
24
|
# ??? which years & what quality flags???
|
25
|
#select ghcn.id, lat,lon, elev, month, avg(value/10) as "TMax", count(*) as "NumDays" from ghcn, stations where ghcn.id in (select id from stations where state=='OR') and ghcn.id==stations.id and value<>-9999 and year>=2000 and element=='TMAX' group by stations.id, month;select ghcn.id, lat,lon, elev, month, avg(value/10) as "TMax", count(*) as "NumDays" from ghcn, stations where ghcn.id in (select id from stations where state=='OR') and ghcn.id==stations.id and value<>-9999 and year>=2000 and element=='TMAX' group by stations.id, month;
|
26
|
#below table from above SQL query
|
27
|
dst<-read.csv('/data/benoit/data_oregon_stations_brian_04242012/station_means.csv',h=T)
|
28
|
modst=dst[dst$month==mo,]
|
29
|
#
|
30
|
# Step 3 - get LST at stations
|
31
|
#
|
32
|
sta_lola=modst[,c("lon","lat")]
|
33
|
library(rgdal)
|
34
|
proj_str="+proj=lcc +lat_1=43 +lat_2=45.5 +lat_0=41.75 +lon_0=-120.5 +x_0=400000 +y_0=0 +ellps=GRS80 +units=m +no_defs";
|
35
|
lookup<-function(r,lat,lon) {
|
36
|
xy<-project(cbind(lon,lat),proj_str);
|
37
|
cidx<-cellFromXY(r,xy);
|
38
|
return(r[cidx])
|
39
|
}
|
40
|
sta_tmax_from_lst=lookup(themolst,sta_lola$lat,sta_lola$lon)
|
41
|
#
|
42
|
# Step 4 - bias at stations
|
43
|
#
|
44
|
sta_bias=sta_tmax_from_lst-modst$TMax;
|
45
|
bias_xy=project(as.matrix(sta_lola),proj_str)
|
46
|
# windows()
|
47
|
plot(modst$TMax,sta_tmax_from_lst,xlab="Station mo Tmax",ylab="LST mo Tmax")
|
48
|
abline(0,1)
|
49
|
#
|
50
|
# Step 5 - interpolate bias
|
51
|
#
|
52
|
# ?? include covariates like elev, distance to coast, cloud frequency, tree height
|
53
|
library(fields)
|
54
|
#windows()
|
55
|
quilt.plot(sta_lola,sta_bias,main="Bias at stations",asp=1)
|
56
|
US(add=T,col="magenta",lwd=2)
|
57
|
#fitbias<-Tps(bias_xy,sta_bias) #use TPS or krige
|
58
|
fitbias<-Krig(bias_xy,sta_bias,theta=1e5) #use TPS or krige
|
59
|
# windows()
|
60
|
surface(fitbias,col=rev(terrain.colors(100)),asp=1,main="Interpolated bias")
|
61
|
#US(add=T,col="magenta",lwd=2)
|
62
|
#
|
63
|
# Step 6 - return to daily station data & calcualate delta=daily T-monthly T from stations
|
64
|
#
|
65
|
library(RSQLite)
|
66
|
m<-dbDriver("SQLite")
|
67
|
con<-dbConnect(m,dbname='c:/data/ghcn_tmintmax.db')
|
68
|
querystr=paste("select ghcn.id, value as 'dailyTmax' from ghcn where ghcn.id in (select id from stations where state=='OR') and value<>-9999",
|
69
|
"and year==",year,"and month==",mo,"and day==",day,
|
70
|
"and element=='TMAX' ")
|
71
|
rs<-dbSendQuery(con,querystr)
|
72
|
d<-fetch(rs,n=-1)
|
73
|
dbClearResult(rs)
|
74
|
dbDisconnect(con)
|
75
|
d$dailyTmax=d$dailyTmax/10 #stored as 1/10 degree C to allow integer storage
|
76
|
dmoday=merge(modst,d,by="id")
|
77
|
# windows()
|
78
|
plot(dailyTmax~TMax,data=dmoday,xlab="Mo Tmax",ylab=paste("Daily for",datelabel),main="across stations in OR")
|
79
|
#
|
80
|
# Step 7 - interpolate delta across space
|
81
|
#
|
82
|
daily_sta_lola=dmoday[,c("lon","lat")] #could be same as before but why assume merge does this - assume not
|
83
|
daily_sta_xy=project(as.matrix(daily_sta_lola),proj_str)
|
84
|
daily_delta=dmoday$dailyTmax-dmoday$TMax
|
85
|
#windows()
|
86
|
quilt.plot(daily_sta_lola,daily_delta,asp=1,main="Station delta for Jan 15")
|
87
|
US(add=T,col="magenta",lwd=2)
|
88
|
#fitdelta<-Tps(daily_sta_xy,daily_delta) #use TPS or krige
|
89
|
fitdelta<-Krig(daily_sta_xy,daily_delta,theta=1e5) #use TPS or krige
|
90
|
# windows()
|
91
|
surface(fitdelta,col=rev(terrain.colors(100)),asp=1,main="Interpolated delta")
|
92
|
#US(add=T,col="magenta",lwd=2)
|
93
|
#
|
94
|
# Step 8 - assemble final answer - T=LST+Bias(interpolated)+delta(interpolated)
|
95
|
#
|
96
|
bias_rast=interpolate(themolst,fitbias)
|
97
|
plot(bias_rast,main="Raster bias")
|
98
|
daily_delta_rast=interpolate(themolst,fitdelta)
|
99
|
plot(daily_delta_rast,main="Raster Daily Delta")
|
100
|
tmax_predicted=themolst+daily_delta_rast-bias_rast
|
101
|
plot(tmax_predicted,main="Predicted daily")
|
102
|
#
|
103
|
# check
|
104
|
#
|
105
|
sta_pred=lookup(tmax_predicted,daily_sta_lola$lat,daily_sta_lola$lon)
|
106
|
RMSE<-function(x,y) {return(mean((x-y)^2)^0.5)}
|
107
|
rmse=RMSE(sta_pred,dmoday$dailyTmax)
|
108
|
plot(sta_pred~dmoday$dailyTmax,xlab=paste("Actual daily for",datelabel),ylab="Pred daily",main=paste("RMSE=",rmse))
|
109
|
abline(0,1)
|
110
|
resid=sta_pred-dmoday$dailyTmax
|
111
|
quilt.plot(daily_sta_lola,resid)
|