1
|
###################################################################################
|
2
|
### R code to aquire and process MOD06_L2 cloud data from the MODIS platform
|
3
|
|
4
|
|
5
|
## connect to server of choice
|
6
|
#system("ssh litoria")
|
7
|
#R
|
8
|
|
9
|
library(sp)
|
10
|
library(spgrass6)
|
11
|
library(rgdal)
|
12
|
library(reshape)
|
13
|
library(ncdf4)
|
14
|
library(geosphere)
|
15
|
library(rgeos)
|
16
|
library(multicore)
|
17
|
library(raster)
|
18
|
library(lattice)
|
19
|
library(rgl)
|
20
|
library(hdf5)
|
21
|
library(rasterVis)
|
22
|
library(heR.Misc)
|
23
|
|
24
|
X11.options(type="Xlib")
|
25
|
ncores=20 #number of threads to use
|
26
|
|
27
|
#setwd("/home/adamw/personal/projects/interp")
|
28
|
setwd("/home/adamw/acrobates/projects/interp")
|
29
|
|
30
|
|
31
|
### Produce shapefile of all stations globally
|
32
|
ghcn=read.fwf("data/ghcn/ghcnd-stations.txt",header=F,widths=c(11,-1,8,-1,9,-1,6,-1,2,-1,29,-1,3,-1,3,-1,5),comment.char="")
|
33
|
colnames(ghcn)=c("id","lat","lon","dem","state","name","gsnflag","hcnflag","wmoid")
|
34
|
ghcn=ghcn[,c("id","dem","lat","lon")]
|
35
|
ghcn$id=as.character(ghcn$id)
|
36
|
ghcn$type="GHCN"
|
37
|
|
38
|
### FAO station normals
|
39
|
fao=read.csv("data/stationdata/FAOCLIM2/world.csv")
|
40
|
fao=fao[,c("ID","Elevation","Lat","Lon")]
|
41
|
colnames(fao)=c("id","dem","lat","lon")
|
42
|
fao$type="FAO"
|
43
|
|
44
|
st=rbind.data.frame(ghcn,fao)
|
45
|
coordinates(st)=c("lon","lat")
|
46
|
projection(st)=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs ")
|
47
|
writeOGR(st,dsn="data/",layer="stationlocations",driver="ESRI Shapefile")
|
48
|
|
49
|
plot(fao,pch=16,cex=.2,col="blue")
|
50
|
plot(st,pch=16,cex=.2,col="red",add=T)
|
51
|
|
52
|
|
53
|
###########################
|
54
|
|
55
|
roi=readOGR("data/regions/Test_sites/Oregon.shp","Oregon")
|
56
|
roi_geo=as(roi,"SpatialLines")
|
57
|
roi=spTransform(roi,CRS(" +proj=sinu +lon_0=0 +x_0=0 +y_0=0"))
|
58
|
roil=as(roi,"SpatialLines")
|
59
|
|
60
|
summarydatadir="data/modis/MOD06_climatologies"
|
61
|
|
62
|
|
63
|
|
64
|
##########################
|
65
|
#### explore the data
|
66
|
|
67
|
months=seq(as.Date("2000-01-15"),as.Date("2000-12-15"),by="month")
|
68
|
|
69
|
|
70
|
## load data
|
71
|
cerfiles=list.files(summarydatadir,pattern="CER_mean_.*tif$",full=T); cerfiles
|
72
|
cer=brick(stack(cerfiles))
|
73
|
setZ(cer,months,name="time")
|
74
|
cer@z=list(months)
|
75
|
cer@zname="time"
|
76
|
layerNames(cer) <- as.character(format(months,"%b"))
|
77
|
#cer=projectRaster(from=cer,crs="+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0",method="ngb")
|
78
|
### TODO: change to bilinear!
|
79
|
|
80
|
cotfiles=list.files(summarydatadir,pattern="COT_mean_.*tif$",full=T); cotfiles
|
81
|
cot=brick(stack(cotfiles))
|
82
|
setZ(cot,months,name="time")
|
83
|
cot@z=list(months)
|
84
|
cot@zname="time"
|
85
|
layerNames(cot) <- as.character(format(months,"%b"))
|
86
|
#cot=projectRaster(from=cot,crs="+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0",method="ngb")
|
87
|
cotm=mean(cot,na.rm=T)
|
88
|
### TODO: change to bilinear!
|
89
|
|
90
|
cldfiles=list.files(summarydatadir,pattern="CLD_mean_.*tif$",full=T); cldfiles
|
91
|
cld=brick(stack(cldfiles))
|
92
|
cld[cld==0]=NA
|
93
|
setZ(cld,months,name="time")
|
94
|
cld@z=list(months)
|
95
|
cld@zname="time"
|
96
|
layerNames(cld) <- as.character(format(months,"%b"))
|
97
|
#cot=projectRaster(from=cot,crs="+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0",method="ngb")
|
98
|
cldm=mean(cld,na.rm=T)
|
99
|
### TODO: change to bilinear if reprojecting!
|
100
|
|
101
|
cer20files=list.files(summarydatadir,pattern="CER_P20um_.*tif$",full=T); cer20files
|
102
|
cer20=brick(stack(cer20files))
|
103
|
setZ(cer20,months,name="time")
|
104
|
cer20@z=list(months)
|
105
|
cer20@zname="time"
|
106
|
layerNames(cer20) <- as.character(format(months,"%b"))
|
107
|
#cot=projectRaster(from=cot,crs="+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0",method="ngb")
|
108
|
cotm=mean(cot,na.rm=T)
|
109
|
### TODO: change to bilinear!
|
110
|
|
111
|
|
112
|
### load PRISM data for comparison
|
113
|
prism=brick("data/prism/prism_climate.nc",varname="ppt")
|
114
|
## project to sinusoidal
|
115
|
projection(prism)=CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")
|
116
|
prism=projectRaster(prism,cer)
|
117
|
prism[prism<0]=NA #for some reason NAvalue() wasn't working
|
118
|
setZ(prism,months,name="time")
|
119
|
prism@z=list(months)
|
120
|
prism@zname="time"
|
121
|
layerNames(prism) <- as.character(format(months,"%b"))
|
122
|
|
123
|
#### build a pixel by variable matrix
|
124
|
vars=c("cer","cer20","cld","cot","prism")
|
125
|
bd=melt(as.matrix(vars[1]))
|
126
|
colnames(bd)=c("cell","month",vars[1])
|
127
|
for(v in vars[-1]) {print(v); bd[,v]=melt(as.matrix(get(v)))$value}
|
128
|
bd=bd[!is.na(bd$cer)|is.na(bd$prism),]
|
129
|
|
130
|
## Summarize annual metrics for full rasters
|
131
|
|
132
|
### get all variables from all months
|
133
|
c01=brick(mclapply(vars,function(v) projectRaster(get(v),crs="+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")))
|
134
|
layerNames(m01)=paste(vars,months,sep="_")
|
135
|
|
136
|
m01=brick(mclapply(vars,function(v) mean(get(v))))#mean(cer),mean(cld),mean(cot),mean(prism))
|
137
|
layerNames(m01)=vars
|
138
|
m01=projectRaster(from=m01,crs="+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")
|
139
|
m01=crop(m01,extent(-125,-115,41,47))
|
140
|
### get station data, subset to stations in region, and transform to sinusoidal
|
141
|
load("data/ghcn/roi_ghcn.Rdata")
|
142
|
load("data/allstations.Rdata")
|
143
|
|
144
|
st2_sin=spTransform(st2,CRS(projection(cer)))
|
145
|
|
146
|
d2=d[d$variable=="ppt"&d$date>=as.Date("2000-01-01"),]
|
147
|
d2=d2[,-grep("variable",colnames(d2)),]
|
148
|
st2=st[st$id%in%d$id,]
|
149
|
#st2=spTransform(st2,CRS(" +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +towgs84=0,0,0"))
|
150
|
d2[,c("lon","lat")]=coordinates(st2)[match(d2$id,st2$id),]
|
151
|
d2$elev=st2$elev[match(d2$id,st2$id)]
|
152
|
d2$month=format(d2$date,"%m")
|
153
|
d2$value=d2$value/10 #convert to mm
|
154
|
|
155
|
### extract MOD06 data for each station
|
156
|
stcer=extract(cer,st2)#;colnames(stcer)=paste("cer_mean_",1:12,sep="")
|
157
|
stcer20=extract(cer20,st2)#;colnames(stcer)=paste("cer_mean_",1:12,sep="")
|
158
|
stcot=extract(cot,st2)#;colnames(stcot)=paste("cot_mean_",1:12,sep="")
|
159
|
stcld=extract(cld,st2)#;colnames(stcld)=paste("cld_mean_",1:12,sep="")
|
160
|
mod06=cbind.data.frame(id=st2$id,lat=st2$lat,lon=st2$lon,stcer,stcer20,stcot,stcld)
|
161
|
mod06l=melt(mod06,id.vars=c("id","lon","lat"))
|
162
|
mod06l[,c("variable","moment","month")]=do.call(rbind,strsplit(as.character(mod06l$variable),"_"))
|
163
|
mod06l=as.data.frame(cast(mod06l,id+lon+lat+month~variable+moment,value="value"))
|
164
|
|
165
|
### Identify stations that have < 10 years of data
|
166
|
cnts=cast(d2,id~.,fun=function(x) length(x[!is.na(x)]),value="count");colnames(cnts)[colnames(cnts)=="(all)"]="count"
|
167
|
summary(cnts)
|
168
|
## drop them
|
169
|
d2=d2[d2$id%in%cnts$id[cnts$count>=365*10],]
|
170
|
|
171
|
|
172
|
### generate monthly means of station data
|
173
|
dc=cast(d2,id+lon+lat+elev~month,value="value",fun=function(x) mean(x,na.rm=T)*30)
|
174
|
dcl=melt(dc,id.vars=c("id","lon","lat","elev"),value="ppt")
|
175
|
colnames(dcl)[colnames(dcl)=="value"]="ppt"
|
176
|
|
177
|
|
178
|
|
179
|
## merge station data with mod06
|
180
|
mod06s=merge(dcl,mod06l)
|
181
|
|
182
|
|
183
|
### draw some plots
|
184
|
gq=function(x,n=10,cut=F) {
|
185
|
if(!cut) return(unique(quantile(x,seq(0,1,len=n+1),na.rm=T)))
|
186
|
if(cut) return(cut(x,unique(quantile(x,seq(0,1,len=n+1),na.rm=T))))
|
187
|
}
|
188
|
|
189
|
### add some additional variables
|
190
|
mod06s$month=factor(mod06s$month,labels=format(as.Date(paste("2000",1:12,"15",sep="-")),"%b"))
|
191
|
mod06s$lppt=log(mod06s$ppt)
|
192
|
mod06s$glon=cut(mod06s$lon,gq(mod06s$lon,n=5),include.lowest=T,ordered=T)#gq(mod06s$lon,n=3))
|
193
|
mod06s$glon2=cut(mod06s$lon,breaks=c(-125,-122,-115),labels=c("Coastal","Inland"),include.lowest=T,ordered=T)#gq(mod06s$lon,n=3))
|
194
|
mod06s$gelev=cut(mod06s$elev,breaks=gq(mod06s$elev,n=3),labels=c("Low","Mid","High"),include.lowest=T,ordered=T)
|
195
|
mod06s$gbin=factor(paste(mod06s$gelev,mod06s$glon2,sep="_"),levels=c("Low_Coastal","Mid_Coastal","High_Coastal","Low_Inland","Mid_Inland","High_Inland"),ordered=T)
|
196
|
mod06s$LWP_mean=(2/3)*mod06s$CER_mean*mod06s$COT_mean
|
197
|
|
198
|
## melt it
|
199
|
mod06sl=melt(mod06s[,!grepl("lppt",colnames(mod06s))],id.vars=c("id","lon","lat","elev","month","ppt","glon","glon2","gelev","gbin"))
|
200
|
levels(mod06sl$variable)=c("Effective Radius (um)","Very Cloudy Days (%)","Cloudy Days (%)","Optical Thickness (%)","Liquid Water Path")
|
201
|
|
202
|
###################################################################
|
203
|
###################################################################
|
204
|
|
205
|
bgyr=colorRampPalette(c("blue","green","yellow","red"))
|
206
|
|
207
|
X11.options(type="cairo")
|
208
|
pdf("output/MOD06_summary.pdf",width=11,height=8.5)
|
209
|
|
210
|
# % cloudy maps
|
211
|
title="Cloudiness (% cloudy days) "
|
212
|
at=unique(quantile(as.matrix(cld),seq(0,1,len=100),na.rm=T))
|
213
|
p=levelplot(cld,xlab.top=title,at=at,col.regions=bgyr(length(at)))+layer(sp.lines(roil, lwd=1.2, col='black'))
|
214
|
print(p)
|
215
|
#bwplot(cer,main=title,ylab="Cloud Effective Radius (microns)")
|
216
|
|
217
|
# CER maps
|
218
|
title="Cloud Effective Radius (microns)"
|
219
|
at=quantile(as.matrix(cer),seq(0,1,len=100),na.rm=T)
|
220
|
p=levelplot(cer,xlab.top=title,at=at,col.regions=bgyr(length(at)))+layer(sp.lines(roil, lwd=1.2, col='black'))
|
221
|
print(p)
|
222
|
#bwplot(cer,main=title,ylab="Cloud Effective Radius (microns)")
|
223
|
|
224
|
# COT maps
|
225
|
title="Cloud Optical Thickness (%)"
|
226
|
at=quantile(as.matrix(cot),seq(0,1,len=100),na.rm=T)
|
227
|
p=levelplot(cot,xlab.top=title,at=at,col.regions=bgyr(length(at)))+layer(sp.lines(roil, lwd=0.8, col='black'))
|
228
|
print(p)
|
229
|
#bwplot(cot,xlab.top=title,ylab="Cloud Optical Thickness (%)")
|
230
|
|
231
|
###########################################
|
232
|
#### compare with PRISM data
|
233
|
|
234
|
at=quantile(as.matrix(subset(m01,subset=1)),seq(0,1,len=100),na.rm=T)
|
235
|
p1=levelplot(subset(m01,subset=1),xlab.top="Effective Radius (um)",at=at,col.regions=bgyr(length(at)),margin=F,
|
236
|
)+layer(sp.lines(roi_geo, lwd=1.2, col='black'))
|
237
|
at=quantile(as.matrix(subset(m01,subset=2)),seq(0,1,len=100),na.rm=T)
|
238
|
p2=levelplot(subset(m01,subset=2),xlab.top="Cloudy days (%)",at=at,col.regions=bgyr(length(at)),margin=F,
|
239
|
)+layer(sp.lines(roi_geo, lwd=1.2, col='black'))
|
240
|
at=quantile(as.matrix(subset(m01,subset=3)),seq(0,1,len=100),na.rm=T)
|
241
|
p3=levelplot(subset(m01,subset=3),xlab.top="Optical Thickness (%)",at=at,col.regions=bgyr(length(at)),margin=F,
|
242
|
)+layer(sp.lines(roi_geo, lwd=1.2, col='black'))
|
243
|
at=quantile(as.matrix(subset(m01,subset=4)),seq(0,1,len=100),na.rm=T)
|
244
|
p4=levelplot(subset(m01,subset=4),xlab.top="PRISM MAP",at=at,col.regions=bgyr(length(at)),margin=F,
|
245
|
)+layer(sp.lines(roi_geo, lwd=1.2, col='black'))
|
246
|
|
247
|
print(p1,split=c(1,1,2,2))
|
248
|
print(p2,split=c(1,2,2,2),new=F)
|
249
|
print(p3,split=c(2,1,2,2),new=F)
|
250
|
print(p4,split=c(2,2,2,2),new=F)
|
251
|
|
252
|
### compare COT and PRISM
|
253
|
print(p3,split=c(1,1,2,1))
|
254
|
print(p4,split=c(2,1,2,1),new=F)
|
255
|
|
256
|
|
257
|
### sample to speed up processing
|
258
|
s=sample(1:nrow(bd),10000)
|
259
|
|
260
|
## melt it to ease comparisons
|
261
|
bdl=melt(bd[s,],measure.vars=c("cer","cld","cot"))
|
262
|
|
263
|
combineLimits(useOuterStrips(xyplot(prism~value|variable+month,data=bdl,pch=16,cex=.2,scales=list(y=list(log=T),x=list(relation="free")),
|
264
|
ylab="PRISM Monthly mean precipitation (mm)",xlab="MOD06 metric",main="PRISM vs. MOD06 (mean monthly ppt)")))+
|
265
|
layer(panel.abline(lm(y~x),col="red"))+
|
266
|
layer(panel.text(0,2.5,paste("R2=",round(summary(lm(y~x))$r.squared,2))))
|
267
|
|
268
|
|
269
|
### Comparison at station values
|
270
|
at=quantile(as.matrix(cotm),seq(0,1,len=100),na.rm=T)
|
271
|
p=levelplot(cotm, layers=1,at=at,col.regions=bgyr(length(at)),main="Mean Annual Cloud Optical Thickness",FUN.margin=function(x) 0)+
|
272
|
layer(sp.lines(roil, lwd=1.2, col='black'))+layer(sp.points(st2_sin, pch=16, col='black'))
|
273
|
print(p)
|
274
|
|
275
|
### monthly comparisons of variables
|
276
|
#mod06sl=melt(mod06s,measure.vars=c("ppt","COT_mean","CER_mean","CER_P20um"))
|
277
|
#bwplot(value~month|variable,data=mod06sl,cex=.5,pch=16,col="black",scales=list(y=list(relation="free")),layout=c(1,3))
|
278
|
#splom(mod06s[grep("CER|COT|CLD|ppt",colnames(mod06s))],cex=.2,pch=16,main="Scatterplot matrix of MOD06 products")
|
279
|
|
280
|
### run some regressions
|
281
|
#plot(log(ppt)~COT_mean,data=mod06s)
|
282
|
#summary(lm(log(ppt)~COT_mean*month,data=mod06s))
|
283
|
|
284
|
## ppt~metric with longitude bins
|
285
|
xyplot(ppt~value|variable,groups=glon,data=mod06sl,
|
286
|
scales=list(y=list(log=T),x=list(relation="free",log=F)),
|
287
|
par.settings = list(superpose.symbol = list(col=bgyr(5),pch=16,cex=.5)),auto.key=list(space="top",title="Station Longitude"),
|
288
|
main="Comparison of MOD06_L2 and Precipitation Monthly Climatologies",ylab="Mean Monthly Station Precipitation (mm)",xlab="MOD06_L2 Product",layout=c(5,1))+
|
289
|
layer(panel.text(9,2.5,label="Coastal stations",srt=30,cex=1.3,col="blue"),columns=1)+
|
290
|
layer(panel.text(13,.9,label="Inland stations",srt=10,cex=1.3,col="red"),columns=1)+
|
291
|
layer(panel.abline(lm(y~x),col="red"))+
|
292
|
layer(panel.text(0,0,paste("R2=",round(summary(lm(y~x))$r.squared,2)),pos=4,cex=.5,col="grey"))
|
293
|
|
294
|
|
295
|
## ppt~metric with longitude bins
|
296
|
#CairoPNG("output/COT.png",width=10,height=5,units="in",dpi=300,pointsize=20)
|
297
|
#png("output/COT.png",width=10,height=5,units="in",res=150)
|
298
|
#trellis.par.set("fontsize",12)
|
299
|
at=quantile(as.matrix(subset(m01,subset=3)),seq(0,1,len=100),na.rm=T)
|
300
|
p1=levelplot(subset(m01,subset=3),xlab.top="Optical Thickness (%)",at=at,col.regions=bgyr(length(at)),margin=F,
|
301
|
)+layer(sp.lines(roi_geo, lwd=1.2, col='black'))+layer(sp.points(st2, cex=.5,col='black'))
|
302
|
at=quantile(as.matrix(subset(m01,subset=3)),seq(0,1,len=100),na.rm=T)
|
303
|
|
304
|
at=quantile(as.matrix(subset(m01,subset=4)),seq(0,1,len=100),na.rm=T)
|
305
|
p2=levelplot(subset(m01,subset=4),xlab.top="PRISM MAP",at=at,col.regions=bgyr(length(at)),margin=F,
|
306
|
)+layer(sp.lines(roi_geo, lwd=1.2, col='black'))+layer(sp.points(st2, cex=.5, col='black'))
|
307
|
|
308
|
p3=xyplot(ppt~value,groups=glon,data=mod06sl[mod06sl$variable=="Optical Thickness (%)",],
|
309
|
scales=list(y=list(log=T),x=list(relation="free",log=F)),
|
310
|
par.settings = list(superpose.symbol = list(col=bgyr(5),pch=16,cex=.3)),auto.key=list(space="right",title="Station \n Longitude"),
|
311
|
ylab="Mean Monthly Station Precipitation (mm)",xlab="Cloud Optical Thickness from MOD06_L2 (%)",layout=c(1,1))+
|
312
|
layer(panel.text(9,2.6,label="Coastal stations",srt=10,cex=1.3,col="blue"),columns=1)+
|
313
|
layer(panel.text(13,.95,label="Inland stations",srt=10,cex=1.3,col="red"),columns=1)
|
314
|
|
315
|
p4=xyplot(ppt~value,groups=glon,data=mod06sl[mod06sl$variable=="Very Cloudy Days (%)",],
|
316
|
scales=list(y=list(log=T),x=list(relation="free",log=F)),
|
317
|
par.settings = list(superpose.symbol = list(col=bgyr(5),pch=16,cex=.3)),auto.key=list(space="right",title="Station \n Longitude"),
|
318
|
ylab="Mean Monthly Station Precipitation (mm)",xlab="Proportion days with Cloud Effective Radius >20um from MOD06_L2 (%)",layout=c(1,1))+
|
319
|
layer(panel.text(9,2.6,label="Coastal stations",srt=10,cex=1.3,col="blue"),columns=1)+
|
320
|
layer(panel.text(13,.95,label="Inland stations",srt=10,cex=1.3,col="red"),columns=1)
|
321
|
|
322
|
#save(p1,p2,p3,file="plotdata.Rdata")
|
323
|
#load("plotdata.Rdata")
|
324
|
|
325
|
#CairoPDF("output/MOD06_Summaryfig.pdf",width=11,height=8.5)
|
326
|
print(p3,position=c(0,0,1,.5),save.object=F)
|
327
|
#print(p4,position=c(0,0,1,.5),save.object=F)
|
328
|
print(p1,split=c(1,1,2,2),new=F)
|
329
|
print(p2,split=c(2,1,2,2),new=F)
|
330
|
#dev.off()
|
331
|
#system("convert output/MOD06_Summaryfig.pdf output/MOD06_Summaryfig.png")
|
332
|
#dev.off()
|
333
|
|
334
|
## with elevation
|
335
|
# xyplot(ppt~value|variable,groups=gbin,data=mod06sl,
|
336
|
# scales=list(y=list(log=T),x=list(relation="free")),
|
337
|
# par.settings = list(superpose.symbol = list(col=c(rep("blue",3),rep("red",3)),pch=rep(c(3,4,8),2),cex=.5)),auto.key=list(space="right",title="Station Longitude"),
|
338
|
# main="Comparison of MOD06_L2 and Precipitation Monthly Climatologies",ylab="Precipitation",xlab="MOD06_L2 Product",layout=c(3,1))+
|
339
|
# layer(panel.text(9,2.5,label="Coastal stations",srt=30,cex=1.3,col="blue"),columns=1)+
|
340
|
# layer(panel.text(13,.9,label="Inland stations",srt=10,cex=1.3,col="red"),columns=1)
|
341
|
|
342
|
## with elevation and longitude bins
|
343
|
combineLimits(useOuterStrips(xyplot(ppt~value|variable+gbin,data=mod06sl,
|
344
|
scales=list(y=list(log=T),x=list(relation="free")),col="black",pch=16,cex=.5,type=c("p","r"),
|
345
|
main="Comparison of MOD06_L2 and Precipitation Monthly Climatologies",ylab="Precipitation",xlab="MOD06_L2 Product")))+
|
346
|
layer(panel.xyplot(x,y,type="r",col="red"))
|
347
|
|
348
|
## *** MOD06 vars vs precipitation by month, colored by longitude
|
349
|
combineLimits(useOuterStrips(xyplot(ppt~value|month+variable,groups=glon,data=mod06sl,cex=.5,pch=16,
|
350
|
scales=list(y=list(log=T),x=list(relation="free")),
|
351
|
par.settings = list(superpose.symbol = list(pch =16, col=bgyr(5),cex=1)),auto.key=list(space="top",title="Station Longitude"),
|
352
|
main="Comparison of MOD06_L2 and Precipitation Monthly Climatologies",ylab="Precipitation",xlab="MOD06_L2 Product")),
|
353
|
margin.x=1,adjust.labels=F)+
|
354
|
layer(panel.abline(lm(y~x),col="red"))+
|
355
|
layer(panel.text(0,0,paste("R2=",round(summary(lm(y~x))$r.squared,2)),pos=4,cex=.5,col="grey"))
|
356
|
|
357
|
|
358
|
xyplot(ppt~CLD_mean|id,data=mod06s,panel=function(x,y,group){
|
359
|
panel.xyplot(x,y,type=c("r"),cex=.5,pch=16,col="red")
|
360
|
panel.xyplot(x,y,type=c("p"),cex=.5,pch=16,col="black")
|
361
|
} ,scales=list(y=list(log=T)),strip=F,main="Monthly Mean Precipitation and % Cloudy by station",
|
362
|
sub="Each panel is a station, each point is a monthly mean",
|
363
|
ylab="Precipitation (mm, log axis)",xlab="% of Cloudy Days")+
|
364
|
layer(panel.text(.5,.5,round(summary(lm(y~x))$r.squared,2),pos=4,cex=.75,col="grey"))
|
365
|
|
366
|
xyplot(ppt~CER_mean|id,data=mod06s,panel=function(x,y,group){
|
367
|
panel.xyplot(x,y,type=c("r"),cex=.5,pch=16,col="red")
|
368
|
panel.xyplot(x,y,type=c("p"),cex=.5,pch=16,col="black")
|
369
|
} ,scales=list(y=list(log=T)),strip=F,main="Monthly Mean Precipitation and Cloud Effective Radius by station",sub="Each panel is a station, each point is a monthly mean",ylab="Precipitation (mm, log axis)",xlab="Mean Monthly Cloud Effective Radius (mm)")+
|
370
|
layer(panel.text(10,.5,round(summary(lm(y~x))$r.squared,2),pos=4,cex=.75,col="grey"))
|
371
|
|
372
|
xyplot(ppt~COT_mean|id,data=mod06s,panel=function(x,y,group){
|
373
|
panel.xyplot(x,y,type=c("r"),cex=.5,pch=16,col="red")
|
374
|
panel.xyplot(x,y,type=c("p"),cex=.5,pch=16,col="black")
|
375
|
} ,scales=list(y=list(log=T)),strip=F,main="Monthly Mean Precipitation and Cloud Optical Thickness by station",
|
376
|
sub="Each panel is a station, each point is a monthly mean \n Number in lower right of each panel is R^2",
|
377
|
ylab="Precipitation (mm, log axis)",xlab="Mean Monthly Cloud Optical Thickness (%)")+
|
378
|
layer(panel.text(10,.5,round(summary(lm(y~x))$r.squared,2),pos=4,cex=.75,col="grey"))
|
379
|
|
380
|
xyplot(ppt~LWP_mean|id,data=mod06s,panel=function(x,y,group){
|
381
|
panel.xyplot(x,y,type=c("r"),cex=.5,pch=16,col="red")
|
382
|
panel.xyplot(x,y,type=c("p"),cex=.5,pch=16,col="black")
|
383
|
} ,scales=list(y=list(log=T)),strip=F,main="Monthly Mean Precipitation and Liquid Water Path by station",
|
384
|
sub="Each panel is a station, each point is a monthly mean \n Number in lower right of each panel is R^2",
|
385
|
ylab="Precipitation (mm, log axis)",xlab="Mean Monthly Liquid Water Path")+
|
386
|
layer(panel.text(10,.5,round(summary(lm(y~x))$r.squared,2),pos=4,cex=.75,col="grey"))
|
387
|
|
388
|
### Calculate the slope of each line
|
389
|
mod06s.sl=dapply(mod06s,list(id=mod06s$id),function(x){
|
390
|
lm1=lm(log(x$ppt)~x$CER_mean,)
|
391
|
data.frame(lat=x$lat[1],lon=x$lon[1],elev=x$elev[1],intcpt=coefficients(lm1)[1],cer=coefficients(lm1)[2],r2=summary(lm1)$r.squared)
|
392
|
})
|
393
|
mod06s.sl$cex=gq(mod06s.sl$r2,n=5,cut=T)
|
394
|
mod06s.sl$cer.s=gq(mod06s.sl$cer,n=5,cut=T)
|
395
|
|
396
|
### and plot it on a map
|
397
|
xyplot(lat~lon,group=cer.s,data=mod06s.sl,par.settings = list(superpose.symbol = list(pch =16, col=bgyr(5),cex=1)),auto.key=list(space="right",title="Slope Coefficient"),asp=1,
|
398
|
main="Slopes of linear regressions {log(ppt)~CloudEffectiveRadius}")+
|
399
|
layer(sp.lines(roi_geo, lwd=1.2, col='black'))
|
400
|
|
401
|
### look for relationships with longitude
|
402
|
xyplot(cer~lon,group=cut(mod06s.sl$elev,gq(mod06s.sl$elev,n=5)),data=mod06s.sl,
|
403
|
par.settings = list(superpose.symbol = list(col=bgyr(5),pch=16,cex=1)),auto.key=list(space="right",title="Station Elevation"),
|
404
|
ylab="Slope of lm(ppt~EffectiveRadius)",xlab="Longitude",main="Precipitation~Effective Radius relationship by latitude")
|
405
|
|
406
|
|
407
|
############################################################
|
408
|
### simple regression to get spatial residuals
|
409
|
m="01"
|
410
|
mod06s2=mod06s#[mod06s$month==m,]
|
411
|
|
412
|
lm1=lm(log(ppt)~CER_mean*month*lon,data=mod06s2); summary(lm1)
|
413
|
mod06s2$pred=exp(predict(lm1,mod06s2))
|
414
|
mod06s2$resid=mod06s2$pred-mod06s2$ppt
|
415
|
mod06s2$residg=gq(mod06s2$resid,n=5,cut=T)
|
416
|
mod06s2$presid=mod06s2$resid/mod06s2$ppt
|
417
|
|
418
|
for(l in c(F,T)){
|
419
|
## all months
|
420
|
xyplot(pred~ppt,groups=gelev,data=mod06s2,
|
421
|
par.settings = list(superpose.symbol = list(col=bgyr(3),pch=16,cex=.75)),auto.key=list(space="right",title="Station Elevation"),
|
422
|
scales=list(log=l),
|
423
|
ylab="Predicted Mean Monthly Precipitation (mm)",xlab="Observed Mean Monthly Precipitation (mm)",main="Predicted vs. Observed for Simple Model",
|
424
|
sub="Red line is y=x")+
|
425
|
layer(panel.abline(0,1,col="red"))
|
426
|
|
427
|
## month by month
|
428
|
print(xyplot(pred~ppt|month,groups=gelev,data=mod06s2,
|
429
|
par.settings = list(superpose.symbol = list(col=bgyr(3),pch=16,cex=.75)),auto.key=list(space="right",title="Station Elevation"),
|
430
|
scales=list(log=l),
|
431
|
ylab="Predicted Mean Monthly Precipitation (mm)",xlab="Observed Mean Monthly Precipitation (mm)",main="Predicted vs. Observed for Simple Model",
|
432
|
sub="Red line is y=x")+
|
433
|
layer(panel.abline(0,1,col="red"))
|
434
|
)}
|
435
|
|
436
|
## residuals by month
|
437
|
xyplot(lat~lon|month,group=residg,data=mod06s2,
|
438
|
par.settings = list(superpose.symbol = list(pch =16, col=bgyr(5),cex=.5)),
|
439
|
auto.key=list(space="right",title="Residuals"),
|
440
|
main="Spatial plot of monthly residuals")+
|
441
|
layer(sp.lines(roi_geo, lwd=1.2, col='black'))
|
442
|
|
443
|
|
444
|
dev.off()
|
445
|
|
446
|
|
447
|
####################################
|
448
|
#### build table comparing various metrics
|
449
|
mods=data.frame(
|
450
|
models=c(
|
451
|
"log(ppt)~CER_mean",
|
452
|
"log(ppt)~CLD_mean",
|
453
|
"log(ppt)~COT_mean",
|
454
|
"log(ppt)~CER_mean*month",
|
455
|
"log(ppt)~CLD_mean*month",
|
456
|
"log(ppt)~COT_mean*month",
|
457
|
"log(ppt)~CER_mean*month*lon",
|
458
|
"log(ppt)~CLD_mean*month*lon",
|
459
|
"log(ppt)~COT_mean*month*lon",
|
460
|
"ppt~CER_mean*month*lon",
|
461
|
"ppt~CLD_mean*month*lon",
|
462
|
"ppt~COT_mean*month*lon"),stringsAsFactors=F)
|
463
|
|
464
|
mods$r2=
|
465
|
do.call(rbind,lapply(1:nrow(mods),function(i){
|
466
|
lm1=lm(as.formula(mods$models[i]),data=mod06s2)
|
467
|
summary(lm1)$r.squared}))
|
468
|
|
469
|
mods
|
470
|
|
471
|
|
472
|
|
473
|
|
474
|
|
475
|
|
476
|
|
477
|
|
478
|
load("data/modis/pointsummary.Rdata")
|
479
|
|
480
|
|
481
|
dsl=melt(ds,id.vars=c("id","date","ppt","lon","lat"),measure.vars= c("Cloud_Water_Path","Cloud_Effective_Radius","Cloud_Optical_Thickness"))
|
482
|
|
483
|
dsl=dsl[!is.nan(dsl$value),]
|
484
|
|
485
|
|
486
|
|
487
|
|
488
|
####
|
489
|
## mean annual precip
|
490
|
dp=d[d$variable=="ppt",]
|
491
|
dp$year=format(dp$date,"%Y")
|
492
|
dm=tapply(dp$value,list(id=dp$id,year=dp$year),sum,na.rm=T)
|
493
|
dms=apply(dm,1,mean,na.rm=T)
|
494
|
dms=data.frame(id=names(dms),ppt=dms/10)
|
495
|
|
496
|
dslm=tapply(dsl$value,list(id=dsl$id,variable=dsl$variable),mean,na.rm=T)
|
497
|
dslm=data.frame(id=rownames(dslm),dslm)
|
498
|
|
499
|
dms=merge(dms,dslm)
|
500
|
dmsl=melt(dms,id.vars=c("id","ppt"))
|
501
|
|
502
|
summary(lm(ppt~Cloud_Effective_Radius,data=dms))
|
503
|
summary(lm(ppt~Cloud_Water_Path,data=dms))
|
504
|
summary(lm(ppt~Cloud_Optical_Thickness,data=dms))
|
505
|
summary(lm(ppt~Cloud_Effective_Radius+Cloud_Water_Path+Cloud_Optical_Thickness,data=dms))
|
506
|
|
507
|
|
508
|
#### draw some plots
|
509
|
#pdf("output/MOD06_summary.pdf",width=11,height=8.5)
|
510
|
png("output/MOD06_summary_%d.png",width=1024,height=780)
|
511
|
|
512
|
## daily data
|
513
|
xyplot(value~ppt/10|variable,data=dsl,
|
514
|
scales=list(relation="free"),type=c("p","r"),
|
515
|
pch=16,cex=.5,layout=c(3,1))
|
516
|
|
517
|
|
518
|
densityplot(~value|variable,groups=cut(dsl$ppt,c(0,50,100,500)),data=dsl,auto.key=T,
|
519
|
scales=list(relation="free"),plot.points=F)
|
520
|
|
521
|
## annual means
|
522
|
|
523
|
xyplot(value~ppt|variable,data=dmsl,
|
524
|
scales=list(relation="free"),type=c("p","r"),pch=16,cex=0.5,layout=c(3,1),
|
525
|
xlab="Mean Annual Precipitation (mm)",ylab="Mean value")
|
526
|
|
527
|
densityplot(~value|variable,groups=cut(dsl$ppt,c(0,50,100,500)),data=dmsl,auto.key=T,
|
528
|
scales=list(relation="free"),plot.points=F)
|
529
|
|
530
|
|
531
|
## plot some swaths
|
532
|
|
533
|
nc1=raster(fs$path[3],varname="Cloud_Effective_Radius")
|
534
|
nc2=raster(fs$path[4],varname="Cloud_Effective_Radius")
|
535
|
nc3=raster(fs$path[5],varname="Cloud_Effective_Radius")
|
536
|
|
537
|
nc1[nc1<=0]=NA
|
538
|
nc2[nc2<=0]=NA
|
539
|
nc3[nc3<=0]=NA
|
540
|
|
541
|
plot(roi)
|
542
|
plot(nc3)
|
543
|
|
544
|
plot(nc1,add=T)
|
545
|
plot(nc2,add=T)
|
546
|
|
547
|
|
548
|
dev.off()
|
549
|
|
550
|
|