Project

General

Profile

Download (21.8 KB) Statistics
| Branch: | Revision:
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
library(car)
24

    
25
X11.options(type="Xlib")
26
ncores=20  #number of threads to use
27

    
28
setwd("/home/adamw/acrobates/projects/interp")
29

    
30
psin=CRS("+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs")
31

    
32
## get MODLAND tile information
33
tb=read.table("http://landweb.nascom.nasa.gov/developers/sn_tiles/sn_bound_10deg.txt",skip=6,nrows=648,header=T)
34
tb$tile=paste("h",sprintf("%02d",tb$ih),"v",sprintf("%02d",tb$iv),sep="")
35
save(tb,file="modlandTiles.Rdata")
36

    
37
tile="h11v08"   #can move this to submit script if needed
38
#tile="h09v04"   #oregon
39

    
40
tile_bb=tb[tb$tile==tile,] ## identify tile of interest
41
roi_ll=extent(tile_bb$lon_min,tile_bb$lon_max,tile_bb$lat_min,tile_bb$lat_max) 
42
#roi=spTransform(roi,psin)
43
#roil=as(roi,"SpatialLines")
44

    
45
dmod06="data/modis/mod06/summary"
46

    
47

    
48
##########################
49
#### explore the data
50

    
51
months=seq(as.Date("2000-01-15"),as.Date("2000-12-15"),by="month")
52

    
53
getmod06<-function(variable){
54
  d=brick(list.files(dmod06,pattern=paste("MOD06_",tile,".nc",sep=""),full=T),varname=toupper(variable))
55
  projection(d)=psin
56
  setZ(d,format(as.Date(d@z$Date),"%m"),name="time")
57
#  d@z=list(months)
58
  layerNames(d) <- as.character(format(as.Date(d@z$Date),"%b"))
59
  return(d)
60
}
61

    
62
cer=getmod06("cer")
63
cld=getmod06("cld")
64
cot=getmod06("cot")
65

    
66
pcol=colorRampPalette(c("brown","red","yellow","darkgreen"))
67
#levelplot(cer,col.regions=pcol(20))
68

    
69
## load WorldClim data for comparison (download then uncompress)
70
#system("wget -P data/worldclim/ http://biogeo.ucdavis.edu/data/climate/worldclim/1_4/grid/cur/prec_30s_bil.zip",wait=F)
71
#system("wget -P data/worldclim/ http://biogeo.ucdavis.edu/data/climate/worldclim/1_4/grid/cur/alt_30s_bil.zip",wait=F)
72

    
73
### load WORLDCLIM data for comparison
74
wc=stack(list.files("data/worldclim/prec_30s_bil/",pattern="bil$",full=T)[c(4:12,1:3)])
75
projection(wc)=CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")
76
wc=crop(wc,roi_ll)
77
wc[wc==55537]=NA
78
wc=projectRaster(wc,cer)#crs=projection(psin))
79
setZ(wc,months,name="time")
80
wc@z=list(months)
81
layerNames(wc) <- as.character(format(months,"%b"))
82
writeRaster(wc,file=paste("data/tiles/",tile,"/worldclim_",tile,".tif",sep=""),format="GTiff")
83

    
84
### load WORLDCLIM elevation 
85
dem=raster(list.files("data/worldclim/alt_30s_bil/",pattern="bil$",full=T))
86
projection(dem)=CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")
87
dem=crop(dem,roi_ll)
88
dem[dem>60000]=NA
89
dem=projectRaster(dem,cer)
90
writeRaster(dem,file=paste("data/tiles/",tile,"/dem_",tile,".tif",sep=""),format="GTiff")
91

    
92

    
93
### get station data, subset to stations in region, and transform to sinusoidal
94
dm=readOGR(paste("data/tiles/",tile,sep=""),paste("station_monthly_",tile,"_PRCP",sep=""))
95
xyplot(latitude~longitude|month,data=dm@data)
96
dm2=spTransform(dm,CRS(projection(cer)))
97
dm2@data[,c("x","y")]=coordinates(dm2)
98

    
99
### extract MOD06 data for each station
100
stcer=extract(cer,dm2,fun=mean);colnames(stcer)=paste("cer_mean_",as.numeric(format(as.Date(cer@z$Date),"%m")),sep="")
101
#stcer20=extract(cer20,st2)#;colnames(stcer)=paste("cer_mean_",1:12,sep="")
102
stcot=extract(cot,dm2);colnames(stcot)=paste("cot_mean_",as.numeric(format(as.Date(cot@z$Date),"%m")),sep="")
103
stcld=extract(cld,dm2);colnames(stcld)=paste("cld_mean_",as.numeric(format(as.Date(cld@z$Date),"%m")),sep="")
104
stdem=extract(dem,dm2)
105
mod06=cbind.data.frame(station=dm$station,stcer,stcot,stcld)
106
mod06l=melt(mod06,id.vars=c("station"));colnames(mod06l)[grep("value",colnames(mod06l))]="mod06"
107
mod06l[,c("variable","moment","month")]=do.call(rbind,strsplit(as.character(mod06l$variable),"_"))
108
mod06l=unique(mod06l)
109
mod06l=cast(mod06l,station+moment+month~variable,value="mod06")
110
mod06l=merge(dm2@data,mod06l,by=c("station","month"))
111
mod06l=mod06l[!is.na(mod06l$cer),]
112

    
113
                                        #mod06l=melt(mod06,id.vars=c("station","longitude","latitude","elevation","month","count","value"))
