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

###  R code to aquire and process MOD06_L2 cloud data from the MODIS platform

## connect to server of choice

#system("ssh litoria")

#R

library(sp)

library(spgrass6)

library(rgdal)

library(geosphere)

library(rgeos)

library(multicore)

library(raster)

library(lattice)

library(rgl)

library(hdf5)

X11.options(type="Xlib")

ncores=20  #number of threads to use

### Produce shapefile of all stations globally

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="")

colnames(ghcn)=c("id","lat","lon","dem","state","name","gsnflag","hcnflag","wmoid")

ghcn=ghcn[,c("id","dem","lat","lon")]

ghcn$id=as.character(ghcn$id)

ghcn$type="GHCN"

### FAO station normals

fao=read.csv("data/stationdata/FAOCLIM2/world.csv")

fao=fao[,c("ID","Elevation","Lat","Lon")]

colnames(fao)=c("id","dem","lat","lon")

fao$type="FAO"

st=rbind.data.frame(ghcn,fao)

coordinates(st)=c("lon","lat")

projection(st)=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs ")

writeOGR(st,dsn="data/",layer="stationlocations",driver="ESRI Shapefile")

plot(fao,pch=16,cex=.2,col="blue")

plot(st,pch=16,cex=.2,col="red",add=T)

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

#### explore the data

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

## load data

cer=brick(stack(cerfiles))

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

cer@z=list(months)

cer@zname="time"

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

#cer=projectRaster(from=cer,crs="+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0",method="ngb")

### TODO: change to bilinear!

cotfiles=list.files(summarydatadir,pattern="COT_mean_.*tif$",full=T); cotfiles

cot=brick(stack(cotfiles))

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

cot@z=list(months)

cot@zname="time"

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

#cot=projectRaster(from=cot,crs="+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0",method="ngb")

cotm=mean(cot,na.rm=T)

### TODO: change to bilinear!

cld=brick(stack(cldfiles))

cld[cld==0]=NA

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

cld@z=list(months)

cld@zname="time"

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

#cot=projectRaster(from=cot,crs="+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0",method="ngb")

cldm=mean(cld,na.rm=T)

cer20=brick(stack(cer20files))

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

cer20@z=list(months)

cer20@zname="time"

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

#cot=projectRaster(from=cot,crs="+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0",method="ngb")

cotm=mean(cot,na.rm=T)

### TODO: change to bilinear!

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

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))

load("data/allstations.Rdata")

st2=st[st$id%in%d$id,]

### extract MOD06 data for each station

stcer=extract(cer,st2)#;colnames(stcer)=paste("cer_mean_",1:12,sep="")

mod06l[,c("variable","moment","month")]=do.call(rbind,strsplit(as.character(mod06l$variable),"_"))

### Identify stations that have < 10 years of data

cnts=cast(d2,id~.,fun=function(x) length(x[!is.na(x)]),value="count");colnames(cnts)[colnames(cnts)=="(all)"]="count"

summary(cnts)

## drop them

d2=d2[d2$id%in%cnts$id[cnts$count>=365*10],]

### generate monthly means of station data

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

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

mod06s=merge(dcl,mod06l)

### draw some plots

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

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

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

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

title="Cloudiness (% cloudy days) "

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

p=levelplot(cld,xlab.top=title,at=at,col.regions=bgyr(length(at)))+layer(sp.lines(roil, lwd=1.2, col='black'))

print(p)

title="Cloud Effective Radius (microns)"

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

p=levelplot(cer,xlab.top=title,at=at,col.regions=bgyr(length(at)))+layer(sp.lines(roil, lwd=1.2, col='black'))

print(p)

# COT maps

title="Cloud Optical Thickness (%)"

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

p=levelplot(cot,xlab.top=title,at=at,col.regions=bgyr(length(at)))+layer(sp.lines(roil, lwd=0.8, col='black'))

print(p)

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

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

### Comparison at station values

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

p=levelplot(cotm, layers=1,at=at,col.regions=bgyr(length(at)),main="Mean Annual Cloud Optical Thickness",FUN.margin=function(x) 0)+

### monthly comparisons of variables

#bwplot(value~month|variable,data=mod06sl,cex=.5,pch=16,col="black",scales=list(y=list(relation="free")),layout=c(1,3))

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

### run some regressions

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

## ppt~metric with longitude bins

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

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

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

#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'))

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)

       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="Proportion days with Cloud Effective Radius >20um 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")

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

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

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

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

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

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

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

### Calculate the slope of each line

  lm1=lm(log(x$ppt)~x$CER_mean,)

mod06s.sl$cex=gq(mod06s.sl$r2,n=5,cut=T)

mod06s.sl$cer.s=gq(mod06s.sl$cer,n=5,cut=T)

       main="Slopes of linear regressions {log(ppt)~CloudEffectiveRadius}")+

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

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

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

### simple regression to get spatial residuals

m="01"

mod06s2=mod06s#[mod06s$month==m,]

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

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,