mod35c6=raster("data/MOD35C6_2009.tif")

system("align.sh data/C5MOD35_ProcessPath.tif data/MOD09_2009.tif data/MOD35pp.tif")

r4=Lines(list(

trans=SpatialLines(list(r1,r2,r3,r4),CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs "))

t1=trd1[[1]]

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)

bgray=colorRampPalette(c("purple","blue","deepskyblue4"))

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

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

strip=strip.custom(par.strip.text=list(cex=.7),bg="transparent")

g1=levelplot(stack(mod35c5,mod35c6,mod09),xlab=" ",scales=list(x=list(draw=F),y=list(alternating=1)),

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

                                       ylab="Latitude",xlab="Longitude"),strip.left=strip,strip = strip)+layer(sp.lines(trans,lwd=2))

            ylab="",xlab=" "),strip = strip,strip.left=F)+layer(sp.lines(trans,lwd=2))

            ylab="",xlab=" "),strip = strip,strip.left=F)+layer(sp.lines(trans,lwd=2))

         y=list(at=c(-20,-10,seq(0,100,len=5)),

  strip = strip,strip.left=F)

#print(p4)

levelplot(pp,asp=1,scales=list(draw=T,rot=0),maxpixels=1e6,cuts=3,

          at=(0:3)+.5,col.regions=c("blue","cyan","tan","darkgreen"),margin=F,

          colorkey=list(space="top",title="MOD35 Processing Path",labels=list(labels=c("Water","Coast","Desert","Land"),at=0:3),height=.25))+

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

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

     plot.legend = TRUE,raster_name="Collection 5 MOD35 2009 Cloud Frequency",

library(raster)#;beginCluster(10)

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