*[~]

*[#]*

### set options for Raster Package

setOptions(progress="text",timer=T)

  print(paste("Loading ",vs$type[i]," for month ",vs$month[i]))

  if(vs$type[i]%in%c("CER")) {

    ## Calculate number of days with effective radius > 20um, found to be linked to precipitating clouds

    ## (Kobayashi 2007)

    calc(td,function(x) mean(ifelse(x<20,0,1),na.rm=T),

      filename=paste(summarydatadir,"/",vs$type[i],"_P20um_",vs$month[i],".tif",sep=""),

      format="GTiff")

#    td2[td2<20]=0

#    td2[td2>=20]=1

#    calc(td2,mean,na.rm=T,

#         filename=paste(summarydatadir,"/",vs$type[i],"_P20um_",vs$month[i],".tif",sep=""),

#         format="GTiff")

mod06s$LWP_mean=(2/3)*mod06s$CER_mean*mod06s$COT_mean

levels(mod06sl$variable)=c("Effective Radius (um)","Cloudy Days (%)","Optical Thickness (%)","Liquid Water Path")

X11.options(type="cairo")

       scales=list(y=list(log=T),x=list(relation="free",log=F)),

       par.settings = list(superpose.symbol = list(col=bgyr(5),pch=16,cex=.5)),auto.key=list(space="top",title="Station Longitude"),

       main="Comparison of MOD06_L2 and Precipitation Monthly Climatologies",ylab="Mean Monthly Station Precipitation (mm)",xlab="MOD06_L2 Product",layout=c(4,1))+

## ppt~metric with longitude bins

#CairoPNG("output/COT.png",width=10,height=5,units="in",dpi=300,pointsize=20)

#png("output/COT.png",width=10,height=5,units="in",res=150)

#trellis.par.set("fontsize",12)

at=quantile(as.matrix(subset(m01,subset=3)),seq(0,1,len=100),na.rm=T)

p1=levelplot(subset(m01,subset=3),xlab.top="Optical Thickness (%)",at=at,col.regions=bgyr(length(at)),margin=F,

   )+layer(sp.lines(roi_geo, lwd=1.2, col='black'))+layer(sp.points(st2, cex=.5,col='black'))

at=quantile(as.matrix(subset(m01,subset=4)),seq(0,1,len=100),na.rm=T)

p2=levelplot(subset(m01,subset=4),xlab.top="PRISM MAP",at=at,col.regions=bgyr(length(at)),margin=F,

    )+layer(sp.lines(roi_geo, lwd=1.2, col='black'))+layer(sp.points(st2, cex=.5, col='black'))

p3=xyplot(ppt~value,groups=glon,data=mod06sl[mod06sl$variable=="Optical Thickness (%)",],

       scales=list(y=list(log=T),x=list(relation="free",log=F)),

       par.settings = list(superpose.symbol = list(col=bgyr(5),pch=16,cex=.3)),auto.key=list(space="right",title="Station \n Longitude"),

ylab="Mean Monthly Station Precipitation (mm)",xlab="Cloud Optical Thickness from MOD06_L2 (%)",layout=c(1,1))+

  layer(panel.text(9,2.6,label="Coastal stations",srt=10,cex=1.3,col="blue"),columns=1)+

  layer(panel.text(13,.95,label="Inland stations",srt=10,cex=1.3,col="red"),columns=1)

save(p1,p2,p3,file="plotdata.Rdata")

load("plotdata.Rdata")

CairoPDF("output/MOD06_Summaryfig.pdf",width=11,height=8.5)

print(p3,position=c(0,0,1,.5),save.object=F)

print(p1,split=c(1,1,2,2),new=F)

print(p2,split=c(2,1,2,2),new=F)

dev.off()

system("convert output/MOD06_Summaryfig.pdf output/MOD06_Summaryfig.png")

                                        #dev.off()

 xyplot(ppt~LWP_mean|id,data=mod06s,panel=function(x,y,group){

  panel.xyplot(x,y,type=c("r"),cex=.5,pch=16,col="red")

  panel.xyplot(x,y,type=c("p"),cex=.5,pch=16,col="black")

} ,scales=list(y=list(log=T)),strip=F,main="Monthly Mean Precipitation and Liquid Water Path by station",

        sub="Each panel is a station, each point is a monthly mean \n Number in lower right of each panel is R^2",

        ylab="Precipitation (mm, log axis)",xlab="Mean Monthly Liquid Water Path")+

  layer(panel.text(10,.5,round(summary(lm(y~x))$r.squared,2),pos=4,cex=.75,col="grey"))