### Drop valitation station-months with fewer than 20 years of data for full record or less than 10 years for MODIS-era record

cldm$cld_all[cldm$cldn_all<20]=NA

cldm$cldsd_all[cldm$cldn_all<20]=NA

cldm$cld[cldm$cldn<10]=NA

cldm$cldsd[cldm$cldn<10]=NA

mod09inter=brick("data/cloud_std_inter.nc",varname="CF_sd")

mod09intra=brick("data/cloud_std_intra.nc",varname="CFsd")

## get stratified sample of points from biomes for illustration

 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"

## get stratified sample of points from biomes for illustration

 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"

## set plotting parameters

my.theme = trellis.par.get()

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

trellis.par.set(my.theme)

res=1e6

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

    layer(panel.polygon(x=c(-180,-180,180,180),y=c(-90,90,90,-90),col="black"),under=T)+

    layer(sp.lines(coast,col="black"),under=F)

    layer(panel.polygon(x=c(-180,-180,180,180),y=c(-90,90,90,-90),col="black"),under=T)+

    layer(panel.xyplot(lon,lat,pch=16,cex=.3,col="black"),data=data.frame(coordinates(st)))+

    layer(panel.polygon(x=c(-180,-180,180,180),y=c(-90,90,90,-90),col="black"),under=T)+

    layer(sp.lines(coast,col="black"),under=F)

    layer(panel.polygon(x=c(-180,-180,180,180),y=c(-90,90,90,-90),col="black"),under=T)+

    layer(sp.lines(coast,col="black"),under=F)

p_mac=levelplot(mod09a,col.regions=colr(n),cuts=99,at=seq(0,100,length=100),margin=F,

    ylim=greg$ylim,maxpixels=res/10,colorkey=list(title="Cloud Frequency (%)", space="top",height=.75),xlab="",ylab="",useRaster=T)+

    layer(panel.polygon(x=c(-180,-180,180,180),y=c(-90,90,90,-90),col="black"),under=T)+

    layer(sp.lines(coast,col="black"),under=F)

p_max=levelplot(mod09max,col.regions=colr(n),cuts=99,at=seq(0,100,length=100),margin=F,maxpixels=res/10,

    ylim=greg$ylim,colorkey=list(space="bottom",height=.75),useRaster=T)+

    layer(panel.polygon(x=c(-180,-180,180,180),y=c(-90,90,90,-90),col="black"),under=T)+

    layer(sp.lines(coast,col="black"),under=F)

p_intra=levelplot(mod09intra,col.regions=colr(n),cuts=99,at=seq(0,40,length=100),margin=F,maxpixels=res/100,

    ylim=greg$ylim,colorkey=list(space="bottom",height=.75),useRaster=T)+

    layer(panel.polygon(x=c(-180,-180,180,180),y=c(-90,90,90,-90),col="black"),under=T)+

    layer(sp.lines(coast,col="black"),under=F)

p_inter=levelplot(mod09inter,col.regions=colr(n),cuts=99,at=seq(0,15,length=100),margin=F,maxpixels=res/100,

    ylim=greg$ylim,colorkey=list(title="Cloud Frequency (%)", space="top",height=.75),useRaster=T)+

    layer(panel.polygon(x=c(-180,-180,180,180),y=c(-90,90,90,-90),col="black"),under=T)+

    layer(sp.lines(coast,col="black"),under=F)

p3=c("Mean Cloud Frequency (%)"=p_mac,"Max Cloud Frequency (%)"=p_max,"Interannual Variability (sd)"=p_inter,"Intraannual Variability (sd)"=p_intra,x.same=T,y.same=F,merge.legends=T,layout=c(2,2))

cldm$resid=NA

# get residuals of simple linear model

cldm$resid[!is.na(cldm$cld_all)&!is.na(cldm$mod09)]=residuals(lm(mod09~cld_all,data=cldm))

colat=bgr(cldm$resid)

phist=histogram(cldm$resid,breaks=colat$at,border=NA,col=colat$col,xlim=c(-30,30),type="count",xlab="MODCF Residuals")#,seq(0,1,len=6),na.rm=T)

pmap=xyplot(lat~lon|month2,data=cldm,groups=cut(cldm$resid,rev(colat$at)),

tr=c(0,66)

  panel.text(70,10,bquote(paste(R^2,"=",.(round(summary(lm(y ~ x))$r.squared,2)))),cex=1)

          ylab="NDP Mean Cloud Amount (%)",xlab="MODCF Cloud Frequency (%)",

wc_dem=stack(as.list(paste("/mnt/data/jetzlab/Data/environ/global/worldclim/alt.bil",sep="")))

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

pars=list(layout.heights=list(key.bottom=2,key.top=1),layout.widths = list(axis.key.padding = 3,axis.left=0.6))

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

    cuts=99,margin=F,max.pixels=1e6,par.settings = pars)

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=1e6,

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

    par.settings = pars)

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=1e6,

    colorkey=list(space="bottom",height=.75,width=1),xlab="",ylab="",main=names(regs)[r],useRaster=T,

    par.settings = pars)

