## Figures associated with MOD35 Cloud Mask Exploration

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

library(raster);beginCluster(10)

library(rasterVis)

library(rgdal)

library(plotKML)

library(Cairo)

library(reshape)

## add tags for distribution

## MOD35

tags=c("TIFFTAG_IMAGEDESCRIPTION='Collection 5 MOD35 Cloud Frequency for 2009 extracted from MOD09GA state_1km bits 0-1. The MOD35 bits encode four categories (with associated confidence that the pixel is actually clear): confidently clear (confidence > 0.99), probably clear (0.99 >= confidence > 0.95), probably cloudy (0.95 >= confidence > 0.66), and confidently cloudy (confidence <= 0.66).  Following the advice of the MODIS science team (Frey, 2010), we binned confidently clear and probably clear together as clear and the other two classes as cloudy.  The daily cloud mask time series were summarized to mean cloud frequency (CF) by calculating the proportion of cloudy days during 2009'",

  "TIFFTAG_DOCUMENTNAME='Collection 5 MOD35 Cloud Frequency'",

  "TIFFTAG_DATETIME='20090101'",

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

system(paste("/usr/local/src/gdal-1.10.0/swig/python/scripts/gdal_edit.py data/MOD35_2009.tif ",paste("-mo ",tags,sep="",collapse=" "),sep=""))

## MOD09

tags=c("TIFFTAG_IMAGEDESCRIPTION='Collection 5 MOD09 Cloud Frequency for 2009 extracted from MOD09GA \'PGE11\' internal cloud mask algorithm (embedded in MOD09GA \'state_1km\' bit 10. The daily cloud mask time series were summarized to mean cloud frequency (CF) by calculating the proportion of cloudy days during 2009'",

  "TIFFTAG_DOCUMENTNAME='Collection 5 MOD09 Cloud Frequency'",

  "TIFFTAG_DATETIME='20090101'",

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

system(paste("/usr/local/src/gdal-1.10.0/swig/python/scripts/gdal_edit.py data/MOD09_2009.tif ",paste("-mo ",tags,sep="",collapse=" "),sep=""))

mod09=raster("data/MOD09_2009.tif")

mod35c5=raster("data/MOD35_2009.tif")

names(mod35c5)="C5MOD35CF"

NAvalue(mod35c5)=0

NAvalue(mod35c6)=255

               "data/MCD12Q1_IGBP_2009_051_wgs84_1km.tif -overwrite ",sep=""))}

lulc=raster("data/MCD12Q1_IGBP_2009_051_wgs84_1km.tif")

  data(worldgrids_pal)  #load palette

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

    lulc_levels2=c("Water","Forest","Forest","Forest","Forest","Forest","Shrublands","Shrublands","Savannas","Savannas","Grasslands","Permanent wetlands","Croplands","Urban and built-up","Cropland/Natural vegetation mosaic","Snow and ice","Barren or sparsely vegetated"),stringsAsFactors=F)

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

  levels(lulc)=list(IGBP)

names(lulc)="MCD12Q1"

## make land mask

if(!file.exists("data/land.tif"))

  land=calc(lulc,function(x) ifelse(x==0,NA,1),file="data/land.tif",options=c("COMPRESS=LZW","ZLEVEL=9","PREDICTOR=2"),datatype="INT1U",overwrite=T)

land=raster("data/land.tif")

### compare MOD43 and MOD17 products

## MOD17

system("align.sh ~/acrobates/adamw/projects/interp/data/modis/MOD17/MOD17A3_Science_NPP_mean_00_12.tif data/MOD09_2009.tif data/MOD17.tif")

mod17=raster("data/MOD17.tif",format="GTiff")

NAvalue(mod17)=65535

  system("align.sh ~/acrobates/adamw/projects/interp/data/modis/MOD17/MOD17A3_Science_NPP_Qc_mean_00_12.tif data/MOD09_2009.tif data/MOD17qc.tif")

mod17qc=raster("data/MOD17qc.tif",format="GTiff")

