### Brief script to explore the total size of the 1km daily (1970-2010) dataset

## Time

dates=seq(as.Date("1970-01-01"),as.Date("2010-01-01"),by=1)

n.dates=length(dates)

## Total area

n.globe=510072000 #km2, from http://en.wikipedia.org/wiki/Earth

## Land Area

n.land=148940000 #km2, from http://en.wikipedia.org/wiki/Earth

## total values

n.landtime=n.dates*n.land

n.globetime=n.dates*n.globe

n.landtime

n.globetime

formatC(n.landtime, format="f",digits=0, big.mark=',')

### approximate storage size for various data types

types=as.data.frame(matrix(c(

  "Short integer",-32768, 32767,2,

  "Long integer", -2147483648,2147483647,4,

  "Single-precision floating-point", -3.4e38,1.2e38,4,

  "Double-precision floating-point", -2.2e308,1.8e308,8

                         ),ncol=4,byrow=T,dimnames=list(1:4,c("type","min","max","bytes"))),stringsAsFactors=F)

types$min=as.numeric(types$min);types$max=as.numeric(types$max);types$bytes=as.numeric(types$bytes)

types

### now estimate TB storage for total values (not full grid, just the values)

types$storage.landtime.TB=n.landtime*types$bytes/2^40

types$storage.globetime.TB=n.globetime*types$bytes/2^40; types

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

    ## also produce means that first eliminate 0 values (added above to indicate clear skies)

    ## to capture 'when cloudy, how thick are the clouds' rather than mean thickness...

    td2=td

    td2[td2==0]=NA

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

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

         format="GTiff")

  }  

library(reshape)

### TODO: change to bilinear if reprojecting!

### load PRISM data for comparison

prism=brick("data/prism/prism_climate.nc",varname="ppt")

## project to sinusoidal

projection(prism)=CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")

prism=projectRaster(prism,cer)

prism[prism<0]=NA  #for some reason NAvalue() wasn't working

setZ(prism,months,name="time")

prism@z=list(months)

prism@zname="time"

layerNames(prism) <- as.character(format(months,"%b"))

####  build a pixel by variable matrix

vars=c("cer","cld","cot","prism")

bd=melt(as.matrix(vars[1]))

colnames(bd)=c("cell","month",vars[1])

for(v in vars[-1]) {print(v); bd[,v]=melt(as.matrix(get(v)))$value}

bd=bd[!is.na(bd$cer)|is.na(bd$prism),]

## Summarize annual metrics for full rasters

### get all variables from all months

c01=brick(mclapply(vars,function(v) projectRaster(get(v),crs="+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")))

layerNames(m01)=paste(vars,months,sep="_")

m01=brick(mclapply(vars,function(v) mean(get(v))))#mean(cer),mean(cld),mean(cot),mean(prism))

layerNames(m01)=vars

m01=projectRaster(from=m01,crs="+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")

m01=crop(m01,extent(-125,-115,41,47))

st2_sin=spTransform(st2,CRS(projection(cer)))

d2$elev=st2$elev[match(d2$id,st2$id)]

dc=cast(d2,id+lon+lat+elev~month,value="value",fun=function(x) mean(x,na.rm=T)*30)

dcl=melt(dc,id.vars=c("id","lon","lat","elev"),value="ppt")

colnames(dcl)[colnames(dcl)=="value"]="ppt"

  if(!cut) return(unique(quantile(x,seq(0,1,len=n+1),na.rm=T)))

  if(cut)  return(cut(x,unique(quantile(x,seq(0,1,len=n+1),na.rm=T))))

### add some additional variables

mod06s$month=factor(mod06s$month,labels=format(as.Date(paste("2000",1:12,"15",sep="-")),"%b"))

mod06s$lppt=log(mod06s$ppt)

mod06s$glon=cut(mod06s$lon,gq(mod06s$lon,n=5),include.lowest=T,ordered=T)#gq(mod06s$lon,n=3))

mod06s$glon2=cut(mod06s$lon,breaks=c(-125,-122,-115),labels=c("Coastal","Inland"),include.lowest=T,ordered=T)#gq(mod06s$lon,n=3))

mod06s$gelev=cut(mod06s$elev,breaks=gq(mod06s$elev,n=3),labels=c("Low","Mid","High"),include.lowest=T,ordered=T)

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)

## melt it

mod06sl=melt(mod06s[,!grepl("lppt",colnames(mod06s))],id.vars=c("id","lon","lat","elev","month","ppt","glon","glon2","gelev","gbin"))

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

###################################################################

###################################################################

#bwplot(cer,main=title,ylab="Cloud Effective Radius (microns)")

#bwplot(cer,main=title,ylab="Cloud Effective Radius (microns)")

#bwplot(cot,xlab.top=title,ylab="Cloud Optical Thickness (%)")

###########################################

#### compare with PRISM data

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

p1=levelplot(subset(m01,subset=1),xlab.top="Effective Radius (um)",at=at,col.regions=bgyr(length(at)),margin=F,

  )+layer(sp.lines(roi_geo, lwd=1.2, col='black'))

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

p2=levelplot(subset(m01,subset=2),xlab.top="Cloudy days (%)",at=at,col.regions=bgyr(length(at)),margin=F,

   )+layer(sp.lines(roi_geo, lwd=1.2, col='black'))

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

p3=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'))

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

p4=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'))

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

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

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

print(p4,split=c(2,2,2,2),new=F)

### compare COT and PRISM

print(p3,split=c(1,1,2,1))

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

### sample to speed up processing

s=sample(1:nrow(bd),10000)

## melt it to ease comparisons

bdl=melt(bd[s,],measure.vars=c("cer","cld","cot"))