114
#mod06l[,c("variable","moment","month2")]=do.call(rbind,strsplit(as.character(mod06l$variable),"_"))
115
#mod06l=as.data.frame(cast(mod06l,station+longitude+latitude+month~variable,value="value.1"))
116
#mod06l=mod06l[mod06l$month==mod06l$month2&!is.na(mod06l$value.1),]
117
mod06l=mod06l[order(mod06l$month),]
118

    
119
xyplot(value~cer|month,data=mod06l,scales=list(relation="free"),pch=16,cex=.5)
120
xyplot(value~cer|station,data=mod06l[mod06l$count>400,],pch=16,cex=.5)
121

    
122
xyplot(cot~month,groups=station,data=mod06l,type="l")
123

    
124

    
125
## explore fit of simple model
126
m=11
127
cor(mod06l[mod06l$month==m,c("value","cer","cld","cot")])
128
lm1=lm(value~latitude+longitude+elevation+cer+cld+cot,data=mod06l[mod06l$month==m,])
129
summary(lm1)
130
crPlots(lm1)
131

    
132
plot(mod06l$value[mod06l$month==m],as.vector(predict(lm1,data=mod06l[mod06l$month==m,])))
133

    
134
### draw some plots
135
gq=function(x,n=10,cut=F) {
136
  if(!cut) return(unique(quantile(x,seq(0,1,len=n+1),na.rm=T)))
137
  if(cut)  return(cut(x,unique(quantile(x,seq(0,1,len=n+1),na.rm=T))))
138
}
139

    
140
### add some additional variables
141
mod06s$month=factor(mod06s$month,labels=format(as.Date(paste("2000",1:12,"15",sep="-")),"%b"))
142
mod06s$lppt=log(mod06s$ppt)
143
mod06s$glon=cut(mod06s$lon,gq(mod06s$lon,n=5),include.lowest=T,ordered=T)#gq(mod06s$lon,n=3))
144
mod06s$glon2=cut(mod06s$lon,breaks=c(-125,-122,-115),labels=c("Coastal","Inland"),include.lowest=T,ordered=T)#gq(mod06s$lon,n=3))
145
mod06s$gelev=cut(mod06s$elev,breaks=gq(mod06s$elev,n=3),labels=c("Low","Mid","High"),include.lowest=T,ordered=T)
146
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)
147
mod06s$LWP_mean=(2/3)*mod06s$CER_mean*mod06s$COT_mean
148

    
149
## melt it
150
mod06sl=melt(mod06s[,!grepl("lppt",colnames(mod06s))],id.vars=c("id","lon","lat","elev","month","ppt","glon","glon2","gelev","gbin"))
151
levels(mod06sl$variable)=c("Effective Radius (um)","Very Cloudy Days (%)","Cloudy Days (%)","Optical Thickness (%)","Liquid Water Path")
152

    
153
###################################################################
154
###################################################################
155

    
156
bgyr=colorRampPalette(c("blue","green","yellow","red"))
157

    
158
X11.options(type="cairo")
159
pdf("output/MOD06_summary.pdf",width=11,height=8.5)
160

    
161
# % cloudy maps
162
title="Cloudiness (% cloudy days) "
163
at=unique(quantile(as.matrix(cld),seq(0,1,len=100),na.rm=T))
164
p=levelplot(cld,xlab.top=title,at=at,col.regions=bgyr(length(at)))+layer(sp.lines(roil, lwd=1.2, col='black'))
165
print(p)
166
#bwplot(cer,main=title,ylab="Cloud Effective Radius (microns)")
167

    
168
# CER maps
169
title="Cloud Effective Radius (microns)"
170
at=quantile(as.matrix(cer),seq(0,1,len=100),na.rm=T)
171
p=levelplot(cer,xlab.top=title,at=at,col.regions=bgyr(length(at)))+layer(sp.lines(roil, lwd=1.2, col='black'))
172
print(p)
173
#bwplot(cer,main=title,ylab="Cloud Effective Radius (microns)")
174

    
175
# COT maps
176
title="Cloud Optical Thickness (%)"
177
at=quantile(as.matrix(cot),seq(0,1,len=100),na.rm=T)
178
p=levelplot(cot,xlab.top=title,at=at,col.regions=bgyr(length(at)))+layer(sp.lines(roil, lwd=0.8, col='black'))
179
print(p)
180
#bwplot(cot,xlab.top=title,ylab="Cloud Optical Thickness (%)")
181

    
182
###########################################
183
#### compare with PRISM data
184

    
185
at=quantile(as.matrix(subset(m01,subset=1)),seq(0,1,len=100),na.rm=T)
186
p1=levelplot(subset(m01,subset=1),xlab.top="Effective Radius (um)",at=at,col.regions=bgyr(length(at)),margin=F,
187
  )+layer(sp.lines(roi_geo, lwd=1.2, col='black'))
