|
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
|
library(xtable)
|
|
17
|
library(ellipse) # for correlation matrix
|
|
18
|
library(maptools) # for rgshhs
|
|
19
|
|
|
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
|
|
|
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
|
|
|
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
|
##topo=brick(list(topo,demb))
|
|
81
|
|
|
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
|
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
|
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
|
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
|
### 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
|
save(s,file=paste("output/mod_sample.Rdata",sep=""))
|
|
136
|
|
|
137
|
|
|
138
|
### spatial regression to fit lulc coefficients
|
|
139
|
Sys.setenv(MKL_NUM_THREADS=24)
|
|
140
|
system("export MKL_NUM_THREADS=24")
|
|
141
|
|
|
142
|
niter=12500
|
|
143
|
lapply(months,function(m){
|
|
144
|
print(paste("################################### Starting month ",m))
|
|
145
|
f1=formula(paste(m,"~logDTCkm+y+Shrub+Grass+Crop+Mosaic+Urban+Barren+Snow+Wetland+dem+eastwest+northsouth",sep=""))
|
|
146
|
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
|
n.samples=niter,sub.samples=c(2500,niter,10),verbose=TRUE, n.report=100)
|
|
152
|
assign(paste("mod_",m,sep=""),m1)
|
|
153
|
save(list=paste("mod_",m,sep=""),file=paste("output/mod_",m,".Rdata",sep=""))
|
|
154
|
})
|
|
155
|
|
|
156
|
### 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
|
|
|
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
|
|
|
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
|
#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
|
sub="One point per pixel, per month",ylab="LST",
|
|
220
|
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
|
## 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
|
|
|
232
|
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
|
|
|
239
|
|
|
240
|
|
|
241
|
##########################################
|
|
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")
|
