1 |
f2b2e7d5
|
Adam M. Wilson
|
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 |
|
|
|