188
at=quantile(as.matrix(subset(m01,subset=2)),seq(0,1,len=100),na.rm=T)
189
p2=levelplot(subset(m01,subset=2),xlab.top="Cloudy days (%)",at=at,col.regions=bgyr(length(at)),margin=F,
190
   )+layer(sp.lines(roi_geo, lwd=1.2, col='black'))
191
at=quantile(as.matrix(subset(m01,subset=3)),seq(0,1,len=100),na.rm=T)
192
p3=levelplot(subset(m01,subset=3),xlab.top="Optical Thickness (%)",at=at,col.regions=bgyr(length(at)),margin=F,
193
   )+layer(sp.lines(roi_geo, lwd=1.2, col='black'))
194
at=quantile(as.matrix(subset(m01,subset=4)),seq(0,1,len=100),na.rm=T)
195
p4=levelplot(subset(m01,subset=4),xlab.top="PRISM MAP",at=at,col.regions=bgyr(length(at)),margin=F,
196
    )+layer(sp.lines(roi_geo, lwd=1.2, col='black'))
197

    
198
print(p1,split=c(1,1,2,2))
199
print(p2,split=c(1,2,2,2),new=F)
200
print(p3,split=c(2,1,2,2),new=F)
201
print(p4,split=c(2,2,2,2),new=F)
202

    
203
### compare COT and PRISM
204
print(p3,split=c(1,1,2,1))
205
print(p4,split=c(2,1,2,1),new=F)
206

    
207

    
208
### sample to speed up processing
209
s=sample(1:nrow(bd),10000)
210

    
211
## melt it to ease comparisons
212
bdl=melt(bd[s,],measure.vars=c("cer","cld","cot"))
213

    
214
combineLimits(useOuterStrips(xyplot(prism~value|variable+month,data=bdl,pch=16,cex=.2,scales=list(y=list(log=T),x=list(relation="free")),
215
                                    ylab="PRISM Monthly mean precipitation (mm)",xlab="MOD06 metric",main="PRISM vs. MOD06 (mean monthly ppt)")))+
216
  layer(panel.abline(lm(y~x),col="red"))+
217
  layer(panel.text(0,2.5,paste("R2=",round(summary(lm(y~x))$r.squared,2))))
218

    
219

    
220
### Comparison at station values
221
at=quantile(as.matrix(cotm),seq(0,1,len=100),na.rm=T)
222
p=levelplot(cotm, layers=1,at=at,col.regions=bgyr(length(at)),main="Mean Annual Cloud Optical Thickness",FUN.margin=function(x) 0)+
223
  layer(sp.lines(roil, lwd=1.2, col='black'))+layer(sp.points(st2_sin, pch=16, col='black'))
