Project

General

Profile

« Previous | Next » 

Revision 41d291a6

Added by Adam Wilson almost 11 years ago

Added multiple tries to python GEE script

View differences:

climate/procedures/ee_compile.R
3 3 
setwd("~/acrobates/adamw/projects/cloud")
4 4 
library(raster)
5 5 
library(doMC)
6 

  		




  7
  
  
    
registerDoMC(4)


  		




  8
  
  
    
beginCluster(4)


  		




  
  6
  
    
library(multicore)


  		




  
  7
  
    
library(foreach)


  		




  
  8
  
    
#library(doMPI)


  		




  
  9
  
    
registerDoMC(10)


  		




  
  10
  
    
#beginCluster(4)


  		




  9
  11
  
    

  		




  10
  12
  
    
tempdir="tmp"


  		




  11
  13
  
    
if(!file.exists(tempdir)) dir.create(tempdir)


  		




  ......




  14
  16
  
    
## Load list of tiles


  		




  15
  17
  
    
tiles=read.table("tile_lat_long_10d.txt",header=T)


  		




  16
  18
  
    

  		




  17
  
  
    
jobs=expand.grid(tile=tiles$Tile,year=2000:2012,month=1:12)


  		




  18
  
  
    
jobs[,c("ULX","ULY","LRX","LRY")]=tiles[match(jobs$tile,tiles$Tile),c("ULX","ULY","LRX","LRY")]


  		




  
  19
  
    
### Build Tiles


  		




  
  20
  
    

  		




  
  21
  
    
## bin sizes


  		




  
  22
  
    
ybin=30


  		




  
  23
  
    
xbin=30


  		




  
  24
  
    

  		




  
  25
  
    
tiles=expand.grid(ulx=seq(-180,180-xbin,by=xbin),uly=seq(90,-90+ybin,by=-ybin))


  		




  
  26
  
    
tiles$h=factor(tiles$ulx,labels=paste("h",sprintf("%02d",1:length(unique(tiles$ulx))),sep=""))


  		




  
  27
  
    
tiles$v=factor(tiles$uly,labels=paste("v",sprintf("%02d",1:length(unique(tiles$uly))),sep=""))


  		




  
  28
  
    
tiles$tile=paste(tiles$h,tiles$v,sep="")


  		




  
  29
  
    
tiles$urx=tiles$ulx+xbin


  		




  
  30
  
    
tiles$ury=tiles$uly


  		




  
  31
  
    
tiles$lrx=tiles$ulx+xbin


  		




  
  32
  
    
tiles$lry=tiles$uly-ybin


  		




  
  33
  
    
tiles$llx=tiles$ulx


  		




  
  34
  
    
tiles$lly=tiles$uly-ybin


  		




  
  35
  
    
tiles$cy=(tiles$uly+tiles$lry)/2


  		




  
  36
  
    
tiles$cx=(tiles$ulx+tiles$urx)/2


  		




  
  37
  
    
tiles=tiles[,c("tile","h","v","ulx","uly","urx","ury","lrx","lry","llx","lly","cx","cy")]


  		




  
  38
  
    

  		




  
  39
  
    
jobs=expand.grid(tile=tiles$tile,year=2000:2012,month=1:12)


  		




  
  40
  
    
jobs[,c("ulx","uly","urx","ury","lrx","lry","llx","lly")]=tiles[match(jobs$tile,tiles$tile),c("ulx","uly","urx","ury","lrx","lry","llx","lly")]


  		




  
  41
  
    

  		




  
  42
  
    
## drop Janurary 2000 from list (pre-modis)


  		




  
  43
  
    
jobs=jobs[!(jobs$year==2000&jobs$month==1),]


  		




  
  44
  
    

  		




  
  45
  
    
#jobs=jobs[jobs$month==1,]


  		




  19
  46
  
    

  		




  20
  47
  
    

  		




  21
  
  
    
