names(mod09)="C5MOD09CF"

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

               " /mnt/data/jetzlab/Data/environ/global/MODIS/MCD12Q1/051/MCD12Q1_051_2009_wgs84.tif ",

names(mod17)="MOD17_unscaled"

names(mod11)="MOD11_unscaled"

#r4=Lines(list(

#  Line(matrix(c(

#                -5.164,42.270,

#                -4.948,42.162

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

#r6=Lines(list(

#  Line(matrix(c(

#                -63.353,-10.746,

#                -63.376,-9.310

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

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

    names(tscaled)=sub("_unscaled","",names(tdd[[j]]))

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

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

## exploring various ways to compare cloud products along landcover or processing path edges

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

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

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

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

### Correlations

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

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

#Across all pixels in the four regions analyzed in Figure 3 there is a much larger correlation between mean NPP and the C5 MOD35 CF (Spearman’s ρ = -0.61, n=58,756) than the C6 MOD35 CF (ρ = 0.00, n=58,756) or MOD09 (ρ = -0.07, n=58,756) products.  

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

## table of correlations

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

#nrow(trdw_cor)

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

## set up some graphing parameters

g2$strip=strip.custom(var.name="Difference (C5MOD35-C5MOD09)",style=1,strip.names=T,strip.levels=F)  #update strip text

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

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

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

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

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

         imod09=td$variable=="C5MOD09CF"

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

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

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

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

           lim=c(-25,100)),

         alternating=F),

                      lines=list(type=c("b","b","b","b","b","l","b","l"),pch=16,cex=.5,

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

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

                       text=list(c("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)

#### export KML

kml_open("output/modiscloud.kml")

readAll(mod35c5)

kml_layer.Raster(mod35c5,

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

    z.lim = c(0,100),colour_scale = get("colour_scale_numeric", envir = plotKML.opts),

    altitudeMode = "clampToGround", balloon = FALSE

system(paste("gdal_translate -of KMLSUPEROVERLAY ",mod35c5@file@name," output/mod35c5.kmz -co FORMAT=JPEG"))

kml_close("modiscloud.kml")

kml_compress("modiscloud.kml",files=c(paste(month.name,".png",sep=""),"obj_legend.png"),zip="/usr/bin/zip")

setwd("~/acrobates/adamw/projects/interp/data/modis/mod35/processpath")

getdata=F

if(getdata)

  system(paste("wget -S --recursive --no-parent --no-directories -N -P swath --accept \"hdf\" --accept \"002|003|004\" ",url))

t=raster(paste("../../MOD17/MOD17A3_Science_NPP_mean_00_12.tif",sep=""))

#stitch="/usr/local/heg/2.12/bin/swtif"

files=paste(getwd(),"/",list.files("swath",pattern="hdf$",full=T),sep="")

outdir="~/acrobates/adamw/projects/interp/data/modis/mod35/processpath/gridded/"

   system(paste("sudo ",stitch," -p ",tempdir(),"/",basename(file),"_MODparms.txt -d -log /dev/null ",sep=""),intern=F)#,ignore.stderr=F)

   if(file.exists(outfile)) return(paste(file," already finished"))

   getlc=function(x,y) {ifelse(y<0|y>6000,NA,((x%/%2^6) %% 2^2))}

#system(paste("nohup /usr/local/gdal-1.10.0/bin/gdalwarp -wm 900 -overwrite -co COMPRESS=LZW -co PREDICTOR=2 -multi ",outdir,"/*.tif MOD35_path_gdalwarp.tif &",sep=""))

#system(paste("nohup /usr/local/gdal-1.10.0/bin/gdalwarp -wm 900 -overwrite -co COMPRESS=LZW -co PREDICTOR=2 -multi -r mode ",outdir,"/*.tif MOD35_path_gdalwarp_mode.tif &",sep=""))

#system(paste(" /usr/local/gdal-1.10.0/bin/gdalwarp -wm 900 -overwrite -co COMPRESS=LZW -co PREDICTOR=2 -multi ",outdir,"/MOD35_L2.A2012001*.tif MOD35_path_gdalwarp.tif",sep=""))

#system(paste("/usr/local/src/gdal-1.10.0/swig/python/scripts/gdal_merge.py -v -init 255 -n 255 -a_nodata 255 -o MOD35_ProcessPath_gdalmerge2.tif -co \"ZLEVEL=9\" -co \"COMPRESS=LZW\" -co \"PREDICTOR=2\" `ls -d -1 ",outdir,"/*.tif --sort=size | head -n 20 ` ",sep=""))

system(paste("/usr/local/src/gdal-1.10.0/swig/python/scripts/gdal_merge.py -v -o MOD35_ProcessPath_gdalmerge2.tif -co \"ZLEVEL=9\" -co \"COMPRESS=LZW\" -co \"PREDICTOR=2\" `ls -d -1 ",outdir,"/*.tif --sort=size ` &",sep=""))

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

#bb="-ulx -72 -uly 11 -lrx -59 -lry -1"

#gf2=  grep("2012009[.]03",gfiles,value=T)

#system(paste("pkmosaic ",bb," -co COMPRESS=LZW -co PREDICTOR=2 ",paste("-i",gf2,collapse=" ")," -o h11v08_path_pkmosaic.tif -ot Byte -m 7 -v -t 255"))

  for(f in gfiles[!imported])  {

    print(f)

    execGRASS("r.in.gdal",input=f,output=basename(f),flags="o")

  }

## calculate mode  - can't have more than 1000 open files

execGRASS("r.series",input=paste(system("g.mlist type=rast pattern=MOD*",intern=T)[1:1000],sep="",collapse=","),output="path1",method="mode",range=c(1,5),flags=c("verbose","overwrite"))

execGRASS("r.series",input=paste(system("g.mlist type=rast pattern=MOD*",intern=T)[1001:2000],sep="",collapse=","),output="path2",method="mode",range=c(1,5),flags=c("verbose","overwrite"))

execGRASS("r.series",input=paste(system("g.mlist type=rast pattern=MOD*",intern=T)[2001:3000],sep="",collapse=","),output="path3",method="mode",range=c(1,5),flags=c("verbose","overwrite"))

execGRASS("r.series",input=paste(system("g.mlist type=rast pattern=MOD*",intern=T)[3001:4000],sep="",collapse=","),output="path4",method="mode",range=c(1,5),flags=c("verbose","overwrite"))

execGRASS("r.series",input=paste(system("g.mlist type=rast pattern=MOD*",intern=T)[4001:5000],sep="",collapse=","),output="path5",method="mode",range=c(1,5),flags=c("verbose","overwrite"))

##  Get mode of modes

execGRASS("r.series",input=paste(system("g.mlist type=rast pattern=path*",intern=T),sep="",collapse=","),output="MOD35_path",method="mode",range=c(1,5),flags=c("verbose","overwrite"))

## Short script to test swtif gridding of various versions

setwd(tempdir())

library(sp)

library(raster)

system("wget ftp://ladsweb.nascom.nasa.gov/allData/6/MOD35_L2/2009/029/MOD35_L2.A2009029.0500.006.2012245113542.hdf")

system("wget ftp://ladsweb.nascom.nasa.gov/allData/6/MOD35_L2/2009/029/MOD35_L2.A2009029.0320.006.2012245113606.hdf")

files=list.files(pattern="hdf")

swtifpath="sudo MRTDATADIR=\"/usr/local/heg/2.12/data\" PGSHOME=/usr/local/heg/2.12/TOOLKIT_MTD PWD=/home/adamw /usr/local/heg/2.12/bin/swtif"

swtifpath="sudo MRTDATADIR=\"/usr/local/heg/2.12b/data\" PGSHOME=/usr/local/heg/2.12b/TOOLKIT_MTD PWD=/home/adamw /usr/local/heg/2.12b/bin/swtif"

## global bounding box

   gbb=cbind(c(-180,-180,180,180,-180),c(-90,90,90,-90,-90))

   gpp = SpatialPolygons(list(Polygons(list(Polygon(gbb)),1)))

   proj4string(gpp)="+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs "

vars=as.data.frame(matrix(c(

#  "Cloud_Mask",              "CM",       "NN",    1,

#  "Quality_Assurance",       "QA",       "NN",    1,

#  "Solar_Zenith",            "SolZen",   "NN", 1,

  "Sensor_Zenith",           "SenZen",   "CUBIC", 1

  ),

  byrow=T,ncol=4,dimnames=list(1:1,c("variable","varid","method","band"))),stringsAsFactors=F)

## define function that grids swaths

swtif<-function(file,var){

  outfile=paste(tempdir(),"/",var$varid,"_",basename(file),sep="")  #gridded path

   ## First write the parameter file (careful, heg is very finicky!)

   hdr=paste("NUM_RUNS = 1")

grp=paste("

BEGIN

INPUT_FILENAME=",file,"

OBJECT_NAME=mod35

FIELD_NAME=",var$variable,"|

BAND_NUMBER = ",var$band,"

OUTPUT_PIXEL_SIZE_X=926.6

OUTPUT_PIXEL_SIZE_Y=926.6

# MODIS 1km Resolution

SPATIAL_SUBSET_UL_CORNER = ( ",bbox(gpp)[2,2]," ",bbox(gpp)[1,1]," )

SPATIAL_SUBSET_LR_CORNER = ( ",bbox(gpp)[2,1]," ",bbox(gpp)[1,2]," )

RESAMPLING_TYPE =",var$method,"

OUTPUT_PROJECTION_TYPE = SIN

OUTPUT_PROJECTION_PARAMETERS = ( 6371007.181 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 )

# projection parameters from http://landweb.nascom.nasa.gov/cgi-bin/QA_WWW/newPage.cgi?fileName=sn_gctp

ELLIPSOID_CODE = WGS84

OUTPUT_TYPE = HDFEOS

OUTPUT_FILENAME = ",outfile,"

END

",sep="")

  ## write it to a file

  cat(c(hdr,grp)    , file=paste(tempdir(),"/",basename(file),"_MODparms.txt",sep=""))

  ## now run the swath2grid tool

  ## write the gridded file

  system(paste(swtifpath," -p ",tempdir(),"/",basename(file),"_MODparms.txt -d  -tmpLatLondir ",tempdir(),sep=""),intern=F,ignore.stderr=F)

   print(paste("Finished processing variable",var$variable,"from ",basename(file),"to",outfile))

}

swtifpath="sudo MRTDATADIR=\"/usr/local/heg/2.12/data\" PGSHOME=/usr/local/heg/2.12/TOOLKIT_MTD PWD=/home/adamw /usr/local/heg/2.12/bin/swtif"

swtifpath="sudo MRTDATADIR=\"/usr/local/heg/2.12b/data\" PGSHOME=/usr/local/heg/2.12b/TOOLKIT_MTD PWD=/home/adamw /usr/local/heg/2.12b/bin/swtif"

## run it

  lapply(1:nrow(vars),function(i) swtif(files[1],vars[i,]))

  lapply(1:nrow(vars),function(i) swtif(files[2],vars[i,]))

### make a plot to compare versions

library(rasterVis)

d=stack(list.files(pattern="SenZen.*hdf$"))

png(file="swtifCompareVersions.png",width=2000,height=1000)

levelplot(d,at=seq(0,100,len=100),col.regions=rainbow(100))

dif=d[[1]]-d[[2]]

levelplot(dif,at=seq(min(dif,na.rm=T),100,len=100),col.regions=rainbow(100))

dev.off()