mod35c6=crop(mod35c6,mod09)

if(!file.exists("data/MCD12Q1_IGBP_2009_v5_wgs84_1km.tif")){

               " ~/acrobatesroot/jetzlab/Data/environ/global/landcover/MODIS/MCD12Q1_IGBP_2009_v5_wgs84.tif ",

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

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

  ## read it in

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

lulc=crop(lulc,mod09)

mod17=crop(mod17,mod09)

mod17qc=crop(mod17qc,mod09)

mod11qc=crop(mod11qc,mod09)

pp=crop(pp,mod09)

## comparison of % cloudy days

dif_c5_09=mod35c5-mod09

dif_c6_09=mod35c6-mod09

dif_c5_c6=mod35c5-mod35c6

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

    tdd[[j]]@history=list(range=trange)

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

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

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

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

  "Cloud_Mask",              "CM",

  "Sensor_Azimuth",       "ZA",

  "Sensor_Zenith",        "SZ"),

  byrow=T,ncol=2,dimnames=list(1:3,c("variable","varid"))),stringsAsFactors=F)

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

getpath<- function(file){  

   setwd(tempdir())

   bfile=sub(".hdf","",basename(file))

   tempfile_path=paste(tempdir(),"/path_",basename(file),sep="")  #gridded path

   tempfile_sz=paste(tempdir(),"/sz_",basename(file),sep="")  # gridded sensor zenith

   tempfile2_path=paste(tempdir(),"/",bfile,".tif",sep="")  #gridded/masked/processed path

   outfile=paste("~/acrobates/adamw/projects/interp/data/modis/mod35/gridded/",bfile,".tif",sep="")  #final file

   if(file.exists(outfile)) return(c(file,0))

   ## get bounding coordinates

#   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

   hdr=paste("NUM_RUNS = ",length(vars$varid),"|MULTI_BAND_HDFEOS:",length(vars$varid),sep="")

FIELD_NAME=",vars$variable,"|

BAND_NUMBER = ",1:length(vars$varid),"

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

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

OUTPUT_OBJECT_NAME = mod35|

RESAMPLING_TYPE =NN

OUTPUT_FILENAME= ",tempfile_path,"

   ## if any remnants from previous runs remain, delete them

   if(length(list.files(tempdir(),pattern=bfile)>0))

     file.remove(list.files(tempdir(),pattern=bfile,full=T))

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

##############  Now run the 5km summary

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

   hdr=paste("NUM_RUNS = ",nrow(vars[-1,]),"|MULTI_BAND_HDFEOS:",nrow(vars[-1,]),sep="")

   grp=paste("

BEGIN

INPUT_FILENAME=",file,"

OBJECT_NAME=mod35

FIELD_NAME=",vars$variable[-1],"|

BAND_NUMBER = ",1,"

OUTPUT_PIXEL_SIZE_X=0.008333333

#0.0416666

OUTPUT_PIXEL_SIZE_Y=0.008333333

#0.0416666

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

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

#OUTPUT_OBJECT_NAME = mod35|

RESAMPLING_TYPE =NN

OUTPUT_PROJECTION_TYPE = GEO

OUTPUT_PROJECTION_PARAMETERS = ( 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 0.0 )

#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= ",tempfile_sz,"

END

",sep="")

   ## write it to a file

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

   ## now run the swath2grid tool

   ## write the gridded file

   print(paste("Starting",file))

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

   if(!file.exists(tempfile_path)) {

   d=raster(paste("HDF4_EOS:EOS_GRID:\"",tempfile_path,"\":mod35:Cloud_Mask_0",sep=""))

   sz=raster(paste("HDF4_EOS:EOS_GRID:\"",tempfile_sz,"\":mod35:Sensor_Zenith",sep=""))

   NAvalue(sz)=0

   ## resample sensor angles to 1km grid and mask paths with angles >=30

#   sz2=resample(sz,d,method="ngb",file=sub("hdf","tif",tempfile_sz),overwrite=T)

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

## establish sudo priveleges to run swtif

system("sudo ls"); mclapply(files,getpath,mc.cores=10)

gfiles=list.files("gridded",pattern="tif$",full=T)

system(paste("/usr/local/gdal-1.10.0/bin/gdalwarp -wm 900 -overwrite -co COMPRESS=LZW -co PREDICTOR=2 -multi -r mode gridded/*.tif MOD35_path_gdalwarp.tif"))