jobs=jobs[jobs$month==1,]


  		




  
  48
  
    
#jobs=jobs[jobs$month==7,]


  		




  22
  49
  
    
## Run the python downloading script


  		




  23
  50
  
    
#system("~/acrobates/adamw/projects/environmental-layers/climate/procedures/ee.MOD09.py -projwin -159 20 -154.5 18.5 -year 2001 -month 6 -region test")   


  		




  24
  
  
    
i=6715


  		




  25
  
  
    
testtiles=c("h02v07","h02v06","h02v08","h03v07","h03v06")


  		




  26
  
  
    
todo=which(jobs$tile%in%testtiles)


  		




  27
  
  
    
todo=todo[1:3]


  		




  
  51
  
    
i=1


  		




  
  52
  
    
#testtiles=c("h02v07","h02v06","h02v08","h03v07","h03v06","h03v08")


  		




  
  53
  
    
#todo=which(jobs$tile%in%testtiles)


  		




  
  54
  
    
#todo=todo[1:3]


  		




  
  55
  
    
#todo=1


  		




  28
  56
  
    
todo=1:nrow(jobs)


  		




  29
  
  
    
#todo=todo[500:503]


  		




  30
  
  
    
mclapply(todo,function(i)


  		




  31
  
  
    
         system(paste("~/acrobates/adamw/projects/environmental-layers/climate/procedures/ee.MOD09.py -projwin ",jobs$ULX[i]," ",jobs$ULY[i]," ",jobs$LRX[i]," ",jobs$LRY[i],


  		




  32
  
  
    
       "  -year ",jobs$year[i]," -month ",jobs$month[i]," -region ",jobs$tile[i],sep="")),mc.cores=1)


  		




  
  57
  
    

  		




  
  58
  
    
checkcomplete=T


  		




  
  59
  
    
if(checkcomplete&exists("df")){  #if desired (and "df" exists from below) drop complete date-tiles


  		




  
  60
  
    
todo=which(!paste(jobs$tile,jobs$year,jobs$month)%in%paste(df$region,df$year,df$month))


  		




  
  61
  
    
}


  		




  
  62
  
    

  		




  
  63
  
    
writeLines(paste("Tiling options will produce",nrow(tiles),"tiles and ",nrow(jobs),"tile-months.  Current todo list is ",length(todo)))


  		




  
  64
  
    

  		




  
  65
  
    
foreach(i=todo) %dopar%{


  		




  
  66
  
    
    system(paste("python ~/acrobates/adamw/projects/environmental-layers/climate/procedures/ee.MOD09.py -projwin ",


  		




  
  67
  
    
                      jobs$ulx[i]," ",jobs$uly[i]," ",jobs$urx[i]," ",jobs$ury[i]," ",jobs$lrx[i]," ",jobs$lry[i]," ",jobs$llx[i]," ",jobs$lly[i]," ",


  		




  
  68
  
    
                      "  -year ",jobs$year[i]," -month ",jobs$month[i]," -region ",jobs$tile[i],sep=""))


  		




  
  69
  
    
     }


  		




  33
  70
  
    

  		




  34
  71
  
    

  		




  35
  72
  
    
##  Get list of available files


  		




  36
  
  
    
df=data.frame(path=list.files("/mnt/data2/projects/cloud/mod09",pattern="*.tif$",full=T),stringsAsFactors=F)


  		




  
  73
  
    
df=data.frame(path=list.files("/mnt/data2/projects/cloud/mod09",pattern="*.tif$",full=T,recur=T),stringsAsFactors=F)


  		




  37
  74
  
    
df[,c("region","year","month")]=do.call(rbind,strsplit(basename(df$path),"_|[.]"))[,c(1,2,3)]


  		




  38
  75
  
    
df$date=as.Date(paste(df$year,"_",df$month,"_15",sep=""),"%Y_%m_%d")


  		




  39
  76
  
    

  		




  
  77
  
    
