library(raster)

#require(plotKML); data(worldgrids_pal)  #load IGBP palette

#IGBP=data.frame(ID=0:16,col=worldgrids_pal$IGBP[-c(18,19)],stringsAsFactors=F)

#IGBP$class=rownames(IGBP);rownames(IGBP)=1:nrow(IGBP)

 if(!file.exists("output/biomesamplepoints.csv")){

     biome=readOGR("data/teow/","biomes")

     n_biomesamples=1000

     library(multicore)

     biomesample=do.call(rbind.data.frame,mclapply(1:length(biome),function(i)

         data.frame(code=biome$code[i],coordinates(spsample(biome[i,],n=n_biomesamples,type="stratified",nsig=2)))))

     colnames(biomesample)[2:3]=c("lon","lat")

     biomesample[,c("biome","realm")]=biome@data[match(biomesample$code,biome$code),c("biome","realm")]

     write.csv(biomesample,"output/biomesamplepoints.csv",row.names=F)

 }

     ## Extract data for points

 if(!file.exists("output/biomesamplepoints_cloud.csv")){

     biomesample=read.csv("output/biomesamplepoints.csv")

     biomep=raster::extract(mod09c,biomesample,sp=T)

     biomep$lon=biomesample$lon

     biomep@data[,c("lon","lat")]=coordinates(biomep)

     write.csv(biomep@data,"output/biomesamplepoints_cloud.csv",row.names=F)

}     

biomep=read.csv("output/biomesamplepoints_cloud.csv")

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

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

at=seq(0,100,length=n)

## set plotting parameters

my.theme = trellis.par.get()

my.theme$strip.background=list(col="transparent")

trellis.par.set(my.theme)

png("output/CF_Figures_%03d.png",width=5000,height=4000,res=600,pointsize=42,bg="white")