combineLimits(useOuterStrips(xyplot(prism~value|variable+month,data=bdl,pch=16,cex=.2,scales=list(y=list(log=T),x=list(relation="free")),

                                    ylab="PRISM Monthly mean precipitation (mm)",xlab="MOD06 metric",main="PRISM vs. MOD06 (mean monthly ppt)")))+

  layer(panel.abline(lm(y~x),col="red"))+

  layer(panel.text(0,2.5,paste("R2=",round(summary(lm(y~x))$r.squared,2))))

  layer(sp.lines(roil, lwd=1.2, col='black'))+layer(sp.points(st2_sin, pch=16, col='black'))

splom(mod06s[grep("CER|COT|CLD",colnames(mod06s))],cex=.2,pch=16,main="Scatterplot matrix of MOD06 products")

#plot(log(ppt)~COT_mean,data=mod06s)

#summary(lm(log(ppt)~COT_mean*month,data=mod06s))

## ppt~metric with longitude bins

 xyplot(ppt~value|variable,groups=glon,data=mod06sl,

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

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

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

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

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

## with elevation

# xyplot(ppt~value|variable,groups=gbin,data=mod06sl,

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

#       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"),

#       main="Comparison of MOD06_L2 and Precipitation Monthly Climatologies",ylab="Precipitation",xlab="MOD06_L2 Product",layout=c(3,1))+

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

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

## with elevation and longitude bins

combineLimits(useOuterStrips(xyplot(ppt~value|variable+gbin,data=mod06sl,

       scales=list(y=list(log=T),x=list(relation="free")),col="black",pch=16,cex=.5,type=c("p","r"),

       main="Comparison of MOD06_L2 and Precipitation Monthly Climatologies",ylab="Precipitation",xlab="MOD06_L2 Product")))+

  layer(panel.xyplot(x,y,type="r",col="red"))

## *** MOD06 vars vs precipitation by month, colored by longitude

combineLimits(useOuterStrips(xyplot(ppt~value|month+variable,groups=glon,data=mod06sl,cex=.5,pch=16,

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

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

                      main="Comparison of MOD06_L2 and Precipitation Monthly Climatologies",ylab="Precipitation",xlab="MOD06_L2 Product")),

              margin.x=1,adjust.labels=F)+

  layer(panel.abline(lm(y~x),col="red"))+

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

 xyplot(ppt~CLD_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 % Cloudy by station",

        sub="Each panel is a station, each point is a monthly mean",

        ylab="Precipitation (mm, log axis)",xlab="% of Cloudy Days")+

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

 xyplot(ppt~CER_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 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)")+

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

} ,scales=list(y=list(log=T)),strip=F,main="Monthly Mean Precipitation and Cloud Optical Thickness 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 Cloud Optical Thickness (%)")+

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

### Calculate the slope of each line

  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)

###  and plot it on a map

### look for relationships with longitude

xyplot(cer~lon,group=cut(mod06s.sl$elev,gq(mod06s.sl$elev,n=5)),data=mod06s.sl,

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

       ylab="Slope of lm(ppt~EffectiveRadius)",xlab="Longitude",main="Precipitation~Effective Radius relationship by latitude")

lm1=lm(log(ppt)~CER_mean*month*lon,data=mod06s2); summary(lm1)

mod06s2$pred=exp(predict(lm1,mod06s2))

mod06s2$residg=gq(mod06s2$resid,n=5,cut=T)

mod06s2$presid=mod06s2$resid/mod06s2$ppt

for(l in c(F,T)){

## all months

  xyplot(pred~ppt,groups=gelev,data=mod06s2,

       par.settings = list(superpose.symbol = list(col=bgyr(3),pch=16,cex=.75)),auto.key=list(space="right",title="Station Elevation"),

       scales=list(log=l),

       ylab="Predicted Mean Monthly Precipitation (mm)",xlab="Observed Mean Monthly Precipitation (mm)",main="Predicted vs. Observed for Simple Model",

       sub="Red line is y=x")+

  layer(panel.abline(0,1,col="red"))

## month by month

  print(xyplot(pred~ppt|month,groups=gelev,data=mod06s2,

       par.settings = list(superpose.symbol = list(col=bgyr(3),pch=16,cex=.75)),auto.key=list(space="right",title="Station Elevation"),

       scales=list(log=l),

       ylab="Predicted Mean Monthly Precipitation (mm)",xlab="Observed Mean Monthly Precipitation (mm)",main="Predicted vs. Observed for Simple Model",

       sub="Red line is y=x")+

  layer(panel.abline(0,1,col="red"))

)}

## residuals by month

xyplot(lat~lon|month,group=residg,data=mod06s2,

       par.settings = list(superpose.symbol = list(pch =16, col=bgyr(5),cex=.5)),

       auto.key=list(space="right",title="Residuals"),

       main="Spatial plot of monthly residuals")+

    layer(sp.lines(roi_geo, lwd=1.2, col='black'))

####################################

#### build table comparing various metrics

mods=data.frame(

  models=c(

    "log(ppt)~CER_mean",

    "log(ppt)~CLD_mean",

    "log(ppt)~COT_mean",

    "log(ppt)~CER_mean*month",

    "log(ppt)~CLD_mean*month",

    "log(ppt)~COT_mean*month",

    "log(ppt)~CER_mean*month*lon",

    "log(ppt)~CLD_mean*month*lon",

    "log(ppt)~COT_mean*month*lon",

    "ppt~CER_mean*month*lon",

    "ppt~CLD_mean*month*lon",

    "ppt~COT_mean*month*lon"),stringsAsFactors=F)

mods$r2=

  do.call(rbind,lapply(1:nrow(mods),function(i){

    lm1=lm(as.formula(mods$models[i]),data=mod06s2)

    summary(lm1)$r.squared}))

mods