224
print(p)
225

    
226
### monthly comparisons of variables
227
#mod06sl=melt(mod06s,measure.vars=c("ppt","COT_mean","CER_mean","CER_P20um"))
228
#bwplot(value~month|variable,data=mod06sl,cex=.5,pch=16,col="black",scales=list(y=list(relation="free")),layout=c(1,3))
229
#splom(mod06s[grep("CER|COT|CLD|ppt",colnames(mod06s))],cex=.2,pch=16,main="Scatterplot matrix of MOD06 products")
230

    
231
### run some regressions
232
#plot(log(ppt)~COT_mean,data=mod06s)
233
#summary(lm(log(ppt)~COT_mean*month,data=mod06s))
234

    
235
## ppt~metric with longitude bins
236
 xyplot(ppt~value|variable,groups=glon,data=mod06sl,
237
       scales=list(y=list(log=T),x=list(relation="free",log=F)),
238
       par.settings = list(superpose.symbol = list(col=bgyr(5),pch=16,cex=.5)),auto.key=list(space="top",title="Station Longitude"),
239
       main="Comparison of MOD06_L2 and Precipitation Monthly Climatologies",ylab="Mean Monthly Station Precipitation (mm)",xlab="MOD06_L2 Product",layout=c(5,1))+
240
  layer(panel.text(9,2.5,label="Coastal stations",srt=30,cex=1.3,col="blue"),columns=1)+
241
  layer(panel.text(13,.9,label="Inland stations",srt=10,cex=1.3,col="red"),columns=1)+
242
  layer(panel.abline(lm(y~x),col="red"))+
243
  layer(panel.text(0,0,paste("R2=",round(summary(lm(y~x))$r.squared,2)),pos=4,cex=.5,col="grey"))
244

    
245

    
246
## ppt~metric with longitude bins
247
#CairoPNG("output/COT.png",width=10,height=5,units="in",dpi=300,pointsize=20)
248
#png("output/COT.png",width=10,height=5,units="in",res=150)
249
#trellis.par.set("fontsize",12)
250
at=quantile(as.matrix(subset(m01,subset=3)),seq(0,1,len=100),na.rm=T)
251
p1=levelplot(subset(m01,subset=3),xlab.top="Optical Thickness (%)",at=at,col.regions=bgyr(length(at)),margin=F,
252
   )+layer(sp.lines(roi_geo, lwd=1.2, col='black'))+layer(sp.points(st2, cex=.5,col='black'))
253
at=quantile(as.matrix(subset(m01,subset=3)),seq(0,1,len=100),na.rm=T)
254

    
255
at=quantile(as.matrix(subset(m01,subset=4)),seq(0,1,len=100),na.rm=T)
256
p2=levelplot(subset(m01,subset=4),xlab.top="PRISM MAP",at=at,col.regions=bgyr(length(at)),margin=F,
257
    )+layer(sp.lines(roi_geo, lwd=1.2, col='black'))+layer(sp.points(st2, cex=.5, col='black'))
258

    
259
p3=xyplot(ppt~value,groups=glon,data=mod06sl[mod06sl$variable=="Optical Thickness (%)",],
260
       scales=list(y=list(log=T),x=list(relation="free",log=F)),
261
       par.settings = list(superpose.symbol = list(col=bgyr(5),pch=16,cex=.3)),auto.key=list(space="right",title="Station \n Longitude"),
262
ylab="Mean Monthly Station Precipitation (mm)",xlab="Cloud Optical Thickness from MOD06_L2 (%)",layout=c(1,1))+
263
  layer(panel.text(9,2.6,label="Coastal stations",srt=10,cex=1.3,col="blue"),columns=1)+
264
  layer(panel.text(13,.95,label="Inland stations",srt=10,cex=1.3,col="red"),columns=1)
265

    
266
p4=xyplot(ppt~value,groups=glon,data=mod06sl[mod06sl$variable=="Very Cloudy Days (%)",],
267
       scales=list(y=list(log=T),x=list(relation="free",log=F)),
268
       par.settings = list(superpose.symbol = list(col=bgyr(5),pch=16,cex=.3)),auto.key=list(space="right",title="Station \n Longitude"),
269
ylab="Mean Monthly Station Precipitation (mm)",xlab="Proportion days with Cloud Effective Radius >20um from MOD06_L2 (%)",layout=c(1,1))+
270
  layer(panel.text(9,2.6,label="Coastal stations",srt=10,cex=1.3,col="blue"),columns=1)+
271
  layer(panel.text(13,.95,label="Inland stations",srt=10,cex=1.3,col="red"),columns=1)