## MOD11 via earth engine

  system("align.sh ~/acrobates/adamw/projects/interp/data/modis/mod11/2009/MOD11_LST_2009.tif data/MOD09_2009.tif data/MOD11_2009.tif")

mod11=raster("data/MOD11_2009.tif",format="GTiff")

  system("align.sh ~/acrobates/adamw/projects/interp/data/modis/mod11/2009/MOD11_Pmiss_2009.tif data/MOD09_2009.tif data/MOD11qc_2009.tif")

mod11qc=raster("data/MOD11qc_2009.tif",format="GTiff")

### Processing path

names(pp)="MOD35pp"

###

n=100

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

cols=grey(seq(0,1,len=n))

cols=rainbow(n)

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

cols=bgyr(n)

### Transects

r1=Lines(list(

  Line(matrix(c(

                -59.251,3.430

r2=Lines(list(

  Line(matrix(c(

                133.746,-31.834,

                134.226,-32.143

                ),ncol=2,byrow=T))),"Australia")

r3=Lines(list(

  Line(matrix(c(

                73.943,27.419,

                74.369,26.877

                ),ncol=2,byrow=T))),"India")

                33.802,12.894

                ),ncol=2,byrow=T))),"Sudan")

writeOGR(SpatialLinesDataFrame(trans,data=data.frame(ID=names(trans)),match.ID=F),"output",layer="transects",driver="ESRI Shapefile",overwrite=T)

## buffer transects to get regional values 

## make polygons of bounding boxes

bb0 <- lapply(slot(transb, "polygons"), bbox)

bb1 <- lapply(bb0, bboxx)

# turn these into matrices using a helper function in splancs

bb2 <- lapply(bb1, function(x) rbind(x, x[1,]))

# close the matrix rings by appending the first coordinate

rn <- row.names(transb)

# get the IDs

bb3 <- vector(mode="list", length=length(bb2))

# make somewhere to keep the output

for (i in seq(along=bb3)) bb3[[i]] <- Polygons(list(Polygon(bb2[[i]])),

                   ID=rn[i])

# loop over the closed matrix rings, adding the IDs

bbs <- SpatialPolygons(bb3, proj4string=CRS(proj4string(transb)))

trd1=lapply(1:length(transb),function(x) {

  td=crop(mod11,transb[x])

  tdd=lapply(list(mod35c5,mod35c6,mod09,mod17,mod17qc,mod11,mod11qc,lulc,pp),function(l) resample(crop(l,transb[x]),td,method="ngb"))

  ## normalize MOD11 and MOD17

    tscaled@history=list(range=trange)

  return(brick(tdd))

})

## bind all subregions into single dataframe for plotting

trd=do.call(rbind.data.frame,lapply(1:length(trd1),function(i){

  d=as.data.frame(as.matrix(trd1[[i]]))

  d[,c("x","y")]=coordinates(trd1[[i]])

  d$trans=names(trans)[i]

  d=melt(d,id.vars=c("trans","x","y"))

  return(d)

}))

