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

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

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

names(mod35c5)="C5MOD35CF"

NAvalue(mod35c5)=0

## mod35C6 annual

  system("gdalbuildvrt data/MOD35C6.vrt `find /home/adamw/acrobates/adamw/projects/interp/data/modis/mod35/summary/ -name '*h[1]*_mean.nc'` ")

  system("/usr/local/gdal-1.10.0/bin/gdalbuildvrt -a_srs '+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs' -sd 4 -b 1 data/MOD35C6_CFday_pmiss.vrt `find /home/adamw/acrobates/adamw/projects/interp/data/modis/mod35/summary/ -name '*h[1]*.nc'` ")

  system("gdalwarp data/MOD35C6_CFday_pmiss.vrt data/MOD35C6_CFday_pmiss.tif -r bilinear")

NAvalue(mod35c6)=255

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

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

#  lulc=ratify(lulc)

  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)

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

#mod09l=mask(mod09,land)

#mod35l=mask(mod35,land)

### compare MOD43 and MOD17 products

## MOD17

#extent(mod17)=alignExtent(mod17,mod09)

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"

## draw lulc-stratified random sample of mod35-mod09 differences 

#samp=sampleStratified(lulc, 1000, exp=10)

#save(samp,file="LULC_StratifiedSample_10000.Rdata")

#mean(dif[samp],na.rm=T)

#Stats(dif,function(x) c(mean=mean(x),sd=sd(x)))

###

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

#dif_c5_c6=mod35c5-mod35c6

dif_p=calc(trd1[[1]], function(x) (x[1]-x[3])/(1-x[1]))

edge=calc(edge(subset(t1,"MCD12Q1"),classes=T,type="inner"),function(x) ifelse(x==1,1,NA))

edgeb=buffer(edge,width=5000)

edgeb=calc(edgeb,function(x) ifelse(is.na(x),0,1))

names(edge)="edge"

names(edgeb)="edgeb"

pedge=calc(edge(subset(t1,"MOD35pp"),classes=T,type="inner"),function(x) ifelse(x==1,1,NA))

pedgeb=buffer(pedge,width=3000)

pedgeb=calc(pedgeb,function(x) ifelse(is.na(x),0,1))

names(pedgeb)="pedgeb"

td1=as.data.frame(stack(t1,edgeb,pedgeb))

td1l=melt(td1,id.vars=c("pedgeb","edgeb","MOD35pp","MCD12Q1"),measure.vars=c("C5MOD35CF","C6MOD35CF","C5MOD09CF"))

td1l=td1l[td1l$pedgeb==1|td1l$edgeb==1,]

cast(MOD35pp~MCD12Q1~variable,fun.aggregate="mean",data=td1l)

## Moving window

tiles=expand.grid(xmin=seq(-180,170,by=10),ymin=seq(-60,80,by=10))

tiles$xmax=tiles$xmin+10;tiles$ymax=tiles$ymin+10

registerDoMC(10)

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

writeLines(c(""), "log.txt")

nrw=nrow(mod35c5)

nby=10

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

writeLines(paste("Processing ",length(nrwg)," groups"))

output=mclapply(nrwg,function(ti){

  ## Extract focal areas 

  ngb=5

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

  pp_ind=getValuesFocal(pp,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

    round(kruskal.test(tval1,tind1)$p.value*100)         # if it works, return p.value*100

  })),nrow=nby,ncol=ncol(mod35c5),byrow=T))     # turn it back into a raster

  ## udpate raster and write it

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

  projection(pp_bias)=projection(mod35c5)

  writeRaster(pp_bias,file=paste("data/tiles/pp_bias_",ti,".tif",sep=""),

              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

    round(kruskal.test(tval1,tind1)$p.value*100)         # if it works, return p.value*100

  })),nrow=nby,ncol=ncol(mod35c5),byrow=T))     # turn it back into a raster

  ## udpate raster and write it

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

  projection(lulc_bias)=projection(mod35c5)

  writeRaster(lulc_bias,file=paste("data/tiles/lulc_bias_",ti,".tif",sep=""),

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

  return(ti)

},mc.cores=10)

## original solution

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

#    if(length(unique(na.omit(pp_ind[i,])))<2) return(255)

#    if(sort(table(pp_ind[i,]),dec=T)[2]<5) return(254)

#    kruskal.test(vals[i,],pp_ind[i,])$p.value*100

#  })),nrow=nrow(t_mod35c5),ncol=ncol(t_mod35c5),byrow=T))

#  extent(pp_bias)=extent(t_mod35c5)

#  writeRaster(pp_bias,file=paste("data/tiles/pp_bias_",ti,".tif",sep=""),

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

  ## Run kruskal test for processing lulc bias

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

#    if(length(unique(na.omit(lulc_ind[i,])))<2) return(NA)

#    if(sort(table(lulc_ind[i,]),dec=T)[2]<5) return(255)

#    kruskal.test(vals[i,],lulc_ind[i,])$p.value*100

#  })),nrow=nrow(t_mod35c5),ncol=ncol(t_mod35c5),byrow=T))

#  extent(lulc_bias)=extent(t_mod35c5)

#  writeRaster(lulc_bias,dataType="INT1U",file=paste("data/tiles/lulc_bias_",ti,".tif",sep=""),

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

### check raster temporary files

showTmpFiles()

#removeTmpFiles(h=1)

## merge all the files back again

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

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

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

  system("gdalwarp -co 'COMPRESS=LZW' -co 'ZLEVEL=9' data/pp_bias.vrt data/pp_bias.tif -r near")

#plot(stack(foc,x1))

pat=c(-0.02,seq(0,0.1,len=50),seq(0.1,1,len=50))

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

levelplot(stack(pp_bias,lulc_bias),col.regions=c("cyan",grayr2(100)),at=pat,colorkey=list(at=pat,cuts=100),margin=F)

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)

### Correlations

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

                        lty=c(0,1,1,1,1,5,1,5),

                        col=c("transparent","red","blue","black","darkgreen","darkgreen","orange","orange")),

                       text=list(

                         c("MODIS Products","C5 MOD09 % Cloudy","C5 MOD35 % Cloudy","C6 MOD35 % Cloudy","MOD17 % Missing","MOD17 (scaled)","MOD11 % Missing","MOD11 (scaled)")),

                       text=list(c("C5 MOD35 Processing Path","Desert","Land")),

                       text=list(c("MCD12Q1 IGBP Land Cover",IGBP$class[sort(unique(transd$value[transd$variable=="MCD12Q1"]+1))])))))),

#print(p4)

CairoPDF("output/mod35compare.pdf",width=11,height=7)

#CairoPNG("output/mod35compare_%d.png",units="in", width=11,height=8.5,pointsize=4000,dpi=1200,antialias="subpixel")

### Global Comparison

print(g2,position=c(0,0,1,0.415),more=F)

#print(g3,position=c(0.31,0.06,.42,0.27),more=F)