## specify some working directories

tifdir="/media/data/MOD06_L2_tif"

summarydatadir="data/modis/MOD06_climatologies"

  byrow=T,ncol=2,dimnames=list(1:10,c("variable","varid"))),stringsAsFactors=F)

swath2grid=function(i=1,files,vars,outdir,upleft="47 -125",lowright="41 -115"){

#RESAMPLING_TYPE =",ifelse(grepl("Flag|Mask|Quality",vars),"NN","CUBIC"),"

projection(td)="+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs +datum=WGS84 +ellps=WGS84 "

## temporary objects to test function below

### Function to extract various SDSs from a single gridded HDF file and use QA data to throw out 'bad' observations

    ## Identify which files to process

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

   ## keep only positive CWP values where quality is 'useful' and 'very good' & scale to real units

   ## keep only positive CER values where quality is 'useful' and 'very good' & scale to real units

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

   ## set CER to 0 in clear-sky pixels

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

#### Now generate daily averages (or maximum in case of cloud flag)

         CLD_daily=max(",paste("if(isnull(CM_cloud_",1:nfs,"),0,CM_cloud_",1:nfs,")",sep="",collapse=","),") 

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

tdates=sort(unique(fs$date))

  path=list.files(tifdir,full=T,recursive=T,pattern="tif$"),

  file=basename(list.files(tifdir,full=F,recursive=T,pattern="tif$")))

fs2$type=substr(fs2$file,1,3)

fs2$date=as.Date(substr(fs2$file,5,12),"%Y%m%d")

fs2$month=format(fs2$date,"%m")

fs2$year=format(fs2$date,"%Y")

fs2$path=as.character(fs2$path)

fs2$file=as.character(fs2$file)

# Define type/month products

vs=expand.grid(type=unique(fs2$type),month=c("01","02","03","04","05","06","07","08","09","10","11","12"))

## identify which have been completed

#done=

#  do.call(rbind,strsplit(list.files(summarydatadir),"_|[.]"))[,3]

#table(done)

#tdates=tdates[!done]

## process the summaries using the raster package

mclapply(1:nrow(vs),function(i){

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

  td=stack(fs2$path[which(fs2$month==vs$month[i]&fs2$type==vs$type[i])])

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

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

       format="GTiff")

  calc(td,sd,na.rm=T,

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

       format="GTiff")

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

  calc(td,function(i)

       sum(!is.na(i)),filename=paste(summarydatadir,"/",vs$type[i],"_count_",vs$month[i],".tif",sep=""),

       format="GTiff")

  calc(td,function(i) sum(ifelse(i==0,0,1)),

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

)

library(rasterVis)

roil=as(roi,"SpatialLines")

summarydatadir="data/modis/MOD06_climatologies"

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

#### explore the data

## load data

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

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

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!

d2=d[d$variable=="ppt",]

st2=spTransform(st2,CRS(" +proj=sinu +lon_0=0 +x_0=0 +y_0=0"))

d2$month=format(d2$date,"%m")

### extract MOD06 data for each station

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

stcot=extract(cot,st2)#;colnames(stcot)=paste("cot_mean_",1:12,sep="")

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

mod06l=melt(mod06,id.vars=c("id","lon","lat"))

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

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

### generate monthly means of station data

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

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

## merge station data with mod06

mod06s=merge(dcl,mod06l)

mod06s$lvalue=log(mod06s$value+1)

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

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

### draw some plots

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

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

# CER maps

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)

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

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

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

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

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

print(p)

### monthly comparisons of variables

mod06sl=melt(mod06s,measure.vars=c("value","COT_mean","CER_mean"))

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",colnames(mod06s))],cex=.2,pch=16)

xyplot(value~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(value~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(value~CER_mean|month,data=mod06s,cex=.5,pch=16,col="black",scales=list(log=T,relation="free"))

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

#xyplot(value~CER_mean|month+cut(COT_mean,breaks=c(0,10,20)),data=mod06s,cex=.5,pch=16,col="black")#,scales=list(relation="fixed"))

summary(lm(lvalue~COT_mean*month,data=mod06s))

dev.off()