Project

General

Profile

« Previous | Next » 

Revision b6e392da

Added by Adam Wilson almost 12 years ago

Identified problem in handling of missing values of MOD06 product. In some areas limiting quality control to the highest quality results in ~90% NAs.

View differences:

climate/procedures/Pleiades.R
1 1
#### Script to facilitate processing of MOD06 data
2 2

  
3 3
setwd("/nobackupp1/awilso10/mod06")
4

  
4 5
library(rgdal)
5 6
library(raster)
6 7

  
......
9 10
tb$tile=paste("h",sprintf("%02d",tb$ih),"v",sprintf("%02d",tb$iv),sep="")
10 11
save(tb,file="modlandTiles.Rdata")
11 12

  
13
## delete temporary log file that can grow to GB
14
system("rm /nobackupp1/awilso10/software/heg/TOOLKIT_MTD/runtime/LogStatus")
15

  
16
tile="h11v08"  # Venezuela
17
tile="h11v07"  # Venezuela coast
18
tile="h09v04"  # Oregon
19

  
12 20
outdir="2_daily"  #directory for separate daily files
13
outdir2="3_summary" #directory for combined daily files and summarized files
21
outdir=paste("daily/",tile,"/",sep="")  #directory for separate daily files
22
       if(!file.exists(outdir)) dir.create(outdir)
23
outdir2="summary" #directory for combined daily files and summarized files
24
       if(!file.exists(outdir2)) dir.create(outdir2)
14 25

  
15 26
  ## load a MOD11A1 file to define grid
