Project

General

Profile

« Previous | Next » 

Revision d39ab57e

Added by Adam Wilson over 11 years ago

Submitted MOD35-landcover bias paper with code from this commit. Also added short script to test swtif program.

View differences:

climate/procedures/MOD35_ExtractProcessPath.r
1 1 
############################
2 2 
####  Extract MOD35 C6 processing path
3 3 

  		




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


  		




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


  		




  5
  5
  
    
library(multicore)


  		




  6
  6
  
    
library(raster)


  		




  7
  7
  
    
library(spgrass6)


  		




  ......




  11
  11
  
    
url="ftp://ladsweb.nascom.nasa.gov/allData/51/MOD35_L2/2012/"


  		




  12
  12
  
    
dir.create("swath")


  		




  13
  13
  
    

  		




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


  		




  
  14
  
    
getdata=F


  		




  
  15
  
    
if(getdata)


  		




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


  		




  15
  17
  
    

  		




  16
  18
  
    

  		




  17
  19
  
    
### make global raster that aligns with MODLAND tiles


  		




  18
  20
  
    
## get MODLAND tile to serve as base


  		




  19
  21
  
    
#system("wget http://e4ftl01.cr.usgs.gov/MOLT/MOD13A3.005/2000.02.01/MOD13A3.A2000032.h00v08.005.2006271174446.hdf")


  		




  20
  22
  
    
#t=raster(paste("HDF4_EOS:EOS_GRID:\"",getwd(),"/MOD13A3.A2000032.h00v08.005.2006271174446.hdf\":MOD_Grid_monthly_1km_VI:1 km monthly NDVI",sep=""))


  		




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


  		




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


  		




  22
  24
  
    
projection(t)


  		




  23
  25
  
    

  		




  24
  26
  
    
## make global extent


  		




  ......




  40
  42
  
    
#### Grid and mosaic the swath data


  		




  41
  43
  
    

  		




  42
  44
  
    
stitch="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"


  		




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


  		




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


  		




  44
  46
  
    

  		




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


  		




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


  		




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


  		




  47
  49
  
    

  		




  48
  50
  
    
## vars to process


  		




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


  		




  ......




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


  		




  58
  60
  
    
   proj4string(gpp)=projection(glb)


  		




  59
  61
  
    

  		




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


  		




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


  		




  61
  63
  
    

  		




  62
  64
  
    
swtif<-function(file,var){


  		




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


  		




  ......




  89
  91
  
    
   ## now run the swath2grid tool


  		




  90
  92
  
    
   ## write the gridded file


  		




  91
  93
  
    
   print(paste("Starting",file))


  		




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


  		




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


  		




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


  		




  94
  96
  
    
}  


  		




  95
  97
  
    

  		




  ......




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


  		




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


  		




  100
  102
  
    
   outfile=paste(outdir,"/",bfile,".tif",sep="")  #final file


  		




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


  		




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


  		




  102
  104
  
    
   ppc=gpp


  		




  103
  105
  
    
#######


  		




  104
  106
  
    
## run swtif for each band


  		




  ......




  113
  115
  
    
   d=raster(paste(tempdir(),"/CM_",basename(file),sep=""))


  		




  114
  116
  
    
   sz=raster(paste(tempdir(),"/SZ_",basename(file),sep=""))


  		




  115
  117
  
    
   NAvalue(sz)=-9999


  		




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


  		




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


  		




  117
  119
  
    
   path=  overlay(d,sz,fun=getlc,filename=tempfile2_path,options=c("COMPRESS=LZW", "LEVEL=9","PREDICTOR=2"),datatype="INT1U",overwrite=T)


  		




  118
  120
  
    
### warp them to align all pixels


  		




  119
  121
  
    
   system(paste("gdalwarp -overwrite -srcnodata 255 -dstnodata 255 -tap -tr 0.008333333 0.008333333 -co COMPRESS=LZW -co ZLEVEL=9 -co PREDICTOR=2 -s_srs \"",projection(t),"\" ",tempfile2_path," ",outfile,sep=""))


  		




  ......




  141
  143
  
    
file.remove(gfiles[check==0])


  		




  142
  144
  
    

  		




  143
  145
  
    
## use new gdal


  		




  144
  
  
    
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.tif &",sep=""))


  		




  
  146
  
    
#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=""))


  		




  
  147
  
    
#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=""))


  		




  
  148
  
    
#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=""))


  		




  145
  149
  
    

  		




  146
  150
  
    

  		




  147
  151
  
    
###  Merge them into a geotiff


  		




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


  		




  
  152
  
    
#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=""))


  		




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


  		




  149
  154
  
    

  		




  150
  
  
    
#  origin(raster(gfiles[5]))


  		




  151
  
  
    
  


  		




  
  155
  
    
 


  		




  152
  156
  
    
  ## try with pktools


  		




  153
  157
  
    
  ## global


  		




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


  		




  
  158
  
    
#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 &"))


  		




  155
  159
  
    
#bb="-ulx -180 -uly 90 -lrx 180 -lry -90"


  		




  156
  160
  
    
#bb="-ulx -180 -uly 90 -lrx 170 -lry 80"


  		




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


  		




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


  		




  158
  162
  
    

  		




  159
  163
  
    

  		




  160
  164
  
    
#expand.grid(x=seq(-180,170,by=10),y=seq(-90,80))


  		




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


  		




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


  		




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


  		




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


  		




  163
  167
  
    

  		




  164
  168
  
    
                                        #  bounding box?  


  		




  165
  169
  
    

  		




  ......




  184
  188
  
    
table(imported)


  		




  185
  189
  
    

  		




  186
  190
  
    
## read in all tifs


  		




  187
  
  
    
  for(f in gfiles[!imported])  execGRASS("r.in.gdal",input=f,output=basename(f),flags="o")


  		




  
  191
  
    
  for(f in gfiles[!imported])  {


  		




  
  192
  
    
    print(f)


  		




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


  		




  
  194
  
    
  }


  		




  
  195
  
    

  		




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


  		




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


  		




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


  		




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


  		




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


  		




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


  		




  
  202
  
    

  		




  
  203
  
    
##  Get mode of modes


  		




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


  		




  188
  205
  
    

  		




  189
  
  
    
## calculate mode


  		




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


  		




  191
  206
  
    
## add colors


  		




  192
  207
  
    
execGRASS("r.colors",map="MOD35_path",rules="MOD35_path_grasscolors.txt")


  		




  193
  208
  
    
## write to disk


  		












Also available in:  Unified diff