levelplot(mod09a,col.regions=colr(n),cuts=100,at=seq(0,100,len=100),colorkey=list(title="Cloud Frequency (%)",space="bottom",adj=1),

levelplot(mod09a,col.regions=colr(n),cuts=100,at=seq(0,100,len=100),colorkey=list(title="Cloud Frequency (%)",space="bottom",adj=1),

levelplot(mod09s,col.regions=colr(n),cuts=100,at=seq(0,100,len=100),colorkey=list(title="Cloud Frequency (%)", space="bottom",adj=2),

levelplot(mod09metrics,col.regions=colr(n),cuts=100,at=seq(0,100,len=100),colorkey=list(title="Cloud Frequency (%)", space="bottom",adj=2),

levelplot(mod09c,col.regions=colr(n),cuts=100,at=seq(0,100,len=100),colorkey=list(title="Cloud Frequency (%)", space="bottom",adj=1), 

p_mac=levelplot(mod09a,col.regions=colr(n),cuts=99,at=seq(0,100,length=100),margin=F,maxpixels=res/10,colorkey=list(title="Cloud Frequency (%)", space="bottom",height=.75),xlab="",ylab="",useRaster=T)+

p_min=levelplot(mod09min,col.regions=colr(n),cuts=99,margin=F,maxpixels=res/10,colorkey=list(title="Cloud Frequency (%)", space="bottom",height=.75),useRaster=T)+

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

#for(r in 1:length(regs)){

#tmap=crop(wc_map,regs[[r]])

#p_map=levelplot(tmap,col.regions=grey(seq(0,1,len=100)),cuts=100,at=seq(tmap@data@min,tmap@data@max,len=100),margin=F,maxpixels=1e5,colorkey=list(space="bottom",height=.75),xlab="",ylab="",main=names(regs)[r],useRaster=T)

#tmac=crop(mod09a,regs[[r]])

#p_mac=levelplot(tmac,col.regions=grey(seq(0,1,len=100)),cuts=100,at=seq(tmac@data@min,tmac@data@max,len=100),margin=F,maxpixels=1e5,colorkey=list(space="bottom",height=.75),useRaster=T)

#p3=c("MOD09 Cloud Frequency (%)"=p_mac,"WorldClim Mean Annual Precip (mm)"=p_map,x.same=T,y.same=T,merge.legends=T,layout=c(2,1))

#print(p3)

#}

#dev.off()

names(gewex)="PATMOS-x GEWEX AVHRR"

## calculate 1-degree means of MODCF data

#MOD_gewex=gewex

#MOD_gewex@data@values=1:length(MOD_gewex@data@values)

#MOD_gewex2=zonal(mod09a,MOD_gewex,fun='mean')

png("output/Resolution_Figures_%03d.png",width=5500,height=2000,res=600,pointsize=36,bg="white")

trellis.par.set(my.theme)

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

r="Venezuela"

# ylab.right = "Cloud Frequency (%)",par.settings = list(layout.widths = list(axis.key.padding = 0.1,axis.left=0.6,ylab.right = 3,right.padding=2)),

p1=levelplot(crop(mod09a,regs[[r]]),col.regions=grey(seq(0,1,len=100)),at=seq(45,100,len=99),

    colorkey=list(space="bottom",width=1,height=1,labels=list(labels=c(50,75,100),at=c(50,75,100))),

    cuts=99,margin=F,max.pixels=1e5,par.settings = list(layout.heights=list(key.bottom=4),layout.widths = list(axis.key.padding = 0,axis.left=0.6)))

p2=levelplot(crop(gewex,regs[[r]]),col.regions=grey(seq(0,1,len=100)),at=seq(.45,1,len=99),cuts=99,margin=F,max.pixels=1e5,colorkey=F)

tmap=crop(wc_map,regs[[r]])

p3=levelplot(tmap,col.regions=grey(seq(0,1,len=100)),cuts=100,at=seq(tmap@data@min,tmap@data@max,len=100),margin=F,maxpixels=1e5,

    colorkey=list(space="bottom",height=.33,width=1,labels=list(labels=c(1000,4000),at=c(1000,4000))),xlab="",ylab="",main=names(regs)[r],useRaster=T)

print(c("MODCF (%)"=p1,"PATMOS-x GEWEX (%)"=p2,"WorldClim Precip (mm)"=p3,x.same=T,y.same=T,merge.legends=T,layout=c(3,1)))

dev.off()

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

### Some stats for paper

# approximate size of mod09ga archive - get total size for one day from the USGS website

size=scan("http://e4ftl01.cr.usgs.gov/MOLT/MOD09GA.005/2000.04.30/",what="char")

## extract all filesizes in MB (all the HDFs) and sum them and covert to TB for the length of the full record

sum(as.numeric(sub("M","",grep("[0-9]*M$",size,value=T))))/1024/1024*as.numeric(as.Date("2013-12-31")-as.Date("2000-02-24"))

## seasonal variability

cellStats(mod09sd,"mean")

### summary by biome

biomep$id=1:nrow(biomep)

biomepl=melt(biomep@data,id.vars=c("id","code","biome","realm"))

colnames(biomepl)[grep("variable",colnames(biomepl))]="month"

biomepl$month=factor(biomepl$month,ordered=T,levels=month.name)

biomepl$realm=factor(biomepl$realm,ordered=T,levels=c("Antarctic","Australasia","Oceania", "IndoMalay", "Neotropics","Palearctic","Nearctic" ))

biomepl$value[biomepl$value<0]=NA

png("output/Biome_Figures_%03d.png",width=5500,height=4000,res=600,pointsize=36,bg="white")

trellis.par.set(my.theme)

#[biomepl$id%in%sample(biomep$id,10000),]

p1=useOuterStrips(xyplot(value~month|realm+biome,groups=id,data=biomepl,panel=function(x,y,groups = groups, subscripts = subscripts){

    panel.xyplot(x,y,col=grey(0.6,0.2),type="l",lwd=.5,groups=groups,subscripts=subscripts)

    panel.smoother(y ~ s(x), method = "gam",lwd=2,subscripts=subscripts,n=24)

},scales=list(y=list(at=c(0,100),lim=c(-20,120),cex=.75,alternating=2,tck=c(0,1)),x=list(at=c(1,7),rot=45,alternating=1)),ylab="Biome",xlab.top="Geographic Realm",ylab.right="MODCF (%)", xlab="Month"),

    strip=strip.custom(par.strip.text = list(cex = .7)),strip.left=strip.custom(horizontal=TRUE,par.strip.text = list(cex = .75)))

p1$par.settings$layout.widths$strip.left[1]=13

p1$par.strip.text$lines=.65

print(p1)

dev.off()

library(maptools)

library(rgeos)

cld=cld[!is.na(cld$Amt),]

## table of stations with > 20 observations per month

cast(cld,StaID~YR,value="Nobs")

mtab=ddply(cld,c('StaID','month'),function(df){ data.frame(count=sum(df$Nobs>20,na.rm=T))})

mtab2=mtab[

    table(mtab$count>10)

stem(mtab$count)

      cldn_all=length(x$Amt[x$Nobs>=Nobsthresh]),

      cldsd=sd(x$Amt[x$YR>=2000&x$Nobs>=Nobsthresh],na.rm=T),

      cldn=length(x$Amt[x$YR>=2000&x$Nobs>=Nobsthresh]))}))

    cldm[,c("lat","lon")]=coordinates(st)[match(cldm$StaID,st$id),c("lat","lon")]

    biome=unionSpatialPolygons(teow,paste(teow$REALM,teow$BIOME,sep="_"), threshold=5)

    realmid=read.csv("/mnt/data/jetzlab/Data/environ/global/teow/official/realm.csv",stringsAsFactors=F,na.strings = "TTTT")

    dt=data.frame(code=row.names(biome),stringsAsFactors=F)

    dt[,c("realmid","biomeid")]=do.call(rbind,strsplit(sub(" ","",dt$code),"_"))

    dt$realm=realmid$realm[match(dt$realmid,realmid$realmid)]

    dt$biome=biomeid$BiomeID[match(dt$biomeid,biomeid$Biome)]

    row.names(dt)=row.names(biome)

    biome=SpatialPolygonsDataFrame(biome,data=dt)

st$biomec=biome$code[dists]

st$realm=biome$realm[dists]

st$biome=biome$biome[dists]

cldm$biomec=st$biomec[match(cldm$StaID,st$id)]

cldm$realm=st$relam[match(cldm$StaID,st$id)]

system(paste("cdo  -f nc4c -O inttime,2012-01-15,12:00:00,7day  -sellonlatbox,",