## add stats to test for missing data


  		




  
  78
  
    
addstats=F


  		




  
  79
  
    
if(addstats){


  		




  
  80
  
    
    df[,c("max","min","mean","sd")]=do.call(rbind.data.frame,mclapply(1:nrow(df),function(i) as.numeric(sub("^.*[=]","",grep("STATISTICS",system(paste("gdalinfo -stats",df$path[i]),inter=T),value=T)))))


  		




  
  81
  
    
    table(df$sd==0)


  		




  
  82
  
    
}


  		




  
  83
  
    

  		




  40
  84
  
    
## subset to testtiles?


  		




  41
  85
  
    
#df=df[df$region%in%testtiles,]


  		




  42
  
  
    
df=df[df$month==1,]


  		




  
  86
  
    
#df=df[df$month==1,]


  		




  43
  87
  
    
table(df$year,df$month)


  		




  44
  88
  
    

  		




  45
  89
  
    
## drop some if not complete


  		




  46
  
  
    
#df=df[df$year<=2009,]


  		




  
  90
  
    
df=df[df$month==1&df$year%in%c(2001:2002,2005),]


  		




  47
  91
  
    
rerun=T  # set to true to recalculate all dates even if file already exists


  		




  48
  92
  
    

  		




  49
  93
  
    
## Loop over existing months to build composite netcdf files


  		




  ......




  51
  95
  
    
## get date


  		




  52
  96
  
    
  print(date)


  		




  53
  97
  
    
  ## Define output and check if it already exists


  		




  
  98
  
    
  tffile=paste(tempdir,"/mod09_",date,".tif",sep="")


  		




  54
  99
  
    
  ncfile=paste(tempdir,"/mod09_",date,".nc",sep="")


  		




  55
  100
  
    
  if(!rerun&file.exists(ncfile)) next


  		




  56
  101
  
    
  ## merge regions to a new netcdf file


  		




  57
  
  
    
  system(paste("gdal_merge.py -o ",ncfile," -n -32768 -of netCDF -ot Int16 ",paste(df$path[df$date==date],collapse=" ")))


  		




  
  102
  
    
#  system(paste("gdal_merge.py -tap -init 5 -o ",ncfile," -n -32768.000 -of netCDF -ot Int16 ",paste(df$path[df$date==date],collapse=" ")))


  		




  
  103
  
    
  system(paste("gdal_merge.py -o ",tffile," -init -32768  -n -32768.000 -ot Int16 ",paste(df$path[df$date==date],collapse=" ")))


  		




  
  104
  
    
  system(paste("gdal_translate -of netCDF ",tffile," ",ncfile," -ot Int16 "))


  		




  
  105
  
    
  file.remove(tffile)


  		




  58
  106
  
    
  system(paste("ncecat -O -u time ",ncfile," ",ncfile,sep=""))


  		




  59
  107
  
    
## create temporary nc file with time information to append to MOD06 data


  		




  60
  108
  
    
  cat(paste("


  		




  ......




  102
  150
  
    

  		




  103
  151
  
    
### generate the monthly mean and sd


  		




  104
  152
  
    
#system(paste("cdo -P 10 -O merge -ymonmean data/mod09.nc -chname,CF,CF_sd -ymonstd data/mod09.nc data/mod09_clim.nc"))


  		




  105
  
  
    
system(paste("cdo  -O -ymonmean data/mod09.nc data/mod09_clim_mean.nc"))


  		




  
  153
  
    
xosystem(paste("cdo  -O -ymonmean data/mod09.nc data/mod09_clim_mean.nc"))


  		




  106
  154
  
    
system(paste("cdo  -O -chname,CF,CF_sd -ymonstd data/mod09.nc data/mod09_clim_sd.nc"))


  		




  107
  155
  
    

  		




  108
  156
  
    
#  Overall mean


  		












Also available in:  Unified diff