library(rgeos)

## read it in

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)

## get % cloudy

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

names(mod09)="MOD09CF"

NAvalue(mod09)=0

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

mod35c5=crop(mod35c5,mod09)

names(mod35c5)="C5MOD35CF"

NAvalue(mod35c5)=0

system("gdalinfo /home/adamw/acrobates/adamw/projects/interp/data/modis/mod35/summary/.nc")

## mod35C6 annual

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 1 -b 1 data/MOD35C6.vrt /home/adamw/acrobates/adamw/projects/interp/data/modis/mod35/summary/*2009mean.nc ")

system("align.sh data/MOD35C6.vrt data/MOD35_2009.tif data/MOD35C6_2009.tif")

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

mod35c6=crop(mod35c6,mod09)

names(mod35c6)="C6MOD35CF"

NAvalue(mod35c6)=255

## Monthly

#system("find /home/adamw/acrobates/adamw/projects/interp/data/modis/mod35/summary/  -name '*[0-9].nc' > data/summaryfiles.txt")

#system("/usr/local/gdal-1.10.0/bin/gdalbuidvrt -allow_projection_difference -b 1 -sd 1 data/MOD35C6_monthly.vrt  ")

##system("/usr/local/gdal-1.10.0/bin/gdalwarp --optfile data/summaryfiles.txt data/MOD35C6_monthly.tif  ")#

#system("/usr/local/gdal-1.10.0/bin/gdal_translate -stats -co \"COMPRESS=LZW\" -of GTiff data/MOD35C6_monthly.vrt data/MOD35C6_monthly.tif ")

## nobs

#system("/usr/local/gdal-1.10.0/bin/gdalbuildvrt -allow_projection_difference -sd 4 data/MOD35C6_monthlypmiss.vrt `find /home/adamw/acrobates/adamw/projects/interp/data/modis/mod35/summary/  -name '*[0-9].nc'` ")

#system("/usr/local/gdal-1.10.0/bin/gdal_translate -stats -co \"COMPRESS=LZW\" -of GTiff data/MOD35C6_monthlypmiss.vrt data/MOD35C6_monthlypmiss.tif ")

                                        #mod35c6=raster("")

## mask cloud masks to land pixels

#mod09l=mask(mod09,land)

#mod35l=mask(mod35,land)

names(mod17qc)="MOD17CF"

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

NAvalue(mod11)=0

names(mod11qc)="MOD11CF"

#mod43=raster("data/mod43_2009.tif",format="GTiff")

#mod43qc=raster("data/mod43_2009.tif",format="GTiff")

#rclmat=matrix( c(0, 0, 0, 1, 1, 11, 2, 2, 16, 3, 3, 1), ncol=3, byrow=TRUE)

#reclassify(pp,rclmat,file="data/MOD35_ProcessPath_reclass.tif",options=c("COMPRESS=LZW","ZLEVEL=9","PREDICTOR=2"),datatype="INT1U",overwrite=T)

#pp=raster("data/MOD35_ProcessPath_reclass.tif")

NAvalue(pp)=255

names(pp)="MOD35pp"

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

#pp=crop(pp,mod09)

#modprod=stack(mod35c5,mod09,pp,lulc)#,mod43,mod43qc)

#names(modprod)=c("MOD17","MOD17qc")#,"MOD43","MOD43qc")

dif_c5_09=mod35c5-mod09

dif_c6_09=mod35c6-mod09

dif_c5_c6=mod35c5-mod35c6

#hist(dif,maxsamp=1000000)

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

#r1=Lines(list(

#  Line(matrix(c(

#                -61.183,1.165,

#                -60.881,0.825

#                -60.106, 3.809,

#                -60.847,3.240

#                ),ncol=2,byrow=T))),"Venezuela")

                -61.688,4.098,

                -59.251,3.430

                33.195,12.512,

                33.802,12.894

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

r6=Lines(list(

  Line(matrix(c(

                -63.353,-10.746,

                -63.376,-9.310

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

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

### write out shapefiles of transects

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 

transb=gBuffer(trans,0.4)

## make polygons of bounding boxes

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

library(splancs)

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

  for(j in which(do.call(c,lapply(tdd,function(i) names(i)))%in%c("MOD11","MOD17"))){

    trange=cellStats(tdd[[j]],range)

    tdd[[j]]=100*(tdd[[j]]-trange[1])/(trange[2]-trange[1])

    tdd[[j]]@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")

#rondonia=trd[trd$trans=="Brazil",]

#trd=trd[trd$trans!="Brazil",]

at=seq(0,100,leng=100)

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

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

cols=bgyr(100)

## global map

library(maptools)

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

g1=levelplot(stack(mod35c5,mod35c6,mod09),xlab=" ",scales=list(x=list(draw=F),y=list(alternating=1)),col.regions=cols,at=at)+layer(sp.polygons(bbs[1:4],lwd=2))+layer(sp.lines(coast,lwd=.5))

#g2=levelplot(dif,col.regions=bgrayr(100),at=seq(-70,70,len=100),margin=F,ylab=" ",colorkey=list("right"))+layer(sp.polygons(bbs[1:4],lwd=2))+layer(sp.lines(coast,lwd=.5))

#trellis.par.set(background=list(fill="white"),panel.background=list(fill="white"))

#g3=histogram(dif,bg="white",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

p1=useOuterStrips(levelplot(value~x*y|variable+trans,data=trd[!trd$variable%in%c("MCD12Q1","MOD35pp"),],asp=1,scales=list(draw=F,rot=0,relation="free"),

                                       at=at,col.regions=cols,maxpixels=7e6,

                                       ylab="Latitude",xlab="Longitude"),strip = strip.custom(par.strip.text=list(cex=.75)))+layer(sp.lines(trans,lwd=2))

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

            ylab="",xlab=" "),strip = strip.custom(par.strip.text=list(cex=.75)),strip.left=F)+layer(sp.lines(trans,lwd=2))

p3=useOuterStrips(

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

            ylab="",xlab=" "),strip = strip.custom(par.strip.text=list(cex=.75)),strip.left=F)+layer(sp.lines(trans,lwd=2))

## transects

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

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

         td=transd[subscripts,]

         imod09=td$variable=="MOD09CF"

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

         imod35=td$variable=="C5MOD35CF"

         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)

         imod17=td$variable=="MOD17"

         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)

         imod11=td$variable=="MOD11"

         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)

         path=td[td$variable=="MOD35pp",]

         panel.segments(path$dist,-5,c(path$dist[-1],max(path$dist,na.rm=T)),-5,col=c("blue","cyan","tan","darkgreen")[path$value+1],subscripts=1:nrow(path),lwd=15,type="l")

         land=td[td$variable=="MCD12Q1",]

         panel.segments(land$dist,-10,c(land$dist[-1],max(land$dist,na.rm=T)),-10,col=IGBP$col[land$value+1],subscripts=1:nrow(land),lwd=15,type="l")

        },subscripts=T,par.settings = list(grid.pars = list(lineend = "butt")),

         x=list(alternating=1,relation="free"),#, lim=c(0,70)),

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

           labels=c("IGBP","MOD35",seq(0,100,len=5)),

           lim=c(-15,100))),

       xlab="Distance Along Transect (km)", ylab="% Missing Data / % of Maximum Value",

       legend=list(

         #inside=list(fun=draw.key(list(

         #                        lines=list(type="b",pch=16,cex=.5,lty=c(1,1,1,5,1,5),col=c("red","blue","darkgreen","darkgreen","orange","orange")),

         #                        text=list(c("MOD09 % Cloudy","MOD35 % Cloudy","MOD17 % Missing","MOD17 (scaled)","MOD11 % Missing","MOD11 (scaled)"),

         #       lwd=1,col=c("red","blue","darkgreen","darkgreen","orange","orange")))),corner=c(1,.34)) ,

        right=list(fun=draw.key(list(columns=1,#title="MCD12Q1",

                   lines=list(type=c("b","b","b","b","l","b","l"),pch=16,cex=.5,lty=c(1,1,1,1,5,1,5),col=c("red","blue","black","darkgreen","darkgreen","orange","orange",rep(NA,12))),

                     rectangles=list(border=NA,col=c(rep(NA,8),"tan","darkgreen",NA,IGBP$col[sort(unique(transd$value[transd$variable=="MCD12Q1"]+1))]) ),

                         text=list(c("MOD09 % Cloudy","C5 MOD35 % Cloudy","C6 MOD35 % Cloudy","MOD17 % Missing","MOD17 (scaled)","MOD11 % Missing","MOD11 (scaled)",

                           "\n \n \n C5 MOD35 Processing Path","Desert","Land","\n \n \n MCD12Q1 IGBP Land Cover",IGBP$class[sort(unique(transd$value[transd$variable=="MCD12Q1"]+1))])))))),

         strip = strip.custom(par.strip.text=list(cex=.75)))

#

print(p4)

#p5=levelplot(value~x*y|variable,data=rondonia,asp=1,scales=list(draw=F,rot=0,relation="free"),colorkey=T)#,

#print(p5)

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

#print(g1,position=c(0,.33,1,1),more=T)

#print(g2,position=c(0,0,1,0.394),more=T)

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

### MOD35 Desert Processing path

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

          at=c(-1:3),col.regions=c("blue","cyan","tan","darkgreen"),margin=F,

          colorkey=list(space="bottom",title="MOD35 Processing Path",labels=list(labels=c("Water","Coast","Desert","Land"),at=0:4-.5)))+layer(sp.polygons(bbs,lwd=2))+layer(sp.lines(coast,lwd=.5))

print(p1,position=c(0,0,.62,1),more=T)

print(p2,position=c(0.6,0.21,0.78,0.79),more=T)

print(p3,position=c(0.76,0.21,1,0.79))

### profile plots

print(p4)

dev.off()

### summary stats for paper

  td=cast(transect+loc+dist~variable,value="value",data=transd)

  td2=melt.data.frame(td,id.vars=c("transect","dist","loc","MOD35pp","MCD12Q1"))

## function to prettyprint mean/sd's

msd= function(x) paste(round(mean(x,na.rm=T),1),"% ±",round(sd(x,na.rm=T),1),sep="")

cast(td2,transect+variable~MOD35pp,value="value",fun=msd)

cast(td2,transect+variable~MOD35pp+MCD12Q1,value="value",fun=msd)

cast(td2,transect+variable~.,value="value",fun=msd)

cast(td2,transect+variable~.,value="value",fun=msd)

cast(td2,variable~MOD35pp,value="value",fun=msd)

cast(td2,variable~.,value="value",fun=msd)

td[td$transect=="Venezuela",]

## scatterplot of MOD35 vs MOD09

trdl=cast(trans+x+y~variable,value="value",data=trd)

xyplot(MOD35C5qc~MOD09qc|trans+as.factor(MOD35pp),pch=16,cex=.2,data=trdl,auto.key=T)+layer(panel.abline(0,1))

### LANDCOVER

levelplot(lulc,col.regions=IGBP$col,scales=list(cex=2),colorkey=list(space="right",at=0:17,labels=list(at=seq(0.5,16.5,by=1),labels=levels(lulc)[[1]]$class,cex=2)),margin=F)

#   glat=as.numeric(do.call(c,strsplit(sub("GRINGPOINTLATITUDE=","",system(paste("gdalinfo ",file," | grep GRINGPOINTLATITUDE"),intern=T)),split=",")))

#   glon=as.numeric(do.call(c,strsplit(sub("GRINGPOINTLONGITUDE=","",system(paste("gdalinfo ",file," | grep GRINGPOINTLONGITUDE"),intern=T)),split=",")))

#   bb=cbind(c(glon,glon[1]),c(glat,glat[1]))

#   pp = SpatialPolygons(list(Polygons(list(Polygon(bb)),1)))

#   proj4string(pp)=projection(glb)

#   ppc=gIntersection(pp,gpp)

#   ppc=gBuffer(ppc,width=0.3)  #buffer a little to remove gaps between images

ppc=gpp

###  Merge them into a geotiff

    system(paste("gdal_merge.py -v -init 255 -n 255 -o MOD35_ProcessPath_gdalmerge2.tif -co \"ZLEVEL=9\" -co \"COMPRESS=LZW\" -co \"PREDICTOR=2\" `ls -d -1 gridded/*.tif --sort=size `",sep=""))

system(paste("pkmosaic -co COMPRESS=LZW -co PREDICTOR=2 ",paste("-i",list.files("gridded",full=T,pattern="tif$"),collapse=" ")," -o MOD35_path_pkmosaic_mode.tif  -m 6 -v -t 255 -t 0 &"))

setwd("~/acrobates/adamw/projects/interp/data")

library(reshape)

library(sp);library(rgdal);

library(reshape)

library(ncdf4)

library(geosphere)

library(rgeos)

library(multicore);library(RSQLite)

stInfo=function(file,inventory,vars=vars){

  ### read in station location data

    cs=c(11,-1,8,-1,9,-1,6,-4,30,-1,3,-1,3,-1,5) #column widths, negative means skip

      print("Importing station data")

      st=read.fwf(file,widths=cs,colClasses="character",comment.char = "",header=F,skip=0)

      colnames(st)=c("id","lat","lon","elev","name","gsnflag","hcnflag","wmoid")

      print("Updating factors")

      for(i in c(5,6,7,8)) st[,i]=as.factor(st[,i])

      for(i in 2:4) st[,i]=as.numeric(st[,i])

    ### inventory

      print("Getting station start-stop years from station inventory file")

      cs=c(11,-1,8,-1,9,-1,4,-1,4,-1,4) #column widths, negative means skip

      stinv=read.fwf(inventory,widths=cs,colClasses="character",comment.char = "",header=F,skip=0)

      colnames(stinv)=c("id","lat","lon","variable","yearstart","yearstop")

      stinv=stinv[stinv$variable%in%vars,]

    ### reshape yearstart

      stinvw1=cast(stinv,id+lat+lon~variable,value="yearstart")

      colnames(stinvw1)[!grepl("id|lat|lon",colnames(stinvw1))]=

            paste(colnames(stinvw1)[!grepl("id|lat|lon",colnames(stinvw1))],".start",sep="")

    ### reshape yearstop

      stinvw2=cast(stinv,id+lat+lon~variable,value="yearstop")

      colnames(stinvw2)[!grepl("id|lat|lon",colnames(stinvw2))]=

            paste(colnames(stinvw2)[!grepl("id|lat|lon",colnames(stinvw2))],".stop",sep="")

      stinvw=merge(stinvw1,stinvw2)

      for(i in colnames(stinvw[,!grepl("id|lat|lon",colnames(stinvw))])) stinvw[,i]=as.numeric(as.character( stinvw[,i]))

      colnames(stinvw)=sub("TMAX","tmax",colnames(stinvw))

      colnames(stinvw)=sub("TMIN","tmin",colnames(stinvw))

      colnames(stinvw)=sub("PRCP","ppt",colnames(stinvw))

    ### merge the data

      print("Merging start-stop dates with station locations")

      st2=merge(st,stinvw[,!grepl("lat|lon",colnames(stinvw))],all.x=T,by=c("id"))

    ## add station type field

    st2$type=ifelse((!is.na(st2$tmax.start)|!is.na(st2$tmin.start))&!is.na(st2$ppt.start),"T&P",

      ifelse((!is.na(st2$tmax.start)|!is.na(st2$tmin.start))&is.na(st2$ppt.start),"T",

             ifelse((is.na(st2$tmax.start)|is.na(st2$tmin.start))&!is.na(st2$ppt.start),"P","None")))

    ### Convert to spatial points and return

      print("Convert to spatial points")

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

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

      proj4string(st2)=CRS("+proj=longlat +datum=WGS84")

      return(st2)

  }

### Process the station data

st=stInfo(file="ghcn/ghcnd-stations.txt",inventory="ghcn/ghcnd-inventory.txt",

    vars=c("PRCP","TMAX","TMIN"))

writeOGR(st,".","stationlocations",driver="ESRI Shapefile",overwrite=T)