transd=do.call(rbind.data.frame,lapply(1:length(trans),function(l) {

  td=as.data.frame(extract(trd1[[l]],trans[l],along=T,cellnumbers=F)[[1]])

  td$loc=extract(trd1[[l]],trans[l],along=T,cellnumbers=T)[[1]][,1]

  td[,c("x","y")]=xyFromCell(trd1[[l]],td$loc)

  td$dist=spDistsN1(as.matrix(td[,c("x","y")]), as.matrix(td[1,c("x","y")]),longlat=T)

  td$transect=names(trans[l])

  td2=melt(td,id.vars=c("loc","x","y","dist","transect"))

  td2=td2[order(td2$variable,td2$dist),]

  # get per variable ranges to normalize

  tr=cast(melt.list(tapply(td2$value,td2$variable,function(x) data.frame(min=min(x,na.rm=T),max=max(x,na.rm=T)))),L1~variable)

  td2$min=tr$min[match(td2$variable,tr$L1)]

  td2$max=tr$max[match(td2$variable,tr$L1)]

  print(paste("Finished ",names(trans[l])))

transd$type=ifelse(grepl("MOD35|MOD09|CF",transd$variable),"CF","Data")

## comparison of % cloudy days

  overlay(mod35c5,mod09,fun=function(x,y) {return(x-y)},file="data/dif_c5_09.tif",format="GTiff",options=c("COMPRESS=LZW","ZLEVEL=9"),overwrite=T)

dif_c5_09=raster("data/dif_c5_09.tif",format="GTiff")

## Identify problematic areas with hard edges in cloud frequency

## set up processing chunks

nby=10

nrwg=seq(1,nrw,by=nby)

  ## Extract focal areas 

  ngb=5

  vals=getValuesFocal(mod35c5,ngb=ngb,row=ti,nrows=nby)

  lulc_ind=getValuesFocal(lulc,ngb=ngb,row=ti,nrows=nby)

  ## processing path

  pp_bias=raster(matrix(do.call(rbind,lapply(1:nrow(vals),function(i) {

    tind1=pp_ind[i,]  #vector of indices

    tval1=vals[i,]    # vector of values

    tind2=tind1[!is.na(tind1)&!is.na(tval1)]  #which classes exist without NAs?

    if(length(unique(tind2))<2) return(255)  #if only one class, return 255

    if(sort(table(tind2),dec=T)[2]<5) return(254) # if too few observations of class 2, return 254

    m=mean(tval1,na.rm=T)

    dif=round(diff(range(tapply(tval1,tind1,mean)))/m*100)

    return(dif)         # if it works, return p.value*100

  projection(pp_bias)=projection(mod35c5)

              format="GTiff",dataType="INT1U",overwrite=T,NAflag=255) #,options=c("COMPRESS=LZW","ZLEVEL=9")

  ## landcover

  lulc_bias=raster(matrix(do.call(rbind,lapply(1:nrow(vals),function(i) {

    tind1=lulc_ind[i,]  #vector of indices

    tval1=vals[i,]    # vector of values

    tind2=tind1[!is.na(tind1)&!is.na(tval1)]  #which classes exist without NAs?

    if(length(unique(tind2))<2) return(255)  #if only one class, return 255

    if(sort(table(tind2),dec=T)[2]<5) return(254) # if too few observations of class 2, return 254

    m=mean(tval1,na.rm=T)

    dif=round(diff(range(tapply(tval1,tind1,mean)))/m*100)

    return(dif)         # if it works, return the normalized difference

  ## udpate raster and write it

  extent(lulc_bias)=extent(mod35c5[ti:(ti+nby-1),1:ncol(mod35c5),drop=F])

  projection(lulc_bias)=projection(mod35c5)

    ## MOD09

  mod09_lulc_bias=raster(matrix(do.call(rbind,lapply(1:nrow(vals_mod09),function(i) {

    tind1=lulc_ind[i,]  #vector of indices

    tval1=vals_mod09[i,]    # vector of values

    tind2=tind1[!is.na(tind1)&!is.na(tval1)]  #which classes exist without NAs?

    if(length(unique(tind2))<2) return(255)  #if only one class, return 255

    if(sort(table(tind2),dec=T)[2]<5) return(254) # if too few observations of class 2, return 254

    m=mean(tval1,na.rm=T)

    dif=round(diff(range(tapply(tval1,tind1,mean)))/m*100)

    return(dif)         # if it works, return normalized difference

  ## udpate raster and write it

  extent(mod09_lulc_bias)=extent(mod09[ti:(ti+nby-1),1:ncol(mod09),drop=F])

  projection(mod09_lulc_bias)=projection(mod09)

  NAvalue(mod09_lulc_bias) <- 255

  writeRaster(mod09_lulc_bias,file=paste("data/tiles/mod09_lulc_bias_",ti,".tif",sep=""),

              format="GTiff",dataType="INT1U",overwrite=T,NAflag=255)#,options=c("COMPRESS=LZW","ZLEVEL=9")

### check raster temporary files

showTmpFiles()

#removeTmpFiles(h=1)

## merge all the files back again

  system("gdalwarp -srcnodata 255 -dstnodata 255 -multi -r bilinear -co 'COMPRESS=LZW' -co 'ZLEVEL=9' data/lulc_bias.vrt data/lulc_bias.tif -r near")

#  system("align.sh data/lulc_bias.vrt data/MOD09_2009.tif data/lulc_bias.tif")

  system("gdalbuildvrt data/pp_bias.vrt `find data/tiles -name 'pp_bias*tif'` ")

  system("align.sh -srcnodata 255 -dstnodata 255 -multi -r bilinear data/pp_bias.vrt data/MOD09_2009.tif data/pp_bias_align.tif &")

  system("gdalwarp -srcnodata 255 -dstnodata 255 -multi -r bilinear -co 'COMPRESS=LZW' -co 'ZLEVEL=9' data/mod09_lulc_bias.vrt data/mod09_lulc_bias.tif -r near")

pp_bias=raster("data/pp_bias.tif")

names(pp_bias)="Processing Path"

names(lulc_bias)="Land Use Land Cover"

## read in WWF biome data to summarize by biome

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

  system("wget -O data/teow.zip http://assets.worldwildlife.org/publications/15/files/original/official_teow.zip?1349272619")

  system("unzip data/teow.zip -d data/teow/")

   biome=readOGR("data/teow/wwf_terr_ecos.shp","wwf_terr_ecos")

   biome=biome[biome$BIOME<50,]

   biome2=gUnaryUnion(biome,id=biome$BIOME)

  ## create biome.csv using names in html file   

  biomeid=read.csv("data/teow/biome.csv",stringsAsFactors=F)

  biome2=SpatialPolygonsDataFrame(biome2,data=biomeid)

  writeOGR(biome2,"data/teow","biomes",driver="ESRI Shapefile")

}

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

biome2=extract(pp_bias,biome[biome$Biome==12,],df=T,fun=function(x) data.frame(mean=mean(x,na.rm=T),sd=sd(x,na.rm=T),prop=(sum(!is.na(x))/length(x))))

grayr2=colorRampPalette(c("grey","green","red"))#

#grayr2=colorRampPalette(grey(c(.75,.5,.25)))

          colorkey=F,margin=F,maxpixels=1e6)+layer(sp.lines(coast,lwd=.5))

          colorkey=T,margin=F,maxpixels=1e4)+layer(sp.lines(coast,lwd=.5))

histogram(lulc_bias)

cor(td1$MOD17,td1$C6MOD35,use="complete",method="spearman")

cor(td1$MOD17[td1$edgeb==1],td1$C5MOD35[td1$edgeb==1],use="complete",method="spearman")

bwplot(value~MOD35pp|variable,data=td1l[td1l$pedgeb==1,],horizontal=F)

#trdw=cast(trd,trans+x+y~variable,value="value")

#cor(trdw$MOD17,trdw$C5MOD35,use="complete",method="spearman")

#by(trdw,trdw$trans,function(x) cor(as.data.frame(na.omit(x[,c("C5MOD35CF","C6MOD35CF","C5MOD09CF","MOD17","MOD11")])),use="complete",method="spearman"))

#trdw_cor=as.data.frame(na.omit(trdw[,c("C5MOD35CF","C6MOD35CF","C5MOD09CF","MOD17","MOD11")]))

#nrow(trdw_cor)

#round(cor(trdw_cor,method="spearman"),2)

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

grayr=colorRampPalette(c("grey","red","darkred"))

bgrayr=colorRampPalette(c("darkblue","blue","grey","red","purple"))

library(maptools)

coast=map2SpatialLines(map("world", interior=FALSE, plot=FALSE),proj4string=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))

  col.regions=cols,at=at,cuts=length(at),maxpixels=1e6,

  colorkey=list(at=at))+

#  layer(sp.polygons(bbs,lwd=5,col="black"))+

  layer(sp.lines(coast,lwd=.5))+

  layer(sp.points(coordinates(bbs),col="black",cex=2,pch=13,lwd=2))

### Plot of differences between MOD09 adn MOD35 masks

#system("gdalinfo -stats /home/adamw/acrobates/adamw/projects/MOD35C5/data/dif_c5_09.tif")

## get quantiles for color bar of differences

#qs=unique(quantile(as.vector(as.matrix(dif_c5_09)),seq(0,1,len=100),na.rm=T))

#c(bgray(sum(qs<0)),grayr(sum(qs>=0)+1))

qs=seq(-80,80,len=100)

g2=levelplot(dif_c5_09,col.regions=bgrayr(100),cuts=100,at=qs,margin=F,ylab=" ",colorkey=list("right",at=qs),maxpixels=1e6)+

  layer(sp.points(coordinates(bbs),col="black",cex=2,pch=13,lwd=2))+

  #layer(sp.polygons(bbs,lwd=2))+

  layer(sp.lines(coast,lwd=.5))

#g3=histogram(dif_c5_09,col="black",border=NA,scales=list(x=list(at=c(-50,0,50)),y=list(draw=F),cex=1))+layer(panel.abline(v=0,col="red",lwd=2))

### regional plots

p2=useOuterStrips(

  levelplot(value~x*y|variable+trans,data=trd[trd$variable%in%c("MCD12Q1"),],

            asp=1,scales=list(draw=F,rot=0,relation="free"),colorkey=F,

            at=c(-1,IGBP$ID),col.regions=IGBP$col,maxpixels=7e7,

            legend=list(

              right=list(fun=draw.key(list(columns=1,#title="MCD12Q1 \n IGBP Land \n Cover",

                           rectangles=list(col=IGBP$col,size=1),

                           text=list(as.character(IGBP$ID),at=IGBP$ID-.5))))),

  levelplot(value~x*y|variable+trans,data=trd[trd$variable%in%c("MOD35pp"),],

            asp=1,scales=list(draw=F,rot=0,relation="free"),colorkey=F,

            at=c(-1:4),col.regions=c("blue","cyan","tan","darkgreen"),maxpixels=7e7,

            legend=list(

              right=list(fun=draw.key(list(columns=1,#title="MOD35 \n Processing \n Path",

                           rectangles=list(col=c("blue","cyan","tan","darkgreen"),size=1),

                           text=list(c("Water","Coast","Desert","Land")))))),

## transects

p4=xyplot(value~dist|transect,groups=variable,type=c("smooth","p"),

       data=transd,panel=function(...,subscripts=subscripts) {

         td=transd[subscripts,]

         panel.xyplot(td$dist[imod35],td$value[imod35],type=c("p","smooth"),span=0.09,subscripts=1:sum(imod35),col="blue",pch=16,cex=.25)

         ## mod35C6

         imod35c6=td$variable=="C6MOD35CF"

         panel.xyplot(td$dist[imod35c6],td$value[imod35c6],type=c("p","smooth"),span=0.09,subscripts=1:sum(imod35c6),col="black",pch=16,cex=.25)

         panel.xyplot(td$dist[imod17],100*((td$value[imod17]-td$min[imod17][1])/(td$max[imod17][1]-td$min[imod17][1])),

                      type=c("smooth"),span=0.09,subscripts=1:sum(imod17),col="darkgreen",lty=5,pch=1,cex=.25)

         imod17qc=td$variable=="MOD17CF"

         panel.xyplot(td$dist[imod17qc],td$value[imod17qc],type=c("p","smooth"),span=0.09,subscripts=1:sum(imod17qc),col="darkgreen",pch=16,cex=.25)

         panel.xyplot(td$dist[imod11],100*((td$value[imod11]-td$min[imod11][1])/(td$max[imod11][1]-td$min[imod11][1])),

                      type=c("smooth"),span=0.09,subscripts=1:sum(imod17),col="orange",lty="dashed",pch=1,cex=.25)

         imod11qc=td$variable=="MOD11CF"

         qcspan=ifelse(td$transect[1]=="Australia",0.2,0.05)

         panel.xyplot(td$dist[imod11qc],td$value[imod11qc],type=c("p","smooth"),npoints=100,span=qcspan,subscripts=1:sum(imod11qc),col="orange",pch=16,cex=.25)

           lim=c(-25,100)),

         alternating=F),

       legend=list(