Project

General

Profile

« Previous | Next » 

Revision 7239cefc

Added by Adam Wilson almost 11 years ago

updated 3_compile and uploaded first g3 dataset to maps engine for sharing

View differences:

climate/research/cloud/MOD09/2_bias.R
114 114 
foreach( i=1:nrow(jobs)) %dopar% {
115 115 
    file=jobs$path[i]
116 116 
    toutfile=paste(datadir,"mcd09bias/", sub(".tif","",basename(file)),"_",jobs$tile[i],".tif",sep="")
117 
    if(file.exists(toutfile)) {writeLines(paste(toutfile,"Exists, moving on"));next}
117 
#    if(file.exists(toutfile)) {writeLines(paste(toutfile,"Exists, moving on"));next}
118 118 
    writeLines(paste("Starting: ",toutfile," tile:",jobs$tile[i]," ( ",i," out of ",nrow(jobs),")"))
119 119 
    ## set sensor-specific parameters
120 120 
    ## add extra region for correction depending on which sensor is being processed
......
154 154 

  		




  155
  155
  
    
foreach( i=1:nrow(df2)) %dopar% {


  		




  156
  156
  
    
    ifile=df2$path[i]


  		




  157
  
  
    
    outfile=paste(datadir,"/mcd09ctif/",df2$sensor[i],"_",df2$month[i],".tif",sep="")


  		




  158
  
  
    
    if(file.exists(outfile)) {print(paste(outfile," exists, moving on...")); return(NULL)}


  		




  
  157
  
    
    outfile=paste(datadir,"/mcd09ctif/",df2$sensor[i],"_",df2$month[i],"_uncompressed.tif",sep="")


  		




  
  158
  
    
    outfile2=paste(datadir,"/mcd09ctif/",df2$sensor[i],"_",df2$month[i],".tif",sep="")


  		




  
  159
  
    
    ##    if(file.exists(outfile)) {print(paste(outfile," exists, moving on...")); return(NULL)}


  		




  159
  160
  
    
    ## create VRT of first band of the full image 


  		




  160
  161
  
    
    fvrt=sub("[.]tif",".vrt",ifile)


  		




  161
  162
  
    
    system(paste("gdalbuildvrt -b 1 ",fvrt," ",ifile))


  		




  162
  163
  
    
    ## mosaic the tiles with the original data (keeping the new data when available)


  		




  163
  
  
    
    tfiles=paste(c(fvrt,list.files(paste(datadir,"/mcd09bias",sep=""),pattern=paste(sub("[.]tif","",basename(outfile)),"_[0-9]*[.]tif",sep=""),full=T)),collapse=" ")


  		




  164
  
  
    
    if(file.exists(outfile)) file.remove(outfile)


  		




  165
  
  
    
    system(paste("gdal_merge.py --config GDAL_CACHEMAX 10000 -init -32768 -n -32768 -co COMPRESS=LZW -co PREDICTOR=2 -co BIGTIFF=yes  -o ",outfile," ",tfiles,sep="")) 


  		




  
  164
  
    
    tfiles=paste(c(fvrt,list.files(paste(datadir,"/mcd09bias",sep=""),pattern=paste(sub("[.]tif","",basename(outfile2)),"_[0-9]*[.]tif",sep=""),full=T)),collapse=" ")


  		




  
  165
  
    
    if(file.exists(outfile)) file.remove(outfile2,outfile)


  		




  
  166
  
    
    system(paste("gdal_merge.py --config GDAL_CACHEMAX 10000 -init -32768 -n -32768 -co BIGTIFF=yes  -o ",outfile," ",tfiles,sep="")) 


  		




  
  167
  
    
    system(paste("gdal_translate -co COMPRESS=LZW -co ZLEVEL=9 -co PREDICTOR=2 ",outfile," ",outfile2,sep=""))


  		




  
  168
  
    
    file.remove(fvrt,outfile)


  		




  166
  169
  
    
    writeLines(paste("Finished ",outfile))


  		




  167
  170
  
    
}


  		




  168
  171
  
    
 


  		












  

    
      climate/research/cloud/MOD09/3_combine.R
    
  





  11
  11
  
    
        tmcd2=paste(datadir,"/mcd09ctif/MCD09_",sprintf("%02d", i),".tif",sep="")


  		




  12
  12
  
    
        ## check if output already exists


  		




  13
  13
  
    
#        if(file.exists(tmcd2)){print(paste(tmcd2,"Exists, moving on..."));return(NULL)}


  		




  
  14
  
    
        if(file.exists(tmcd2)){print(paste(tmcd2,"Exists, deleting it..."));file.remove(tmcd,tmcd2)}


  		




  14
  15
  
    
        ## Take average between images


  		




  15
  16
  
    
        ## switch NA values to 32768 to facilitate recasting to 8-bit below, otherwise they are confounded with 0 cloud values


  		




  16
  17
  
    
        ops=paste("-t_srs 'EPSG:4326' -multi -srcnodata -32768 -dstnodata 32767 -r bilinear -te -180 -90 180 90 -tr 0.008333333333333 -0.008333333333333",


  		




  ......




  34
  35
  
    
        f=list.files(paste(datadir,"/mcd09tif",sep=""),pattern=paste(".*[O|Y].*_",sprintf("%02d",i),"[.]tif$",sep=""),full=T)


  		




  35
  36
  
    
        ## Define output and check if it already exists


  		




  36
  37
  
    
        tmcd=paste(datadir,"/mcd09ctif/MCD09sd_",sprintf("%02d", i),"_uncompressed.tif",sep="")


  		




  
  38
  
    
        tmcd=paste(datadir,"/mcd09ctif/MCD09sd_",sprintf("%02d", i),"_uncompressed.tif",sep="")


  		




  37
  39
  
    
        tmcd2=paste(datadir,"/mcd09ctif/MCD09sd_",sprintf("%02d", i),".tif",sep="")


  		




  38
  40
  
    
        ## check if output already exists


  		




  39
  41
  
    
        if(file.exists(tmcd2)){print(paste(tmcd2,"Exists, moving on..."));return(NULL)}


  		




  ......




  96
  98
  
    
}


  		




  97
  99
  
    

  		




  98
  100
  
    

  		




  
  101
  
    
################


  		




  
  102
  
    
### calculate inter vs. intra annual variability


  		




  
  103
  
    
f3=list.files(paste(datadir,"/mcd09ctif",sep=""),pattern=paste(".*MCD09_[0-9].[.]tif$",sep=""),full=T)


  		




  
  104
  
    
f3sd=list.files(paste(datadir,"/mcd09ctif",sep=""),pattern=paste(".*MCD09sd_[0-9].[.]tif$",sep=""),full=T)


  		




  
  105
  
    

  		




  
  106
  
    
dmean=stack(as.list(f3))


  		




  
  107
  
    
dsd=stack(as.list(f3sd))


  		




  
  108
  
    

  		




  
  109
  
    
dinter=calc(dmean,sd,file=paste(datadir,"/mcd09ctif/inter.tif",sep=""),options=c("COMPRESS=LZW","ZLEVEL=9"))


  		




  
  110
  
    
dintra=calc(dsd,mean,file=paste(datadir,"/mcd09ctif/intra.tif",sep=""),options=c("COMPRESS=LZW","ZLEVEL=9"))


  		




  
  111
  
    

  		




  99
  112
  
    
tplot=F


  		




  100
  113
  
    
if(tplot){


  		




  101
  114
  
    
    gcol=colorRampPalette(c("blue","yellow","red"))


  		












Also available in:  Unified diff