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 |
c7192b5a
|
Adam M. Wilson
|
library(xtable)
|
17 |
|
|
library(ellipse) # for correlation matrix
|
18 |
|
|
library(maptools) # for rgshhs
|
19 |
f2b2e7d5
|
Adam M. Wilson
|
|
20 |
|
|
X11.options(type="X11")
|
21 |
|
|
|
22 |
|
|
ncores=20 #number of threads to use
|
23 |
|
|
|
24 |
|
|
|
25 |
|
|
### copy lulc data to litoria
|
26 |
|
|
setwd("data/lulc")
|
27 |
|
|
system("scp atlas:/home/parmentier/data_Oregon_stations/W_Layer* .")
|
28 |
|
|
|
29 |
|
|
|
30 |
|
|
setwd("/home/adamw/acrobates/projects/interp")
|
31 |
|
|
|
32 |
|
|
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")
|
33 |
|
|
months=format(ISOdate(2004,1:12,1),"%b")
|
34 |
|
|
|
35 |
|
|
### load station data, subset to stations in region, and transform to sinusoidal
|
36 |
|
|
load("data/ghcn/roi_ghcn.Rdata")
|
37 |
|
|
load("data/allstations.Rdata")
|
38 |
|
|
|
39 |
|
|
d2=d[d$variable=="tmax"&d$date>=as.Date("2000-01-01"),]
|
40 |
|
|
d2=d2[,-grep("variable",colnames(d2)),]
|
41 |
|
|
st2=st[st$id%in%d$id,]
|
42 |
|
|
st2=spTransform(st2,projs)
|
43 |
|
|
d2[,c("lon","lat")]=coordinates(st2)[match(d2$id,st2$id),]
|
44 |
|
|
d2$elev=st2$elev[match(d2$id,st2$id)]
|
45 |
|
|
d2$month=format(d2$date,"%m")
|
46 |
|
|
#d2$value=d2$value/10 #convert to mm
|
47 |
|
|
|
48 |
|
|
|
49 |
|
|
## load topographical data
|
50 |
|
|
topo=brick(as.list(list.files("data/topography",pattern="rst$",full=T)))
|
51 |
|
|
topo=calc(topo,function(x) ifelse(x<0,NA,x))
|
52 |
|
|
names(topo)=c("aspect","dem","slope")
|
53 |
|
|
colnames(topo@data@values)=c("aspect","dem","slope")
|
54 |
|
|
projection(topo)=projs
|
55 |
|
|
NS=sin(subset(topo,subset="slope")*pi/180)*cos(subset(topo,subset="aspect")*pi/180);names(NS)="northsouth"
|
56 |
|
|
EW=sin(subset(topo,subset="slope")*pi/180)*sin(subset(topo,subset="aspect")*pi/180);names(EW)="eastwest"
|
57 |
|
|
slope=sin(subset(topo,subset="slope")*pi/180);names(slope)="slope"
|
58 |
c7192b5a
|
Adam M. Wilson
|
|
59 |
|
|
## load coastline data
|
60 |
|
|
roill=bbox(projectExtent(topo,CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")))
|
61 |
|
|
coast=getRgshhsMap("/home/adamw/acrobates/Global/GSHHS_shp/gshhs/gshhs_f.b",xlim=c(roill[1,1],roill[1,2]),ylim=c(roill[2,1],roill[2,2]))
|
62 |
|
|
coast=as(as(coast,"SpatialLines"),"SpatialLinesDataFrame")
|
63 |
|
|
coast@data=data.frame(id=1:length(coast@lines)) #convert to dataframe
|
64 |
|
|
coast=spTransform(coast,projs)
|
65 |
|
|
rcoast=rasterize(coast,topo)
|
66 |
|
|
dist=distance(rcoast)/1000 #get distance to coast and convert to km
|
67 |
|
|
### transform distance to coast using b-log(dist), where b \approx log(DTC of the farthest point on earth from the sea)
|
68 |
|
|
### Garcia-Castellanos, Daniel, and Umberto Lombardo. 2007. “Poles of Inaccessibility: A Calculation Algorithm for the Remotest Places on Earth.” Scottish Geographical Journal 123 (3): 227–233. doi:10.1080/14702540801897809.
|
69 |
|
|
### farthest point is ~2510km from ocean
|
70 |
|
|
log(2510) #round up to 8 to be sure we won't get negative numbers
|
71 |
|
|
dist2=brick(dist,calc(dist,function(x) (8-log(x+1))))
|
72 |
|
|
layerNames(dist2)=c("DTCkm","logDTCkm")
|
73 |
|
|
colnames(dist2@data@values)=c("DTCkm","logDTCkm")
|
74 |
|
|
topo2=stack(subset(topo,"dem"),EW,NS,slope,dist2)
|
75 |
f2b2e7d5
|
Adam M. Wilson
|
|
76 |
|
|
## create binned elevation for stratified sampling
|
77 |
|
|
demc=quantile(subset(topo,subset="dem"))
|
78 |
|
|
demb=calc(subset(topo,subset="dem"),function(x) as.numeric(cut(x,breaks=demc)))
|
79 |
|
|
names(demb)="demb"
|
80 |
c7192b5a
|
Adam M. Wilson
|
##topo=brick(list(topo,demb))
|
81 |
f2b2e7d5
|
Adam M. Wilson
|
|
82 |
|
|
|
83 |
|
|
### load the lulc data as a brick
|
84 |
|
|
lulc=brick(as.list(list.files("data/lulc",pattern="rst$",full=T)))
|
85 |
|
|
#projection(lulc)=
|
86 |
|
|
#plot(lulc)
|
87 |
|
|
|
88 |
|
|
## Enter lulc types (from Nakaegawa 2011)
|
89 |
|
|
lulct=as.data.frame(matrix(c(
|
90 |
|
|
"Forest",1,
|
91 |
|
|
"Shrub",2,
|
92 |
|
|
"Grass",3,
|
93 |
|
|
"Crop",4,
|
94 |
|
|
"Mosaic",5,
|
95 |
|
|
"Urban",6,
|
96 |
|
|
"Barren",7,
|
97 |
|
|
"Snow",8,
|
98 |
|
|
"Wetland",9,
|
99 |
|
|
"Water",10),byrow=T,ncol=2,dimnames=list(1:10,c("class","id"))),stringsAsFactors=F)
|
100 |
|
|
colnames(lulc@data@values)=lulct$class[as.numeric(gsub("[a-z]|[A-Z]|[_]|83","",layerNames(lulc)))]
|
101 |
|
|
layerNames(lulc)=lulct$class[as.numeric(gsub("[a-z]|[A-Z]|[_]|83","",layerNames(lulc)))]
|
102 |
|
|
lulc=calc(lulc,function(x) ifelse(is.na(x),0,x))
|
103 |
|
|
projection(lulc)=projs
|
104 |
|
|
|
105 |
|
|
### load the LST data
|
106 |
|
|
lst=brick(as.list(list.files("data/lst",pattern="rescaled.rst$",full=T)[c(4:12,1:3)]))
|
107 |
|
|
lst=lst-273.15
|
108 |
|
|
colnames(lst@data@values)=format(as.Date(paste("2000-",as.numeric(gsub("[a-z]|[A-Z]|[_]|83","",layerNames(lst))),"-15",sep="")),"%b")
|
109 |
|
|
layerNames(lst)=format(as.Date(paste("2000-",as.numeric(gsub("[a-z]|[A-Z]|[_]|83","",layerNames(lst))),"-15",sep="")),"%b")
|
110 |
|
|
projection(lst)=projs
|
111 |
|
|
|
112 |
|
|
|
113 |
|
|
######################################
|
114 |
|
|
## compare LULC with station data
|
115 |
c7192b5a
|
Adam M. Wilson
|
st2=SpatialPointsDataFrame(st2,data=cbind.data.frame(st2@data,demb=extract(demb,st2),extract(topo,st2),extract(topo2,st2),extract(lulc,st2),extract(lst,st2)))
|
116 |
f2b2e7d5
|
Adam M. Wilson
|
stlulc=extract(lulc,st2) #overlay stations and LULC values
|
117 |
|
|
st2$lulc=do.call(c,lapply(apply(stlulc,1,function(x) which.max(x)),function(x) ifelse(is.null(names(x)),NA,names(x))))
|
118 |
|
|
|
119 |
|
|
|
120 |
|
|
### generate sample of points to speed processing
|
121 |
|
|
n=10000/length(unique(demb))
|
122 |
|
|
n2=30 #number of knots
|
123 |
|
|
s=sampleStratified(demb,size=n,sp=T)
|
124 |
|
|
#s=spsample(as(topo,"SpatialGrid"),n=n,type="regular")
|
125 |
|
|
s2=spsample(as(topo,"SpatialGrid"),n=n2,type="regular")
|
126 |
|
|
|
127 |
|
|
|
128 |
c7192b5a
|
Adam M. Wilson
|
s=SpatialPointsDataFrame(s,data=cbind.data.frame(x=coordinates(s)[,1],y=coordinates(s)[,2],demb=extract(demb,s),extract(topo,s),extract(topo2,s),extract(lulc,s),extract(lst,s)))
|
129 |
f2b2e7d5
|
Adam M. Wilson
|
### drop areas with no LST data
|
130 |
|
|
#s=s[!is.na(s$Aug),]
|
131 |
|
|
s=s[apply(s@data,1,function(x) all(!is.na(x))),]
|
132 |
|
|
|
133 |
|
|
### add majority rules lulc for exploration
|
134 |
|
|
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)])
|
135 |
c7192b5a
|
Adam M. Wilson
|
save(s,file=paste("output/mod_sample.Rdata",sep=""))
|
136 |
f2b2e7d5
|
Adam M. Wilson
|
|
137 |
|
|
|
138 |
|
|
### spatial regression to fit lulc coefficients
|
139 |
c7192b5a
|
Adam M. Wilson
|
Sys.setenv(MKL_NUM_THREADS=24)
|
140 |
|
|
system("export MKL_NUM_THREADS=24")
|
141 |
f2b2e7d5
|
Adam M. Wilson
|
|
142 |
c7192b5a
|
Adam M. Wilson
|
niter=12500
|
143 |
|
|
lapply(months,function(m){
|
144 |
f2b2e7d5
|
Adam M. Wilson
|
print(paste("################################### Starting month ",m))
|
145 |
c7192b5a
|
Adam M. Wilson
|
f1=formula(paste(m,"~logDTCkm+y+Shrub+Grass+Crop+Mosaic+Urban+Barren+Snow+Wetland+dem+eastwest+northsouth",sep=""))
|
146 |
f2b2e7d5
|
Adam M. Wilson
|
m1=spLM(f1,data=s@data,coords=coordinates(s),knots=coordinates(s2),
|
147 |
|
|
starting=list("phi"=0.6,"sigma.sq"=1, "tau.sq"=1),
|
148 |
|
|
sp.tuning=list("phi"=0.01, "sigma.sq"=0.05, "tau.sq"=0.05),
|
149 |
|
|
priors=list("phi.Unif"=c(0.3, 3), "sigma.sq.IG"=c(2, 1), "tau.sq.IG"=c(2, 1)),
|
150 |
|
|
cov.model="exponential",
|
151 |
c7192b5a
|
Adam M. Wilson
|
n.samples=niter,sub.samples=c(2500,niter,10),verbose=TRUE, n.report=100)
|
152 |
f2b2e7d5
|
Adam M. Wilson
|
assign(paste("mod_",m,sep=""),m1)
|
153 |
|
|
save(list=paste("mod_",m,sep=""),file=paste("output/mod_",m,".Rdata",sep=""))
|
154 |
|
|
})
|
155 |
|
|
|
156 |
c7192b5a
|
Adam M. Wilson
|
### Read in results
|
157 |
|
|
load(paste("output/mod_sample.Rdata",sep=""))
|
158 |
|
|
ms=lapply(months,function(m) {print(m); load(paste("output/mod_",m,".Rdata",sep="")) ; get(paste("mod_",m,sep="")) })
|
159 |
|
|
|
160 |
|
|
### generate summaries
|
161 |
|
|
ms1=do.call(rbind.data.frame,lapply(1:length(ms),function(i){
|
162 |
|
|
mi=ms[[i]]
|
163 |
|
|
m1.s=mcmc(t(mi$p.samples),start=round(niter/4),thin=1)
|
164 |
|
|
m1.s2=as.data.frame(t(apply(m1.s,1,quantile,c(0.025,0.5,0.975))))
|
165 |
|
|
colnames(m1.s2)=c("Q2.5","Q50","Q97.5")
|
166 |
|
|
m1.s2$parm=rownames(m1.s)
|
167 |
|
|
m1.s2$month=factor(months[i],levels=months,ordered=T)
|
168 |
|
|
m1.s2$type=ifelse(m1.s2$parm%in%c("tau.sq","sigma.sq","phi"),"Spatial",ifelse(m1.s2$parm%in%c("(Intercept)","eastwest","northsouth","dem","logDTCkm","y"),"Topography","LULC"))
|
169 |
|
|
return(m1.s2)
|
170 |
|
|
}))
|
171 |
|
|
## drop spatial parameters
|
172 |
|
|
#ms1=ms1[!ms1$type%in%"Spatial",]
|
173 |
|
|
## Convert dem to degrees/km to make it comparable to other topographical parameters
|
174 |
|
|
#ms1[ms1$parm%in%c("x","y","dem"),c("Q2.5","Q50","Q97.5")]=ms1[ms1$parm%in%c("x","y","dem"),c("Q2.5","Q50","Q97.5")]
|
175 |
f2b2e7d5
|
Adam M. Wilson
|
|
176 |
|
|
|
177 |
|
|
### load oregon boundary for comparison
|
178 |
|
|
roi=spTransform(as(readOGR("data/regions/Test_sites/Oregon.shp","Oregon"),"SpatialLines"),projs)
|
179 |
|
|
|
180 |
|
|
|
181 |
|
|
bgyr=colorRampPalette(c("blue","green","yellow","red"))
|
182 |
|
|
pdf("output/lst_lulc.pdf",width=11,height=8.5)
|
183 |
|
|
|
184 |
|
|
### Summary plots of covariates
|
185 |
|
|
## LULC
|
186 |
|
|
at=seq(0.1,100,length=100)
|
187 |
|
|
levelplot(lulc,at=at,col.regions=bgyr(length(at)),
|
188 |
|
|
main="Land Cover Classes",sub="Sub-pixel %")+
|
189 |
|
|
layer(sp.lines(roi, lwd=1.2, col='black'))
|
190 |
|
|
|
191 |
|
|
## TOPO
|
192 |
|
|
at=unique(quantile(as.matrix(subset(topo2,c("eastwest","northsouth","slope"))),seq(0,1,length=100),na.rm=T))
|
193 |
|
|
levelplot(subset(topo2,c("eastwest","northsouth","slope")),at=at,col.regions=bgyr(length(at)),
|
194 |
|
|
main="Topographic Variables",sub="")+
|
195 |
|
|
layer(sp.lines(roi, lwd=1.2, col='black'))
|
196 |
|
|
|
197 |
c7192b5a
|
Adam M. Wilson
|
|
198 |
|
|
## DEM
|
199 |
|
|
at=unique(quantile(as.matrix(subset(topo2,c("dem"))),seq(0,1,length=100),na.rm=T))
|
200 |
|
|
levelplot(subset(topo2,c("dem")),at=at,col.regions=bgyr(length(at)),
|
201 |
|
|
main="Elevation",sub="",margin=F)+
|
202 |
|
|
layer(sp.lines(roi, lwd=1.2, col='black'))
|
203 |
|
|
|
204 |
|
|
## DTC
|
205 |
|
|
at=unique(quantile(as.matrix(subset(topo2,c("logDTCkm"))),seq(0,1,length=100),na.rm=T))
|
206 |
|
|
levelplot(subset(topo2,c("logDTCkm")),at=at,col.regions=bgyr(length(at)),
|
207 |
|
|
main="Elevation",sub="",margin=F)+
|
208 |
|
|
layer(sp.lines(roi, lwd=1.2, col='black'))
|
209 |
|
|
|
210 |
f2b2e7d5
|
Adam M. Wilson
|
#LST
|
211 |
|
|
at=quantile(as.matrix(lst),seq(0,1,length=100),na.rm=T)
|
212 |
|
|
levelplot(lst,at=at,col.regions=bgyr(length(at)),
|
213 |
|
|
main="MOD11A1 Mean Monthly LST")+
|
214 |
|
|
layer(sp.lines(roi, lwd=1.2, col='black'))
|
215 |
|
|
|
216 |
|
|
|
217 |
|
|
xyplot(value~dem|variable,groups=lulc,melt(s@data[,c("dem","lulc",months)],id.vars=c("dem","lulc")),pch=16,cex=.5,
|
218 |
|
|
main="Month-by-month scatterplots of Elevation and LST, grouped by LULC",
|
219 |
c7192b5a
|
Adam M. Wilson
|
sub="One point per pixel, per month",ylab="LST",
|
220 |
f2b2e7d5
|
Adam M. Wilson
|
par.settings = list(superpose.symbol = list(pch=16,cex=.5)),auto.key=list(space="right",title="LULC"))+
|
221 |
|
|
layer(panel.abline(lm(y~x),col="red"))+
|
222 |
|
|
layer(panel.text(500,50,bquote(beta==.(round(coefficients(lm(y~x))[2],4)))))
|
223 |
|
|
|
224 |
c7192b5a
|
Adam M. Wilson
|
## just forest, grouped by distance to coast
|
225 |
|
|
xyplot(value~dem|variable,groups=cut(DTCkm,breaks=c(0,300,1000)),melt(s@data[s$lulc=="Forest",c("DTCkm","dem","lulc",months)],id.vars=c("DTCkm","dem","lulc")),pch=16,cex=.5,
|
226 |
|
|
main="Month-by-month scatterplots of Elevation and LST, grouped by Distance to Coast",
|
227 |
|
|
sub="Showing only forest class",ylab="LST",
|
228 |
|
|
par.settings = list(superpose.symbol = list(pch=16,cex=.5)),auto.key=list(space="right",title="Distance to \n Coast (km)"))+
|
229 |
|
|
layer(panel.abline(lm(y~x),col="red"))+
|
230 |
|
|
layer(panel.text(500,50,bquote(beta==.(round(coefficients(lm(y~x))[2],4)))))
|
231 |
f2b2e7d5
|
Adam M. Wilson
|
|
232 |
c7192b5a
|
Adam M. Wilson
|
useOuterStrips(xyplot(value~dem|variable+lulc,groups=cut(DTCkm,breaks=c(0,300,1000)),
|
233 |
|
|
melt(s@data[s$lulc!="Snow",c("DTCkm","dem","lulc",months)],id.vars=c("DTCkm","dem","lulc")),pch=16,cex=.5,
|
234 |
|
|
main="Month-by-month scatterplots of Elevation and LST, grouped by LULC and colored by distance to coast",
|
235 |
|
|
sub="One point per pixel, per month",ylab="LST",
|
236 |
|
|
par.settings = list(superpose.symbol = list(pch=16,cex=.5)),auto.key=list(space="right",title="Distance to \n Coast (km)")))+
|
237 |
|
|
layer(panel.abline(lm(y~x),col="red"))
|
238 |
f2b2e7d5
|
Adam M. Wilson
|
|
239 |
|
|
|
240 |
|
|
|
241 |
c7192b5a
|
Adam M. Wilson
|
##########################################
|
242 |
|
|
### Model output
|
243 |
|
|
|
244 |
|
|
### show fitted values
|
245 |
|
|
spplot(s,zcol="demb",col.regions=terrain.colors(5),cex=.5)+layer(sp.lines(roi,lwd=1.2,col="black"))+
|
246 |
|
|
layer(sp.points(s2,col="blue",pch=13,cex=2,lwd=2))+
|
247 |
|
|
layer(sp.points(st2,col="black",lwd=1.5))
|
248 |
|
|
|
249 |
|
|
|
250 |
|
|
## look at correlation of LULC variables
|
251 |
|
|
lulct=s@data[,colnames(s@data)[colnames(s@data)%in%unique(s$lulc)]]
|
252 |
|
|
lulcc=cor(lulct)[order(cor(lulct)[1,]),order(cor(lulct)[1,])]
|
253 |
|
|
round(cor(lulct),2)
|
254 |
|
|
#print(xtable(round(cor(lulct),2),include.row.names=F))
|
255 |
|
|
plotcorr(lulcc,col=colors <- c("#A50F15","#DE2D26","#FB6A4A","#FCAE91","#FEE5D9","white","#EFF3FF","#BDD7E7","#6BAED6","#3182BD","#08519C")[5*lulcc+6],
|
256 |
|
|
main="Correlation Matrix")
|
257 |
|
|
|
258 |
|
|
t1=as.data.frame(table(s$lulc));colnames(t1)=c("class","samplefreq")
|
259 |
|
|
t1$sampleprop=t1$samplefreq/sum(t1$samplefreq)
|
260 |
|
|
t2=as.data.frame(table(st2$lulc))
|
261 |
|
|
t1$stationfreq=0
|
262 |
|
|
t1$stationfreq[t1$class%in%t2$Var1]=t2$Freq[match(t2$Var1,t1$class[t1$class%in%t2$Var1])]
|
263 |
|
|
t1$stationprop=t1$stationfreq/sum(t1$stationfreq)
|
264 |
|
|
print(xtable(t1[,c(1,3,5)]),include.rownames=F)
|
265 |
|
|
|
266 |
|
|
### two plots of model parameters by month to show the effects of LULC and topography on LST
|
267 |
|
|
## Spatial
|
268 |
|
|
xyplot(Q50~month|parm,data=ms1[ms1$type=="Spatial",],panel=function(x,y,subscripts){
|
269 |
|
|
tms1=ms1[ms1$type=="Spatial",][subscripts,]
|
270 |
|
|
panel.segments(tms1$month,tms1$Q2.5,tms1$month,tms1$Q97.5,groups=tms1$parm,subscripts=1:nrow(tms1))
|
271 |
|
|
panel.xyplot(tms1$month,tms1$Q50,pch=16,type="l")
|
272 |
|
|
panel.abline(h=0,lty="dashed",col="grey")
|
273 |
|
|
},auto.key=list(space="right"),par.settings = list(superpose.symbol = list(pch=16,cex=.5,col=1:15),superpose.line=list(col=1:15)),
|
274 |
|
|
subscripts=T,scales=list(y=list(relation="free",rot=90)),ylab="Parameter Coefficient (slope) with 95% Credible Intervals",xlab="Month",
|
275 |
|
|
main="Spatial Parameters",
|
276 |
|
|
sub="")
|
277 |
|
|
|
278 |
|
|
## LULC
|
279 |
|
|
xyplot(Q50~month,groups=parm,data=ms1[ms1$type=="LULC",],panel=function(x,y,subscripts){
|
280 |
|
|
tms1=ms1[ms1$type=="LULC",][subscripts,]
|
281 |
|
|
sig=ifelse(tms1$Q2.5<0&tms1$Q97.5<0|tms1$Q2.5>0&tms1$Q97.5>0,"red","black")
|
282 |
|
|
panel.segments(tms1$month,tms1$Q2.5,tms1$month,tms1$Q97.5,groups=tms1$parm,subscripts=1:nrow(tms1))
|
283 |
|
|
panel.xyplot(tms1$month,tms1$Q50,groups=tms1$parm,pch=16,type="l",subscripts=1:nrow(tms1))
|
284 |
|
|
panel.abline(h=0,lty="dashed",col="grey")
|
285 |
|
|
},auto.key=list(space="right"),par.settings = list(superpose.symbol = list(pch=16,cex=.5,col=1:15),superpose.line=list(col=1:15)),
|
286 |
|
|
subscripts=T,scales=list(y=list(relation="free",rot=90)),ylab="Parameter Coefficient (slope) with 95% Credible Intervals",xlab="Month",
|
287 |
|
|
main="Effects of LULC on LST",
|
288 |
|
|
sub="Coefficients are unstandardized and represent the change in LST expected with a 1% increase in that class from 100% Forest")
|
289 |
|
|
|
290 |
|
|
## Topography
|
291 |
|
|
xyplot(Q50~month|parm,data=ms1[ms1$type=="Topography",],panel=function(x,y,subscripts){
|
292 |
|
|
tms1=ms1[ms1$type=="Topography",][subscripts,]
|
293 |
|
|
panel.segments(tms1$month,tms1$Q2.5,tms1$month,tms1$Q97.5,subscripts=1:nrow(tms1))
|
294 |
|
|
panel.xyplot(tms1$month,tms1$Q50,groups=tms1$parm,pch=16,type="l",subscripts=1:nrow(tms1))
|
295 |
|
|
panel.abline(h=0,lty="dashed",col="grey")
|
296 |
|
|
},auto.key=list(space="right"),par.settings = list(superpose.symbol = list(pch=16,cex=.5,col=1:15),superpose.line=list(col=1:15)),
|
297 |
|
|
subscripts=T,scales=list(y=list(relation="free",rot=90)),ylab="Parameter Coefficient (slope) with 95% Credible Intervals",xlab="Month",
|
298 |
|
|
main="Effects of Topography on LST",
|
299 |
|
|
sub="Coefficients are unstandardized. Intercept is degrees C, eastwest and northsouth range \n from -1 (90 degree slope) to 0 (flat) to 1 (90 degree slope), and dem is m \n logDTCkm is in log(km), and y (lat) is m")
|
300 |
|
|
|
301 |
|
|
### Capture same pattern using only station data?
|
302 |
|
|
|
303 |
|
|
|
304 |
|
|
dev.off()
|
305 |
|
|
|
306 |
|
|
shrinkpdf<-function(pdf,maxsize=1,suffix="_small",verbose=T){
|
307 |
|
|
require(multicore)
|
308 |
|
|
wd=getwd()
|
309 |
|
|
td=paste(tempdir(),"/pdf",sep="")
|
310 |
|
|
if(!file.exists(td)) dir.create(td)
|
311 |
|
|
if(verbose) print("Performing initial compression")
|
312 |
|
|
setwd(td)
|
313 |
|
|
if(verbose) print("Splitting the PDF to parallelize the processing")
|
314 |
|
|
system(paste("pdftk ",wd,"/",pdf," burst",sep=""))
|
315 |
|
|
mclapply(list.files(td,pattern="pdf$"),function(f){
|
316 |
|
|
## loop through all pages, perform compression with with ps2pdf
|
317 |
|
|
if(verbose) print("Performing initial compression")
|
318 |
|
|
system(paste("ps2pdf -dUseFlateCompression=true ",td,"/",f," ",td,"/compressed_",f,sep=""))
|
319 |
|
|
file.rename(paste(td,"/compressed_",f,sep=""),paste(td,"/",f,sep=""))
|
320 |
|
|
## get sysmte size
|
321 |
|
|
size=file.info(paste(td,"/",f,sep=""))$size*0.000001 #get sizes of individual pages
|
322 |
|
|
toobig=size>=maxsize
|
323 |
|
|
if(verbose&toobig) print(paste("Resizing ",basename(paste(td,"/",f,sep="")),sep=""))
|
324 |
|
|
system(paste("gs -dBATCH -dTextAlphaBits=4 -dNOPAUSE -r300 -q -sDEVICE=png16m -sOutputFile=- -q ",paste(td,"/",f,sep="")," | convert -background transparent -quality 100 -density 300 - ",f,sep=""))
|
325 |
|
|
})
|
326 |
|
|
if(verbose) print("Compiling the final pdf")
|
327 |
|
|
setwd(wd)
|
328 |
|
|
system(paste("gs -dBATCH -dNOPAUSE -q -sDEVICE=pdfwrite -sOutputFile=",strsplit(pdf,".",fixed=T)[[1]][1],suffix,".pdf ",td,"/*.pdf",sep=""))
|
329 |
|
|
file.remove(list.files(td,full=T))
|
330 |
|
|
if(verbose) print("Finished!!")
|
331 |
|
|
}
|
332 |
|
|
|
333 |
|
|
shrinkpdf("output/lst_lulc.pdf")
|