p4=levelplot(crop(wc_dem,regs[[r]]),col.regions=grey(seq(0,1,len=100)),cuts=99,margin=F,max.pixels=1e6,

    colorkey=list(space="bottom",height=.75,width=1),

    par.settings = pars)#,labels=list(labels=c(1000,4000),at=c(1000,4000))))

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

## number of stations retained

length(unique(cldm$StaID[!is.na(cldm$cld_all)]))

length(unique(cldm$StaID[!is.na(cldm$cld)]))

### Compare time periods

     n <- as.integer(nobs(model))

lm_all1=lm(cld_all~mod09,data=cldm[!is.na(cldm$cld)&cldm$cldn_all>=10,])

lm_mod=lm(cld~mod09,data=cldm[cldm$cldn==10,])

mods=list("1970-2009"=lm_all1,"2000-2009"=lm_mod)

screenreg(mods,digits=2,single.row=T,custom.model.names=names(mods),custom.coef.names = c("Intercept", "MODCF"))

## assess latitude bias

cldm$abslat=abs(cldm$lat)

cldm$absdif=abs(cldm$difm)

abslm=lm(absdif~abslat*I(abslat^2),data=cldm[cldm$cldn_all>30,])

xyplot(absdif~abslat|month2,type=c("p","smooth"),data=cldm,cex=.25,pch=16)

plot(absdif~abslat,data=cldm[cldm$cldn_all>30,],cex=.25,pch=16)

##

## Summarize clouds by biome

setwd("~/acrobates/adamw/projects/cloud/")

library(rgdal)

library(raster)

library(geosphere)

library(rgeos)

library(maptools)

library(multicore)

## create Biome shapefile

