### set up grass session

  tf=paste(tempdir(),"/grass", Sys.getpid(),"/", sep="")

  ## set up tempfile for this PID

  initGRASS(gisBase="/usr/lib/grass64",gisDbase=tf,SG=td,override=T,location="mod06",mapset="PERMANENT",home=tf,pid=Sys.getpid())

  system(paste("g.proj -c proj4=\"+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +datum=WGS84 +units=m +no_defs\"",sep=""))

#system("g.mapset PERMANENT")

  execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",outdir,"/",fs$file[1],"\":mod06:Cloud_Mask_1km_0",sep=""),

            output="modisgrid",flags=c("quiet","overwrite","o"))

  system("g.region rast=modisgrid save=roi --overwrite")

  system("g.region roi")

  system("g.region -p")

#  getLocationProj()

  ## Identify which files to process

  ### print some summary info

#                 QA_CWP_",i,"=   ((QA_",i," / 2^8) % 2^1 )==1

#                 QA_CWP2_",i,"=  ((QA_",i," / 2^9) % 2^2 )==3

#   execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod06:Cloud_Water_Path",sep=""),

#            output=paste("CWP_",i,sep=""),title="cloud_water_path",

#            flags=c("overwrite","o")) ; print("")

#   execGRASS("r.null",map=paste("CWP_",i,sep=""),setnull="-9999")

#   system(paste("r.mapcalc \"CWP_",i,"=if(QA_CWP_",i,"&&QA_CWP2_",i,"&&CWP_",i,">=0,CWP_",i,"*0.009999999776482582,null())\"",sep=""))   

#   system(paste("r.mapcalc \"CWP2_",i,"=if(CM_clear_",i,"==0,CWP_",i,",0)\"",sep=""))   

  execGRASS("r.out.gdal",input="CER_daily",output=paste(tifdir,"/CER_",format(date,"%Y%m%d"),".tif",sep=""),nodata=-999,flags=c("quiet"))

  execGRASS("r.out.gdal",input="COT_daily",output=paste(tifdir,"/COT_",format(date,"%Y%m%d"),".tif",sep=""),nodata=-999,flags=c("quiet"))

  execGRASS("r.out.gdal",input="CLD_daily",output=paste(tifdir,"/CLD_",format(date,"%Y%m%d"),".tif",sep=""),nodata=99,flags=c("quiet"))

 system(paste("rm -R ",tf,sep=""))

### print update

  print(paste(" ###################################################################               Finished ",date,"

################################################################"))

done=tdates%in%as.Date(substr(list.files(tifdir),5,12),"%Y%m%d")

mclapply(tdates,function(date) loadcloud(date,fs=fs))

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

library(reshape2)

library(heR.Misc)

roi_geo=as(roi,"SpatialLines")

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

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)

### TODO: change to bilinear!

### get station data, subset to stations in region, and transform to sinusoidal

d2=d[d$variable=="ppt"&d$date>=as.Date("2000-01-01"),]

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

d2$value=d2$value/10 #convert to mm

stcld=extract(cld,st2)#;colnames(stcld)=paste("cld_mean_",1:12,sep="")

mod06=cbind.data.frame(id=st2$id,lat=st2$lat,lon=st2$lon,stcer,stcot,stcld)

mod06l=as.data.frame(cast(mod06l,id+lon+lat+month~variable+moment,value="value"))

### 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],]

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

gq=function(x,n=10,cut=F) {

  if(!cut) unique(quantile(x,seq(0,1,len=n+1)))

  if(cut)  cut(x,unique(quantile(x,seq(0,1,len=n+1))))

}

# % cloudy maps

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)

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

splom(mod06s[grep("CER|COT|CLD",colnames(mod06s))],cex=.2,pch=16)

### run some regressions

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

xyplot(ppt~CLD_mean,data=mod06s,cex=.5,pch=16,col="black",scales=list(y=list(log=F)),main="Comparison of monthly mean CLD and precipitation",ylab="Precipitation (log axis)",xlab="% Days Cloudy")

xyplot(ppt~CER_mean,data=mod06s,cex=.5,pch=16,col="black",scales=list(y=list(log=T)),main="Comparison of monthly mean CER and precipitation",ylab="Precipitation (log axis)")

xyplot(ppt~COT_mean,data=mod06s,cex=.5,pch=16,col="black",scales=list(y=list(log=T)),main="Comparison of monthly mean COT and precipitation",ylab="Precipitation (log axis)")

xyplot(ppt~CER_mean|month,data=mod06s,cex=.5,pch=16,col="black",scales=list(log=T,relation="free"))

xyplot(ppt~COT_mean|month,data=mod06s,cex=.5,pch=16,col="black",scales=list(log=T,relation="free"))

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

} ,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",ylab="Precipitation (mm, log axis)",xlab="Mean Monthly Cloud Optical Thickness")

### Calculate the slope of each line and plot it on a map

mod06s.sl=dapply(mod06s,list(id=mod06s$id),function(x){

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

  data.frame(lat=x$lat[1],lon=x$lon[1],intcpt=coefficients(lm1)[1],cer=coefficients(lm1)[2],r2=summary(lm1)$r.squared)

})

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

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

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,

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

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

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

### simple regression to get spatial residuals

m="01"

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

lm1=lm(ppt~CER_mean*month,data=mod06s2)

summary(lm1)

mod06s2$pred=predict(lm1,mod06s2)

mod06s2$resid=mod06s2$pred-mod06s2$ppt

mod06sr=cast(mod06s2,id+lon+lat~month,value="resid",fun=function(x) mean(x,na.rm=T))

mod06sr=melt(mod06sr,id.vars=c("id","lon","lat"),value="resid")

mod06sr$residg=cut(mod06sr$value,quantile(mod06sr$value,seq(0,1,len=11),na.rm=T))

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

         par.settings = list(superpose.symbol = list(pch =16, col=terrain.colors(10),cex=.5)),

           auto.key=T)

plot(pred~value,data=mod06s,log="xy")