272

    
273
#save(p1,p2,p3,file="plotdata.Rdata")
274
#load("plotdata.Rdata")
275

    
276
#CairoPDF("output/MOD06_Summaryfig.pdf",width=11,height=8.5)
277
print(p3,position=c(0,0,1,.5),save.object=F)
278
#print(p4,position=c(0,0,1,.5),save.object=F)
279
print(p1,split=c(1,1,2,2),new=F)
280
print(p2,split=c(2,1,2,2),new=F)
281
#dev.off()
282
#system("convert output/MOD06_Summaryfig.pdf output/MOD06_Summaryfig.png")
283
                                        #dev.off()
284

    
285
## with elevation
286
# xyplot(ppt~value|variable,groups=gbin,data=mod06sl,
287
#       scales=list(y=list(log=T),x=list(relation="free")),
288
#       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"),
289
#       main="Comparison of MOD06_L2 and Precipitation Monthly Climatologies",ylab="Precipitation",xlab="MOD06_L2 Product",layout=c(3,1))+
290
#  layer(panel.text(9,2.5,label="Coastal stations",srt=30,cex=1.3,col="blue"),columns=1)+
291
#  layer(panel.text(13,.9,label="Inland stations",srt=10,cex=1.3,col="red"),columns=1)
292

    
293
## with elevation and longitude bins
294
combineLimits(useOuterStrips(xyplot(ppt~value|variable+gbin,data=mod06sl,
295
       scales=list(y=list(log=T),x=list(relation="free")),col="black",pch=16,cex=.5,type=c("p","r"),
296
       main="Comparison of MOD06_L2 and Precipitation Monthly Climatologies",ylab="Precipitation",xlab="MOD06_L2 Product")))+
297
  layer(panel.xyplot(x,y,type="r",col="red"))
298

    
299
## *** MOD06 vars vs precipitation by month, colored by longitude
300
combineLimits(useOuterStrips(xyplot(ppt~value|month+variable,groups=glon,data=mod06sl,cex=.5,pch=16,
301
                      scales=list(y=list(log=T),x=list(relation="free")),
302
                      par.settings = list(superpose.symbol = list(pch =16, col=bgyr(5),cex=1)),auto.key=list(space="top",title="Station Longitude"),
303
                      main="Comparison of MOD06_L2 and Precipitation Monthly Climatologies",ylab="Precipitation",xlab="MOD06_L2 Product")),
304
              margin.x=1,adjust.labels=F)+
305
  layer(panel.abline(lm(y~x),col="red"))+
306
  layer(panel.text(0,0,paste("R2=",round(summary(lm(y~x))$r.squared,2)),pos=4,cex=.5,col="grey"))
307

    
308

    
309
 xyplot(ppt~CLD_mean|id,data=mod06s,panel=function(x,y,group){
310
  panel.xyplot(x,y,type=c("r"),cex=.5,pch=16,col="red")
311
  panel.xyplot(x,y,type=c("p"),cex=.5,pch=16,col="black")
312
} ,scales=list(y=list(log=T)),strip=F,main="Monthly Mean Precipitation and % Cloudy by station",
313
        sub="Each panel is a station, each point is a monthly mean",
314
        ylab="Precipitation (mm, log axis)",xlab="% of Cloudy Days")+
315
  layer(panel.text(.5,.5,round(summary(lm(y~x))$r.squared,2),pos=4,cex=.75,col="grey"))
316

    
317
 xyplot(ppt~CER_mean|id,data=mod06s,panel=function(x,y,group){
318
  panel.xyplot(x,y,type=c("r"),cex=.5,pch=16,col="red")
319
  panel.xyplot(x,y,type=c("p"),cex=.5,pch=16,col="black")
320
} ,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)")+
321
  layer(panel.text(10,.5,round(summary(lm(y~x))$r.squared,2),pos=4,cex=.75,col="grey"))
322

    
323
 xyplot(ppt~COT_mean|id,data=mod06s,panel=function(x,y,group){
324
  panel.xyplot(x,y,type=c("r"),cex=.5,pch=16,col="red")
325
  panel.xyplot(x,y,type=c("p"),cex=.5,pch=16,col="black")
326
} ,scales=list(y=list(log=T)),strip=F,main="Monthly Mean Precipitation and Cloud Optical Thickness by station",
327
        sub="Each panel is a station, each point is a monthly mean \n Number in lower right of each panel is R^2",
328
        ylab="Precipitation (mm, log axis)",xlab="Mean Monthly Cloud Optical Thickness (%)")+