if(!file.exists("data/teow/biomes.shp")){

    teow=readOGR("/mnt/data/jetzlab/Data/environ/global/teow/official/","wwf_terr_ecos")

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

    biomeid=read.csv("/mnt/data/jetzlab/Data/environ/global/teow/biome.csv",stringsAsFactors=F)

    realmid=read.csv("/mnt/data/jetzlab/Data/environ/global/teow/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)

    ## add area and centroid to each polygon

    biome$areakm=do.call(c,mclapply(1:length(biome),function(i) {print(i); return(areaPolygon(biome[i,])/1000000)}))

    biome@data[,c("lon","lat")]=coordinates(gCentroid(biome,byid=T))

    writeOGR(biome,"data/teow","biomes",driver="ESRI Shapefile",overwrite=T)

}

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

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

### load cloud data

#mod09s=brick("data/cloud_yseasmean.nc",varname="CF");names(mod09s)=c("DJF","MAM","JJA","SON")

#mod09c=brick("data/cloud_ymonmean.nc",varname="CF");names(mod09c)=month.name

#mod09s=brick("data/cloud_ymonstd.nc",varname="CFsd");names(mod09s)=month.name

mod09a=brick("data/cloud_mean.nc",varname="CF");names(mod09a)="Mean Annual Cloud Frequency"

mod09min=brick("data/cloud_min.nc",varname="CF");names(mod09min)="Min Cloud Frequency"

mod09max=brick("data/cloud_max.nc",varname="CF");names(mod09max)="Max Cloud Frequency"

mod09inter=brick("data/cloud_std_inter.nc",varname="CF_sd");names(mod09inter)="Inter-annual SD"

mod09intra=brick("data/cloud_std_intra.nc",varname="CFsd");names(mod09intra)="Intra-annual SD"

## Extract biome-level summaries

fsum=function(x,na.rm=T) paste(mean(x,na.rm=T),sd(x,na.rm=T),min(x,na.rm=T),max(x,na.rm=T),sep="_")

fs=list(min=mod09min,max=mod09max,intersd=mod09inter,intrasd=mod09intra,mean=mod09a)

fs=brick(list(min=mod09min,max=mod09max,intersd=mod09inter,intrasd=mod09intra,mean=mod09a))

biomed=do.call(rbind.data.frame,

    mclapply(1:length(fs),function(i) {

        print(i)

        td=extract(fs[[i]],biome,fun=fsum,df=F,sp=F)

        td2=do.call(rbind.data.frame,strsplit(td,"_"))

        td3=as.data.frame(apply(td2,2,as.numeric))

        colnames(td3)=c("mean","sd","min","max")

        td3$metric=names(fs)[i]

        td3$biome=biome$code

        return(td3)

    }))

write.csv(biomed,file="output/biomesummary.csv",row.names=F)

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

## get stratified sample of points from biomes for illustration

 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"                     

setwd("~/acrobates/adamw/projects/cloud")

df=data.frame(path=list.files("data/mod09cloud",pattern="*.tif$",full=T,recur=T),stringsAsFactors=F)

df[,c("month")]=do.call(rbind,strsplit(basename(df$path),"_|[.]|-"))[,c(6)]

df$date=as.Date(paste(2013,"_",df$month,"_15",sep=""),"%Y_%m_%d")

## use ramdisk?

tmpfs=tempdir()

ramdisk=T

if(ramdisk) {

    system("sudo mkdir -p /mnt/ram")

    system("sudo mount -t ramfs -o size=30g ramfs /mnt/ram")

    system("sudo chmod a+w /mnt/ram")

    tmpfs="/mnt/ram"

}

#i=2

for(i in 1:nrow(df)){

## Define output and check if it already exists

    temptffile=paste(tmpfs,"/modcf_",df$month[i],".tif",sep="")

    tffile=paste("data/mod09cloud2/modcf_",df$month[i],".tif",sep="")

    ncfile=paste("data/mod09cloud2/modcf_",df$month[i],".nc",sep="")

                                        #    if(file.exists(tffile)) next

    ## warp to wgs84

    ops=paste("-t_srs 'EPSG:4326' -multi -r bilinear -te -180 -90 180 90 -tr 0.008333333333333 -0.008333333333333",

        "-co BIGTIFF=YES -ot Byte --config GDAL_CACHEMAX 300000 -wm 300000 -wo NUM_THREADS:10 -wo SOURCE_EXTRA=5")

    system(paste("gdalwarp -overwrite ",ops," ",df$path[i]," ",temptffile))

    ## update metadata

    tags=c(paste("TIFFTAG_IMAGEDESCRIPTION='Cloud Frequency extracted from Collection 5 MYD09GA and MOD09GA (PGE11) internal cloud mask algorithm (embedded in M*D09GA state_1km bit 10).",

        " Band 1 represents the overall 2000-2013 mean cloud frequency for month ",df$month[i],

        ".  Band 2 is the four times the standard deviation of the mean monthly cloud frequencies over 2000-2013.'",sep=""),

        "TIFFTAG_DOCUMENTNAME='MODCF: MODIS Cloud Frequency'",

        paste("TIFFTAG_DATETIME='2013-",df$month[i],"-15'",sep=""),

        "TIFFTAG_ARTIST='Adam M. Wilson (adam.wilson@yale.edu)'")

    ## compress

#    trans_ops=paste(" -co COMPRESS=LZW -stats -co BLOCKXSIZE=256 -co BLOCKYSIZE=256 -co PREDICTOR=2 -co FORMAT=NC4C")

#    system(paste("gdal_translate ",trans_ops," ",paste("-mo ",tags,sep="",collapse=" ")," ",temptffile," ",tffile," ",sep=""))

    trans_ops=paste(" -co COMPRESS=DEFLATE -stats -co PREDICTOR=2 -co FORMAT=NC4C -co ZLEVEL=9")

    system(paste("gdal_translate ",trans_ops," ",temptffile," ",ncfile," ",sep=""))

file.remove(temptffile)

if(ramdisk) {

    ## unmount the ram disk

    system(paste("sudo umount ",tmpfs)

}

    ## Loop over existing months to build composite netcdf files

    foreach(date=unique(df$date)) %dopar% {

        ## get date

        print(date)

        ## Define output and check if it already exists

        temptffile=paste(tmpfs,"/modcf_",df$month[i],".tif",sep="")

        ncfile=paste("data/mod09cloud2/modcf_",df$month[i],".nc",sep="")

        ## check if output already exists

        if(!rerun&file.exists(ncfile)) return(NA)

        ## warp to wgs84

        ops=paste("-t_srs 'EPSG:4326' -multi -r bilinear -te -180 -90 180 90 -tr 0.008333333333333 -0.008333333333333",

            "-co BIGTIFF=YES -ot Byte --config GDAL_CACHEMAX 300000 -wm 300000 -wo NUM_THREADS:10 -wo SOURCE_EXTRA=5")

        system(paste("gdalwarp -overwrite ",ops," ",df$path[i]," ",temptffile))

        ## Warp to WGS84 grid and convert to netcdf

        trans_ops=paste(" -co COMPRESS=DEFLATE -stats -co FORMAT=NC4C -co ZLEVEL=9")

        system(paste("gdal_translate -of netCDF ",trans_ops," ",temptffile," ",ncfile))

        ## file.remove(temptffile)

        system(paste("ncecat -O -u time ",ncfile," ",ncfile,sep=""))

        ## create temporary nc file with time information to append to CF data

        cat(paste("

    }"),file=paste(tempdir(),"/",date,"_time.cdl",sep=""))

        system(paste("ncgen -o ",tempdir(),"/",date,"_time.nc ",tempdir(),"/",date,"_time.cdl",sep=""))

        ## add time dimension to ncfile and compress (deflate)

        system(paste("ncks --fl_fmt=netcdf4_classic -L 9 -A ",tempdir(),"/",date,"_time.nc ",ncfile,sep=""))

        ## add other attributes

        system(paste("ncrename -v Band1,CF ",ncfile,sep=""))

        system(paste("ncrename -v Band2,CFsd ",ncfile,sep=""))

        system(paste("ncatted ",

                     ## CF Mean

                     " -a units,CF,o,c,\"%\" ",

                     " -a valid_range,CF,o,b,\"0,100\" ",

                                        #               " -a scale_factor,CF,o,b,\"0.1\" ",

                     " -a _FillValue,CF,o,b,\"255\" ",

                     " -a missing_value,CF,o,b,\"255\" ",

                     " -a long_name,CF,o,c,\"Cloud Frequency (%)\" ",

                     " -a NETCDF_VARNAME,CF,o,c,\"Cloud Frequency (%)\" ",

                     ## CF Standard Deviation

                     " -a units,CFsd,o,c,\"SD\" ",

                     " -a valid_range,CFsd,o,b,\"0,200\" ",

                     " -a scale_factor,CFsd,o,f,\"0.25\" ",

                     " -a _FillValue,CFsd,o,b,\"255\" ",

                     " -a missing_value,CFsd,o,b,\"255\" ",

                     " -a long_name,CFsd,o,c,\"Cloud Frequency (%) Intra-month (2000-2013) Standard Deviation\" ",

                     " -a NETCDF_VARNAME,CFsd,o,c,\"Cloud Frequency (%) Intra-month (2000-2013) Standard Deviation\" ",

                     ## global

                     " -a title,global,o,c,\"Cloud Climatology from MOD09 Cloud Mask\" ",

                     " -a institution,global,o,c,\"Jetz Lab, EEB, Yale University, New Haven, CT\" ",

                     " -a source,global,o,c,\"Derived from MOD09GA Daily Data\" ",

                     " -a comment,global,o,c,\"Developed by Adam M. Wilson (adam.wilson@yale.edu / http://adamwilson.us)\" ",

                     ncfile,sep=""))

        ## add the fillvalue attribute back (without changing the actual values)

                                        #system(paste("ncatted -a _FillValue,CF,o,b,-32768 ",ncfile,sep=""))

        if(as.numeric(system(paste("cdo -s ntime ",ncfile),intern=T))<1) {

            print(paste(ncfile," has no time, deleting"))

            file.remove(ncfile)

        }

        print(paste(basename(ncfile)," Finished"))

    }

system(paste("cdo  -f nc4c -z zip -O -chname,CF,CF_sd -ymonstd -selmon,1,2,3,4,5,6 data/cloud_monthly.nc data/cloud_ymonsd_1-6.nc"))

system(paste("cdo  -f nc4c -z zip -O -chname,CF,CF_sd -ymonstd -selmon,7,8,9,10,11,12 data/cloud_monthly.nc data/cloud_ymonsd_7-12.nc"))

system(paste("cdo  -f nc4c -z zip -O mergetime  data/cloud_ymonsd_1-6.nc  data/cloud_ymonsd_7-12.nc data/cloud_ymonstd.nc"))

system("cdo -f nc4c -z zip  timmin data/cloud_ymonmean.nc data/cloud_min.nc")

system("cdo -f nc4c -z zip  timmax data/cloud_ymonmean.nc data/cloud_max.nc")

## standard deviation of mean monthly values give intra-annual variability

system("cdo -f nc4c -z zip -chname,CF,CFsd -timstd data/cloud_ymonmean.nc data/cloud_std_intra.nc")

## mean of monthly standard deviations give inter-annual variability 

system("cdo -f nc4c -z zip -chname,CF,CFsd -timmean data/cloud_ymonstd.nc data/cloud_std_inter.nc")