16 27
gridfile=list.files("/nobackupp4/datapool/modis/MOD11A1.005/2006.01.27/",pattern=paste(tile,".*hdf$",sep=""),full=T)[1]
......
20 31

  
21 32
### get list of files to process
22 33
datadir="/nobackupp4/datapool/modis/MOD06_L2.005/"
34
datadir="/nobackupp1/awilso10/mod06/data"   #for data downloaded from 
23 35

  
24 36
fs=data.frame(path=list.files(datadir,full=T,recursive=T,pattern="hdf"),stringsAsFactors=F)
25 37
fs$file=basename(fs$path)
......
38 50

  
39 51
#### Generate submission file
40 52
## identify which have been completed
41
done=alldates%in%substr(list.files(outdir),7,14)
42
table(done)
43
notdone=alldates[!done]  #these are the dates that still need to be processed
53
#fdone=system(paste("ssh -q lou2 \"ls ",pdir," | grep \"nc$\"\"",sep=""),intern=T)
54
fdone=list.files(outdir,pattern="nc$")
55
done=alldates%in%substr(fdone,14,21)
44 56

  
45
if(exists("fdly")){
46
  notdone=notdone[notdone%in%fdly$dateid[is.na(fdly$npar)]]
57
if(exists("fdly")){  #update using table from below
58
  done[alldates%in%fdly$dateid[is.na(fdly$npar)]]=F
59
}
47 60

  
61
table(done)
62
notdone=alldates[!done]  #these are the dates that still need to be processed
48 63

  
49
tile="h11v08"   #can move this to submit script if needed
50 64
script="/u/awilso10/environmental-layers/climate/procedures/MOD06_L2_process.r"
51
write.table(paste("--verbose ",script," --date ",notdone," --tile \"",tile,"\"",sep=""),file="notdone.txt",row.names=F,col.names=F,quote=F)
52
#write.table(paste("--verbose ",script," date=",notdone[1:30],sep=""),file="notdone.txt",row.names=F,col.names=F,quote=F)
53
#write.table(notdone[1:30],file="notdone.txt",row.names=F,col.names=F,quote=F)
65

  
66
write.table(paste("--verbose ",script," --date ",notdone," --verbose T --tile \"",tile,"\"",sep=""),file=paste(tile,"_notdone.txt",sep=""),row.names=F,col.names=F,quote=F)
54 67

  
55 68
save(fs,alldates,gridfile,td,file="allfiles.Rdata")
56 69

  
57
## run script
58
#cat(paste("
59
##! /bin/bash
60
#source ~/moduleload
61
#source ~/.bashrc
62
#Rscript --verbose --vanilla /u/awilso10/environmental-layers/climate/procedures/MOD06_L2_process.r date=$1
63
#Rscript --verbose --vanilla rtest
64
#",sep=""),file="MOD06_process2")
65
#system("chmod +x MOD06_process2")
66

  
67
#cat(paste("
68
#library(rgdal)
69
#GDALinfo
70
#",sep=""),file="rtest")
71

  
72

  
73
## Submission script
74

  
75
#cat(paste("
76
##PBS -S /bin/bash
77
##PBS -l select=2:ncpus=16:model=san
78
####PBS -l select=4:ncpus=8:model=neh
79
###PBS -l select=1:ncpus=12:model=wes
80
######## old: select=48:ncpus=8:mpiprocs=8:model=neh
81
##PBS -l walltime=2:00:00
82
##PBS -j oe
83
##PBS -m e
84
##PBS -V
85
##PBS -q devel
86
##PBS -o log/log_^array_index^
87
##PBS -o log/log_DataCompile.log
88
##PBS -M adam.wilson@yale.edu
89
##PBS -N MOD06
90

  
91
## cd to working directory
92
#cd /nobackupp1/awilso10/mod06
93

  
94
## set some memory limits
95
#  ulimit -d 1500000 -m 1500000 -v 1500000  #limit memory usage
96
#  source /usr/local/lib/global.profile
97
#  source /u/awilso10/.bashrc
98
#  source /u/awilso10/moduleload
99
### export a few important variables
100
#  export NNODES=32
101
#  export R_LIBS=\"/u/awilso10/R/x86_64-unknown-linux-gnu-library/2.15/\"
102
### Run the script!
103
## current version not parallelizing across nodes!
104
#  TMPDIR=$TMPDIR Rscript --verbose --vanilla /u/awilso10/environmental-layers/climate/procedures/MOD06_L2_process.r date=20000403
105

  
106
#WORKLIST=notdone.txt
107
##EXE=\"Rscript\"
108
#EXE="./MOD06_process2"
109
#LOG=log/log_DataCompile.log
110
#MQUEUE=/nobackupp4/pvotava/software/share/mqueue-eg/mqueue/mqueue
111

  
112
#TMPDIR=$TMPDIR mpiexec -np $NNODES $MQUEUE -l $WORKLIST -p $EXE -v -v -v --random-starts 2-4 --work-analyze #> $LOG
113
#exit 0
114
#",sep=""),file="MOD06_process")
115

  
116

  
117
### simplified qsub script
70
### qsub script
118 71
cat(paste("
119 72
#PBS -S /bin/bash
120
#PBS -l select=48:ncpus=8:mpiprocs=8
73
#PBS -l select=40:ncpus=8:mpiprocs=8
121 74
##PBS -l select=2:ncpus=4:mpiprocs=4
122
#PBS -l walltime=02:00:00
75
#PBS -l walltime=2:00:00
123 76
#PBS -j n
124 77
#PBS -m be
125 78
#PBS -N mod06
126 79
#PBS -q devel
127 80
#PBS -V
128 81

  
129
CORES=384
82
CORES=320
130 83
HDIR=/u/armichae/pr/
131 84
  source $HDIR/etc/environ.sh
132 85
  source /u/awilso10/.bashrc
133 86
IDIR=/nobackupp1/awilso10/mod06/
134 87
##WORKLIST=$HDIR/var/run/pxrRgrs/work.txt
135
WORKLIST=$IDIR/notdone.txt
88
WORKLIST=$IDIR/",tile,"_notdone.txt
136 89
EXE=Rscript
137 90
LOGSTDOUT=$IDIR/log/log.stdout
138 91
LOGSTDERR=$IDIR/log/log.stderr
139 92
mpiexec -np $CORES pxargs -a $WORKLIST -p $EXE -v -v -v --work-analyze 1> $LOGSTDOUT 2> $LOGSTDERR
140
"),file="mod06_qsub")
93
",sep=""),file=paste(tile,"_mod06_qsub",sep=""))
141 94

  
142 95

  
143 96
### Check the file
144
system("cat mod06_qsub")
97
system(paste("cat ",tile,"_mod06_qsub",sep=""))
145 98
#system("cat ~/environmental-layers/climate/procedures/MOD06_L2_process.r")
146 99

  
147 100
## check queue status
......
149 102
system("/u/scicon/tools/bin/qtop.pl 492352")
150 103

  
151 104
## Submit it (and keep the pid)!
152
system("qsub mod06_qsub")
153
system("/u/scicon/tools/bin/pdsh_gdb -j 568835 -d tmp -s -u awilso10")
105
system(paste("qsub ",tile,"_mod06_qsub",sep=""))
154 106

  
155 107
## work in interactive mode
156 108
# system("qsub -I -l walltime=2:00:00 -lselect=2:ncpus=16:model=san -q devel")
......
161 113
system(paste("/u/scicon/tools/bin/qps ",568835))
162 114
system(paste("qstat -t -x",pid))
163 115

  
164
system("qstat devel ") 
165
#system("qstat | grep awilso10") 
166 116

  
167 117
####################################
168 118

  
119
  ## move to permanent storage on lou?
120
#  system(paste("scp -r ",ncfile," ",pdir,sep=""))
121
file.remove(
122
list.files(pattern="filetable[.]temp|GetAttrtemp|core[.]|MCFWrite[.]temp")
123
)
169 124

  
170 125
################################################################################
171 126
## now generate the climatologies
172 127
fdly=data.frame(
173 128
  path=list.files(outdir,pattern="nc$",full=T),
174 129
  file=list.files(outdir,pattern="nc$"))
175
  fdly$dateid=substr(fdly$file,7,14)
176
  fdly$date=as.Date(substr(fdly$file,7,14),"%Y%m%d")
130
  fdly$dateid=substr(fdly$file,14,21)
131
  fdly$date=as.Date(substr(fdly$file,14,21),"%Y%m%d")
177 132
  fdly$month=format(fdly$date,"%m")
178 133
  fdly$year=format(fdly$date,"%Y")
134
nrow(fdly)
179 135

  
180 136
## check validity (via npar and ntime) of nc files
181 137
for(i in 1:nrow(fdly)){
182
 # fdly$ntime[i]<-as.numeric(system(paste("cdo  sinfo ",fdly$path[i]),intern=T))
138
  fdly$ntime[i]<-as.numeric(system(paste("cdo -s ntime ",fdly$path[i]),intern=T))
183 139
  fdly$npar[i]<-as.numeric(system(paste("cdo -s npar ",fdly$path[i]),intern=T))
184
  print(i)
140
  fdly$fyear[i]<-as.numeric(system(paste("cdo -s showyear ",fdly$path[i]),intern=T))
141
  fdly$fmonth[i]<-as.numeric(system(paste("cdo -s showmon ",fdly$path[i]),intern=T))
142
  fdly$fvar[i]<-system(paste("cdo -s showvar ",fdly$path[i]),intern=T)
143
  print(paste(i," out of ",nrow(fdly)," for year ",  fdly$fyear[i]))
144
}
145

  
146
### table of problematic files
147
table(is.na(fdly$npar))
148
table(fdly$fmonth)
149
table(fdly$fvar)
150

  
151
fdly[is.na(fdly$npar),]
152

  
153
## delete files that fail check?
154
delete=T
155
if(delete) {
156
  drop=is.na(fdly$npar)|fdly$fvar!=" CER COT CLD"
157
  file.remove(as.character(fdly$path[drop]))
158
  fdly=fdly[!drop,]
185 159
}
186 160

  
187 161
## Combine all days within years into a single file (can't mergetime all days at once because this opens too many files)
162
ncopath="/nasa/sles11/nco/4.0.8/gcc/mpt/bin/"
163

  
164
library(multicore)
165

  
166
## create temporary directory to put intermediate files (will be deleted when R quits)
188 167
tsdir=paste(tempdir(),"/summary",sep="")
189 168
dir.create(tsdir)
190
lapply(unique(fdly$year),function(y){
191
  system(paste("cdo -O mergetime ",paste(fdly$path[!is.na(fdly$npar)&fdly$year==y],collapse=" ")," ",tsdir,"/MOD09_",tile,"_",y,"_daily.nc",sep=""))
192
  print(paste("Finished merging daily files for year",y))
193
})
194
## Combine the year-by-year files into a single daily file
195
system(paste("cdo -O mergetime ",paste(list.files(tsdir,full=T,pattern="daily[.]nc$"),collapse=" ")," ",outdir2,"/MOD09_",tile,"_daily.nc",sep=""))
196 169

  
197
system(paste("cdo -O monmean ",outdir2,"/MOD09_",tile,"_daily.nc ",outdir2,"/",tile,"_monmean.nc",sep=""))
198
system(paste("cdo -O ymonmean ",outdir2,"/MOD09_",tile,"_daily.nc ",outdir2,"/",tile,"_ymonmean.nc",sep=""))
199
system(paste("cdo -O ymonstd ",outdir2,"/MOD09_",tile,"_daily.nc ",outdir2,"/",tile,"_ymonstd.nc",sep=""))
170
## create a file of daily scenes for each year
171
#mclapply(unique(fdly$year),function(y){
172
#  system(paste("cdo -O mergetime ",paste(fdly$path[!is.na(fdly$npar)&fdly$year==y],collapse=" ")," ",tsdir,"/MOD06_",tile,"_",y,"_daily.nc",sep=""))
173
#  print(paste("Finished merging daily files for year",y))
174
#},mc.cores=5)
175

  
176

  
177
## merge all daily files to create a single file with all dates
178
  system(paste(ncopath,"ncrcat -O ",outdir,"/*nc ",outdir2,"/MOD06_",tile,"_daily.nc",sep=""))
179

  
180
##############################
181
## Combine the year-by-year files into a single daily file in the summary directory (for archiving)
182
#system(paste("cdo -O mergetime ",paste(list.files(tsdir,full=T,pattern="daily[.]nc$"),collapse=" ")," ",outdir2,"/MOD06_",tile,"_daily.nc",sep=""))
183
system(paste(ncopath,"ncatted ",
184
" -a units,time,o,c,\"days since 2000-1-1 0:0:0\" ",
185
" -a title,global,o,c,\"MODIS Cloud Product (MOD06) Daily Timeseries\" ",
186
" -a institution,global,o,c,\"Yale University\" ",
187
" -a source,global,o,c,\"MODIS Cloud Product (MOD06)\" ",
188
" -a comment,global,o,c,\"Compiled by Adam M. Wilson (adam.wilson@yale.edu)\" ",
189
outdir2,"/MOD06_",tile,"_daily.nc",sep=""))
190

  
191
#system(paste("cdo -O monmean ",outdir2,"/MOD06_",tile,"_daily.nc ",tsdir,"/MOD06_",tile,"_monmean.nc",sep=""))
192

  
193
#############################
194
##  Generate the Climatologies
195
myear=as.integer(max(fdly$year))  #this year will be used in all dates of monthly climatologies (and day will = 15)
196

  
197
## Monthly means
198
system(paste("cdo -O sorttimestamp -setyear,",myear," -setday,15 -ymonmean ",outdir2,"/MOD06_",tile,"_daily.nc ",tsdir,"/MOD06_",tile,"_ymonmean.nc",sep=""))
199

  
200
## Monthly standard deviation
201
system(paste("cdo -O sorttimestamp -setyear,",myear," -setday,15 -ymonstd ",outdir2,"/MOD06_",tile,"_daily.nc ",tsdir,"/MOD06_",tile,"_ymonstd.nc",sep=""))
202
system(paste(ncopath,"ncrename -v CER,CER_sd -v CLD,CLD_sd -v COT,COT_sd ",tsdir,"/MOD06_",tile,"_ymonstd.nc",sep=""))
203
system(paste(ncopath,"ncatted ",
204
" -a long_name,CER_sd,o,c,\"Cloud Particle Effective Radius (standard deviation of daily observations)\" ",
205
" -a long_name,CLD_sd,o,c,\"Cloud Mask (standard deviation of daily observations)\" ",
206
" -a long_name,COT_sd,o,c,\"Cloud Optical Thickness (standard deviation of daily observations)\" ",
207
tsdir,"/MOD06_",tile,"_ymonstd.nc",sep=""))
208

  
209
## cer > 20
210
system(paste("cdo -O  sorttimestamp -setyear,",myear," -setday,15 -ymonmean -gtc,20 -selvar,CER ",outdir2,"/MOD06_",tile,"_daily.nc ",tsdir,"/MOD06_",tile,"_ymoncer20.nc",sep=""))
211
system(paste(ncopath,"ncrename -v CER,CER20 ",tsdir,"/MOD06_",tile,"_ymoncer20.nc",sep=""))
212
system(paste(ncopath,"ncatted ",
213
" -a long_name,CER20,o,c,\"Proportion of Days with Cloud Particle Effective Radius > 20um\" ",
214
" -a units,CER20,o,c,\"Proportion\" ",
215
tsdir,"/MOD06_",tile,"_ymoncer20.nc",sep=""))
216

  
217
## number of observations
218
system(paste("cdo -O sorttimestamp  -setyear,",myear," -setday,15 -nint -mulc,100 -ymonmean -eqc,9999 -setmisstoc,9999   -selvar,CER,CLD ",outdir2,"/MOD06_",tile,"_daily.nc ",tsdir,"/MOD06_",tile,"_ymonmiss.nc",sep=""))
219
system(paste(ncopath,"ncrename -v CER,CER_pmiss -v CLD,CLD_pmiss ",tsdir,"/MOD06_",tile,"_ymonmiss.nc",sep=""))
220
system(paste(ncopath,"ncatted ",
221
" -a long_name,CER_pmiss,o,c,\"Proportion of Days with missing data for CER\" ",
222
" -a long_name,CLD_pmiss,o,c,\"Proportion of Days with missing data for CLD\" ",
223
" -a scale_factor,CER_pmiss,o,d,0.01 ",
224
" -a units,CER_pmiss,o,c,\"Proportion\" ",
225
" -a scale_factor,CLD_pmiss,o,d,0.01 ",
226
" -a units,CLD_pmiss,o,c,\"Proportion\" ",
227
tsdir,"/MOD06_",tile,"_ymonmiss.nc",sep=""))
228

  
229
## append variables to a single file
230
system(paste(ncopath,"ncks -O ",tsdir,"/MOD06_",tile,"_ymonmean.nc  ",tsdir,"/MOD06_",tile,"_ymon.nc",sep=""))
231
system(paste(ncopath,"ncks -A ",tsdir,"/MOD06_",tile,"_ymonstd.nc  ",tsdir,"/MOD06_",tile,"_ymon.nc",sep=""))
232
system(paste(ncopath,"ncks -A ",tsdir,"/MOD06_",tile,"_ymoncer20.nc  ",tsdir,"/MOD06_",tile,"_ymon.nc",sep=""))
233
system(paste(ncopath,"ncks -A ",tsdir,"/MOD06_",tile,"_ymonmiss.nc  ",tsdir,"/MOD06_",tile,"_ymon.nc",sep=""))
234

  
235
## append sinusoidal grid from one of input files as CDO doesn't transfer all attributes
236
system(paste(ncopath,"ncks -A -v sinusoidal ",list.files(tsdir,full=T,pattern="daily[.]nc$")[1],"  ",tsdir,"/MOD06_",tile,"_ymon.nc",sep=""))
237

  
238
## invert latitude so it plays nicely with gdal
239
system(paste(ncopath,"ncpdq -O -a -y ",tsdir,"/MOD06_",tile,"_ymon.nc ",outdir2,"/MOD06_",tile,".nc",sep=""))
240

  
241
## update attributes
242
system(paste(ncopath,"ncatted ",
243
#" -a standard_parallel,sinusoidal,o,c,\"0\" ",
244
#" -a longitude_of_central_meridian,sinusoidal,o,c,\"0\" ",
245
#" -a latitude_of_central_meridian,sinusoidal,o,c,\"0\" ",
246
" -a units,time,o,c,\"days since 2000-1-1 0:0:0\" ",
247
" -a title,global,o,c,\"MODIS Cloud Product (MOD06) Climatology\" ",
248
" -a institution,global,o,c,\"Yale University\" ",
249
" -a source,global,o,c,\"MODIS Cloud Product (MOD06)\" ",
250
" -a comment,global,o,c,\"Compiled by Adam M. Wilson (adam.wilson@yale.edu)\" ",
251
outdir2,"/MOD06_",tile,".nc",sep=""))
252

  
200 253

  
201 254
print("Finished!   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
202 255
## quit R
......
206 259
#################################################################
207 260

  
208 261
### copy the files back to Yale
209
list.files("2_daily")
210
system(paste("scp ",outdir2,"/*_ymonmean.nc adamw@acrobates.eeb.yale.edu:/data/personal/adamw/projects/interp/data/modis/Venezuela/summary",sep=""))
262
system(paste("scp ",outdir2,"/MOD06_",tile,".nc adamw@acrobates.eeb.yale.edu:/data/personal/adamw/projects/interp/data/modis/mod06/summary",sep=""))
211 263

  
212 264
system("scp  /tmp/Rtmp6I6tFn/MOD06_L2.A2000061.1615.051.2010273184629.hdf adamw@acrobates.eeb.yale.edu:/data/personal/adamw/projects/interp/data/modis/Venezuela")
213 265
system("scp 2_daily/MOD06_20000410.nc adamw@acrobates.eeb.yale.edu:/data/personal/adamw/projects/interp/data/modis/Venezuela")

Also available in: Unified diff