329
  layer(panel.text(10,.5,round(summary(lm(y~x))$r.squared,2),pos=4,cex=.75,col="grey"))
330

    
331
 xyplot(ppt~LWP_mean|id,data=mod06s,panel=function(x,y,group){
332
  panel.xyplot(x,y,type=c("r"),cex=.5,pch=16,col="red")
333
  panel.xyplot(x,y,type=c("p"),cex=.5,pch=16,col="black")
334
} ,scales=list(y=list(log=T)),strip=F,main="Monthly Mean Precipitation and Liquid Water Path by station",
335
        sub="Each panel is a station, each point is a monthly mean \n Number in lower right of each panel is R^2",
336
        ylab="Precipitation (mm, log axis)",xlab="Mean Monthly Liquid Water Path")+
337
  layer(panel.text(10,.5,round(summary(lm(y~x))$r.squared,2),pos=4,cex=.75,col="grey"))
338

    
339
### Calculate the slope of each line
340
mod06s.sl=dapply(mod06s,list(id=mod06s$id),function(x){
341
  lm1=lm(log(x$ppt)~x$CER_mean,)
342
  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)
343
})
344
mod06s.sl$cex=gq(mod06s.sl$r2,n=5,cut=T)
345
mod06s.sl$cer.s=gq(mod06s.sl$cer,n=5,cut=T)
346

    
347
###  and plot it on a map
348
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,
349
       main="Slopes of linear regressions {log(ppt)~CloudEffectiveRadius}")+
350
  layer(sp.lines(roi_geo, lwd=1.2, col='black'))
351

    
352
### look for relationships with longitude
353
xyplot(cer~lon,group=cut(mod06s.sl$elev,gq(mod06s.sl$elev,n=5)),data=mod06s.sl,
354
       par.settings = list(superpose.symbol = list(col=bgyr(5),pch=16,cex=1)),auto.key=list(space="right",title="Station Elevation"),
355
       ylab="Slope of lm(ppt~EffectiveRadius)",xlab="Longitude",main="Precipitation~Effective Radius relationship by latitude")
356

    
357

    
358
############################################################
359
### simple regression to get spatial residuals
360
m="01"
361
mod06s2=mod06s#[mod06s$month==m,]
362

    
363
lm1=lm(log(ppt)~CER_mean*month*lon,data=mod06s2); summary(lm1)
364
mod06s2$pred=exp(predict(lm1,mod06s2))
365
mod06s2$resid=mod06s2$pred-mod06s2$ppt
366
mod06s2$residg=gq(mod06s2$resid,n=5,cut=T)
367
mod06s2$presid=mod06s2$resid/mod06s2$ppt
368

    
369
for(l in c(F,T)){
370
## all months
371
  xyplot(pred~ppt,groups=gelev,data=mod06s2,
372
       par.settings = list(superpose.symbol = list(col=bgyr(3),pch=16,cex=.75)),auto.key=list(space="right",title="Station Elevation"),
373
       scales=list(log=l),
374
       ylab="Predicted Mean Monthly Precipitation (mm)",xlab="Observed Mean Monthly Precipitation (mm)",main="Predicted vs. Observed for Simple Model",
375
       sub="Red line is y=x")+
376
  layer(panel.abline(0,1,col="red"))
377

    
378
## month by month
379
  print(xyplot(pred~ppt|month,groups=gelev,data=mod06s2,
380
       par.settings = list(superpose.symbol = list(col=bgyr(3),pch=16,cex=.75)),auto.key=list(space="right",title="Station Elevation"),
381
       scales=list(log=l),
382
       ylab="Predicted Mean Monthly Precipitation (mm)",xlab="Observed Mean Monthly Precipitation (mm)",main="Predicted vs. Observed for Simple Model",
383
       sub="Red line is y=x")+
384
  layer(panel.abline(0,1,col="red"))
385
)}
386

    
387
## residuals by month
388
xyplot(lat~lon|month,group=residg,data=mod06s2,
389
       par.settings = list(superpose.symbol = list(pch =16, col=bgyr(5),cex=.5)),
390
       auto.key=list(space="right",title="Residuals"),
391
       main="Spatial plot of monthly residuals")+
392
    layer(sp.lines(roi_geo, lwd=1.2, col='black'))
