|
1
|
library(sp)
|
|
2
|
library(spgrass6)
|
|
3
|
library(rgdal)
|
|
4
|
library(reshape)
|
|
5
|
library(ncdf4)
|
|
6
|
library(geosphere)
|
|
7
|
library(rgeos)
|
|
8
|
library(multicore)
|
|
9
|
library(raster)
|
|
10
|
library(lattice)
|
|
11
|
library(rgl)
|
|
12
|
library(hdf5)
|
|
13
|
library(rasterVis)
|
|
14
|
library(heR.Misc)
|
|
15
|
library(spBayes)
|
|
16
|
|
|
17
|
|
|
18
|
X11.options(type="X11")
|
|
19
|
|
|
20
|
ncores=20 #number of threads to use
|
|
21
|
|
|
22
|
|
|
23
|
### copy lulc data to litoria
|
|
24
|
setwd("data/lulc")
|
|
25
|
system("scp atlas:/home/parmentier/data_Oregon_stations/W_Layer* .")
|
|
26
|
|
|
27
|
|
|
28
|
setwd("/home/adamw/acrobates/projects/interp")
|
|
29
|
|
|
30
|
projs=CRS("+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 +datum=NAD83 +units=m +no_defs")
|
|
31
|
months=format(ISOdate(2004,1:12,1),"%b")
|
|
32
|
|
|
33
|
### load station data, subset to stations in region, and transform to sinusoidal
|
|
34
|
load("data/ghcn/roi_ghcn.Rdata")
|
|
35
|
load("data/allstations.Rdata")
|
|
36
|
|
|
37
|
d2=d[d$variable=="tmax"&d$date>=as.Date("2000-01-01"),]
|
|
38
|
d2=d2[,-grep("variable",colnames(d2)),]
|
|
39
|
st2=st[st$id%in%d$id,]
|
|
40
|
st2=spTransform(st2,projs)
|
|
41
|
d2[,c("lon","lat")]=coordinates(st2)[match(d2$id,st2$id),]
|
|
42
|
d2$elev=st2$elev[match(d2$id,st2$id)]
|
|
43
|
d2$month=format(d2$date,"%m")
|
|
44
|
#d2$value=d2$value/10 #convert to mm
|
|
45
|
|
|
46
|
|
|
47
|
## load topographical data
|
|
48
|
topo=brick(as.list(list.files("data/topography",pattern="rst$",full=T)))
|
|
49
|
topo=calc(topo,function(x) ifelse(x<0,NA,x))
|
|
50
|
names(topo)=c("aspect","dem","slope")
|
|
51
|
colnames(topo@data@values)=c("aspect","dem","slope")
|
|
52
|
projection(topo)=projs
|
|
53
|
NS=sin(subset(topo,subset="slope")*pi/180)*cos(subset(topo,subset="aspect")*pi/180);names(NS)="northsouth"
|
|
54
|
EW=sin(subset(topo,subset="slope")*pi/180)*sin(subset(topo,subset="aspect")*pi/180);names(EW)="eastwest"
|
|
55
|
slope=sin(subset(topo,subset="slope")*pi/180);names(slope)="slope"
|
|
56
|
topo2=stack(subset(topo,"dem"),EW,NS,slope)
|
|
57
|
|
|
58
|
## create binned elevation for stratified sampling
|
|
59
|
demc=quantile(subset(topo,subset="dem"))
|
|
60
|
demb=calc(subset(topo,subset="dem"),function(x) as.numeric(cut(x,breaks=demc)))
|
|
61
|
names(demb)="demb"
|
|
62
|
#topo=brick(list(topo,demb))
|
|
63
|
|
|
64
|
|
|
65
|
### load the lulc data as a brick
|
|
66
|
lulc=brick(as.list(list.files("data/lulc",pattern="rst$",full=T)))
|
|
67
|
#projection(lulc)=
|
|
68
|
#plot(lulc)
|
|
69
|
|
|
70
|
## Enter lulc types (from Nakaegawa 2011)
|
|
71
|
lulct=as.data.frame(matrix(c(
|
|
72
|
"Forest",1,
|
|
73
|
"Shrub",2,
|
|
74
|
"Grass",3,
|
|
75
|
"Crop",4,
|
|
76
|
"Mosaic",5,
|
|
77
|
"Urban",6,
|
|
78
|
"Barren",7,
|
|
79
|
"Snow",8,
|
|
80
|
"Wetland",9,
|
|
81
|
"Water",10),byrow=T,ncol=2,dimnames=list(1:10,c("class","id"))),stringsAsFactors=F)
|
|
82
|
colnames(lulc@data@values)=lulct$class[as.numeric(gsub("[a-z]|[A-Z]|[_]|83","",layerNames(lulc)))]
|
|
83
|
layerNames(lulc)=lulct$class[as.numeric(gsub("[a-z]|[A-Z]|[_]|83","",layerNames(lulc)))]
|
|
84
|
lulc=calc(lulc,function(x) ifelse(is.na(x),0,x))
|
|
85
|
projection(lulc)=projs
|
|
86
|
|
|
87
|
### load the LST data
|
|
88
|
lst=brick(as.list(list.files("data/lst",pattern="rescaled.rst$",full=T)[c(4:12,1:3)]))
|
|
89
|
lst=lst-273.15
|
|
90
|
colnames(lst@data@values)=format(as.Date(paste("2000-",as.numeric(gsub("[a-z]|[A-Z]|[_]|83","",layerNames(lst))),"-15",sep="")),"%b")
|
|
91
|
layerNames(lst)=format(as.Date(paste("2000-",as.numeric(gsub("[a-z]|[A-Z]|[_]|83","",layerNames(lst))),"-15",sep="")),"%b")
|
|
92
|
projection(lst)=projs
|
|
93
|
|
|
94
|
|
|
95
|
######################################
|
|
96
|
## compare LULC with station data
|
|
97
|
stlulc=extract(lulc,st2) #overlay stations and LULC values
|
|
98
|
st2$lulc=do.call(c,lapply(apply(stlulc,1,function(x) which.max(x)),function(x) ifelse(is.null(names(x)),NA,names(x))))
|
|
99
|
|
|
100
|
|
|
101
|
### generate sample of points to speed processing
|
|
102
|
n=10000/length(unique(demb))
|
|
103
|
n2=30 #number of knots
|
|
104
|
s=sampleStratified(demb,size=n,sp=T)
|
|
105
|
#s=spsample(as(topo,"SpatialGrid"),n=n,type="regular")
|
|
106
|
s2=spsample(as(topo,"SpatialGrid"),n=n2,type="regular")
|
|
107
|
|
|
108
|
|
|
109
|
s=SpatialPointsDataFrame(s,data=cbind.data.frame(extract(topo,s),extract(topo2,s),extract(lulc,s),extract(lst,s)))
|
|
110
|
### drop areas with no LST data
|
|
111
|
#s=s[!is.na(s$Aug),]
|
|
112
|
s=s[apply(s@data,1,function(x) all(!is.na(x))),]
|
|
113
|
|
|
114
|
### add majority rules lulc for exploration
|
|
115
|
s$lulc=apply(s@data[,colnames(s@data)%in%lulct$class],1,function(x) (colnames(s@data)[colnames(s@data)%in%lulct$class])[which.max(x)])
|
|
116
|
|
|
117
|
|
|
118
|
### spatial regression to fit lulc coefficients
|
|
119
|
|
|
120
|
niter=1250
|
|
121
|
mclapply(months,function(m){
|
|
122
|
print(paste("################################### Starting month ",m))
|
|
123
|
f1=formula(paste(m,"~Shrub+Grass+Crop+Mosaic+Urban+Barren+Snow+Wetland+dem+eastwest+northsouth",sep=""))
|
|
124
|
m1=spLM(f1,data=s@data,coords=coordinates(s),knots=coordinates(s2),
|
|
125
|
starting=list("phi"=0.6,"sigma.sq"=1, "tau.sq"=1),
|
|
126
|
sp.tuning=list("phi"=0.01, "sigma.sq"=0.05, "tau.sq"=0.05),
|
|
127
|
priors=list("phi.Unif"=c(0.3, 3), "sigma.sq.IG"=c(2, 1), "tau.sq.IG"=c(2, 1)),
|
|
128
|
cov.model="exponential",
|
|
129
|
n.samples=niter, verbose=TRUE, n.report=100)
|
|
130
|
assign(paste("mod_",m,sep=""),m1)
|
|
131
|
save(list=paste("mod_",m,sep=""),file=paste("output/mod_",m,".Rdata",sep=""))
|
|
132
|
})
|
|
133
|
|
|
134
|
|
|
135
|
m1.s=mcmc(t(m1$p.samples),start=round(niter/4),thin=1)
|
|
136
|
bwplot(value~X2,melt(as.matrix(m1$p.samples)[,!colnames(m1$p.samples)%in%c("(Intercept)","sigma.sq","tau.sq","phi","northsouth","dem","eastwest")]))
|
|
137
|
densityplot(m1$p.samples)
|
|
138
|
|
|
139
|
### load oregon boundary for comparison
|
|
140
|
roi=spTransform(as(readOGR("data/regions/Test_sites/Oregon.shp","Oregon"),"SpatialLines"),projs)
|
|
141
|
|
|
142
|
|
|
143
|
bgyr=colorRampPalette(c("blue","green","yellow","red"))
|
|
144
|
pdf("output/lst_lulc.pdf",width=11,height=8.5)
|
|
145
|
|
|
146
|
### Summary plots of covariates
|
|
147
|
## LULC
|
|
148
|
at=seq(0.1,100,length=100)
|
|
149
|
levelplot(lulc,at=at,col.regions=bgyr(length(at)),
|
|
150
|
main="Land Cover Classes",sub="Sub-pixel %")+
|
|
151
|
layer(sp.lines(roi, lwd=1.2, col='black'))
|
|
152
|
|
|
153
|
## TOPO
|
|
154
|
at=unique(quantile(as.matrix(subset(topo2,c("eastwest","northsouth","slope"))),seq(0,1,length=100),na.rm=T))
|
|
155
|
levelplot(subset(topo2,c("eastwest","northsouth","slope")),at=at,col.regions=bgyr(length(at)),
|
|
156
|
main="Topographic Variables",sub="")+
|
|
157
|
layer(sp.lines(roi, lwd=1.2, col='black'))
|
|
158
|
|
|
159
|
#LST
|
|
160
|
at=quantile(as.matrix(lst),seq(0,1,length=100),na.rm=T)
|
|
161
|
levelplot(lst,at=at,col.regions=bgyr(length(at)),
|
|
162
|
main="MOD11A1 Mean Monthly LST")+
|
|
163
|
layer(sp.lines(roi, lwd=1.2, col='black'))
|
|
164
|
|
|
165
|
|
|
166
|
### show fitted values
|
|
167
|
spplot(s,zcol="dem")+layer(sp.lines(roi,lwd=1.2,col="black"))+layer(sp.points(s2,col="blue",pch=13,cex=2))
|
|
168
|
|
|
169
|
histogram(~value|variable,data=melt(s@data[,colnames(s@data)%in%lulct$class]))
|
|
170
|
|
|
171
|
bwplot(lulc~value|variable,data=melt(s@data,id.vars="lulc",measure.vars=layerNames(lst)),horizontal=T)
|
|
172
|
|
|
173
|
|
|
174
|
xyplot(value~dem|variable,groups=lulc,melt(s@data[,c("dem","lulc",months)],id.vars=c("dem","lulc")),pch=16,cex=.5,
|
|
175
|
main="Month-by-month scatterplots of Elevation and LST, grouped by LULC",
|
|
176
|
sub="One point per (subsetted) pixel, per month"
|
|
177
|
par.settings = list(superpose.symbol = list(pch=16,cex=.5)),auto.key=list(space="right",title="LULC"))+
|
|
178
|
layer(panel.abline(lm(y~x),col="red"))+
|
|
179
|
layer(panel.text(500,50,bquote(beta==.(round(coefficients(lm(y~x))[2],4)))))
|
|
180
|
|
|
181
|
|
|
182
|
|
|
183
|
|
|
184
|
|