393

    
394

    
395
dev.off()
396

    
397

    
398
####################################
399
#### build table comparing various metrics
400
mods=data.frame(
401
  models=c(
402
    "log(ppt)~CER_mean",
403
    "log(ppt)~CLD_mean",
404
    "log(ppt)~COT_mean",
405
    "log(ppt)~CER_mean*month",
406
    "log(ppt)~CLD_mean*month",
407
    "log(ppt)~COT_mean*month",
408
    "log(ppt)~CER_mean*month*lon",
409
    "log(ppt)~CLD_mean*month*lon",
410
    "log(ppt)~COT_mean*month*lon",
411
    "ppt~CER_mean*month*lon",
412
    "ppt~CLD_mean*month*lon",
413
    "ppt~COT_mean*month*lon"),stringsAsFactors=F)
414
  
415
mods$r2=
416
  do.call(rbind,lapply(1:nrow(mods),function(i){
417
    lm1=lm(as.formula(mods$models[i]),data=mod06s2)
418
    summary(lm1)$r.squared}))
419

    
420
mods
421

    
422

    
423

    
424

    
425

    
426

    
427

    
428

    
429
load("data/modis/pointsummary.Rdata")
430

    
431

    
432
dsl=melt(ds,id.vars=c("id","date","ppt","lon","lat"),measure.vars=  c("Cloud_Water_Path","Cloud_Effective_Radius","Cloud_Optical_Thickness"))
433

    
434
dsl=dsl[!is.nan(dsl$value),]
435

    
436

    
437

    
438

    
439
####
440
## mean annual precip
441
dp=d[d$variable=="ppt",]
442
dp$year=format(dp$date,"%Y")
443
dm=tapply(dp$value,list(id=dp$id,year=dp$year),sum,na.rm=T)
444
dms=apply(dm,1,mean,na.rm=T)
445
dms=data.frame(id=names(dms),ppt=dms/10)
446

    
447
dslm=tapply(dsl$value,list(id=dsl$id,variable=dsl$variable),mean,na.rm=T)
448
dslm=data.frame(id=rownames(dslm),dslm)
449

    
450
dms=merge(dms,dslm)
451
dmsl=melt(dms,id.vars=c("id","ppt"))
452

    
453
summary(lm(ppt~Cloud_Effective_Radius,data=dms))
454
summary(lm(ppt~Cloud_Water_Path,data=dms))
455
summary(lm(ppt~Cloud_Optical_Thickness,data=dms))
456
summary(lm(ppt~Cloud_Effective_Radius+Cloud_Water_Path+Cloud_Optical_Thickness,data=dms))
457

    
458

    
459
#### draw some plots
460
#pdf("output/MOD06_summary.pdf",width=11,height=8.5)
461
png("output/MOD06_summary_%d.png",width=1024,height=780)
462

    
463
 ## daily data
464
xyplot(value~ppt/10|variable,data=dsl,
465
       scales=list(relation="free"),type=c("p","r"),
466
       pch=16,cex=.5,layout=c(3,1))
467

    
468

    
469
densityplot(~value|variable,groups=cut(dsl$ppt,c(0,50,100,500)),data=dsl,auto.key=T,
470
            scales=list(relation="free"),plot.points=F)
471

    
472
## annual means
473

    
474
xyplot(value~ppt|variable,data=dmsl,
475
       scales=list(relation="free"),type=c("p","r"),pch=16,cex=0.5,layout=c(3,1),
476
       xlab="Mean Annual Precipitation (mm)",ylab="Mean value")
477

    
478
densityplot(~value|variable,groups=cut(dsl$ppt,c(0,50,100,500)),data=dmsl,auto.key=T,
479
            scales=list(relation="free"),plot.points=F)
480

    
481

    
482
## plot some swaths
483

    
484
nc1=raster(fs$path[3],varname="Cloud_Effective_Radius")
485
nc2=raster(fs$path[4],varname="Cloud_Effective_Radius")
486
nc3=raster(fs$path[5],varname="Cloud_Effective_Radius")
487

    
488
nc1[nc1<=0]=NA
489
nc2[nc2<=0]=NA
490
nc3[nc3<=0]=NA
491

    
492
plot(roi)
493
plot(nc3)
494

    
495
plot(nc1,add=T)
496
plot(nc2,add=T)
497

    
498

    
499
dev.off()
500

    
501

    
(11-11/14)