Project

General

Profile

« Previous | Next » 

Revision 3133718a

Added by Adam Wilson over 11 years ago

Added new script to process daily MOD35 files

View differences:

climate/procedures/MOD35_Explore.r
4 4
library(raster)
5 5
library(rgdal)
6 6

  
7
f=list.files(pattern="Cloud_Mask_1_1")
7
f=list.files(pattern="*.hdf")
8

  
9
Sys.setenv(GEOL_AS_GCPS = "PARTIAL")
8 10

  
9 11
GDALinfo(f[1])
12
system(paste("gdalinfo",f[1]))
13
GDALinfo("HDF4_EOS:EOS_SWATH:\"MOD35_L2.A2000100.1445.006.2012252024758.hdf\":mod35:Cloud_Mask")
14
system("gdalinfo HDF4_EOS:EOS_SWATH:\"MOD35_L2.A2000100.1445.006.2012252024758.hdf\":mod35:Cloud_Mask | tail -n 200")
15

  
16
system("gdalwarp -overwrite -geoloc -order 2 -r near -s_srs \"EPSG:4326\" HDF4_EOS:EOS_SWATH:\"MOD35_L2.A2000100.1445.006.2012252024758.hdf\":mod35:Cloud_Mask cloudmask.tif")
17
system("gdalwarp -overwrite -r near -s_srs \"EPSG:4326\" HDF4_EOS:EOS_SWATH:\"MOD35_L2.A2000100.1445.006.2012252024758.hdf\":mod35:Cloud_Mask:1 cloudmask2.tif")
18

  
10 19

  
11 20
## get tile
12 21
tile=raster("~/acrobates/projects/interp/data/modis/mod06/summary/MOD06_h09v04.nc",varname="CER")
climate/procedures/MOD35_L2_process.r
1
###################################################################################
2
###  R code to aquire and process MOD35_L2 cloud data from the MODIS platform
3

  
4

  
5
# Redirect all warnings to stderr()
6
#options(warn = -1)
7
#write("2) write() to stderr", stderr())
8
#write("2) write() to stdout", stdout())
9
#warning("2) warning()")
10

  
11

  
12
## import commandline arguments
13
library(getopt)
14
## load libraries
15
require(reshape)
16
require(geosphere)
17
require(raster)
18
require(rgdal)
19
require(spgrass6)
20
require(RSQLite)
21

  
22
## get options
23
opta <- getopt(matrix(c(
24
                        'date', 'd', 1, 'character',
25
                        'tile', 't', 1, 'character',
26
                        'verbose','v',1,'logical',
27
                        'help', 'h', 0, 'logical'
28
                        ), ncol=4, byrow=TRUE))
29
if ( !is.null(opta$help) )
30
  {
31
       prg <- commandArgs()[1];
32
          cat(paste("Usage: ", prg,  " --date | -d <file> :: The date to process\n", sep=""));
33
          q(status=1);
34
     }
35

  
36

  
37
## default date and tile to play with  (will be overwritten below when running in batch)
38
date="20010410"
39
tile="h11v08"
40
platform="pleiades" 
41
verbose=T
42

  
43
## now update using options if given
44
date=opta$date  
45
tile=opta$tile 
46
verbose=opta$verbose  #print out extensive information for debugging?
47
outdir=paste("daily/",tile,"/",sep="")  #directory for separate daily files
48
## get year and doy from date
49
year=format(as.Date(date,"%Y%m%d"),"%Y")
50
doy=format(as.Date(date,"%Y%m%d"),"%j")
51

  
52
if(platform=="pleiades"){
53
  ## location of MOD06 files
54
  datadir=paste("/nobackupp4/datapool/modis/MOD06_L2.005/",year,"/",doy,"/",sep="")
55
  ## path to some executables
56
  ncopath="/nasa/sles11/nco/4.0.8/gcc/mpt/bin/"
57
  swtifpath="/nobackupp1/awilso10/software/heg/bin/swtif"
58
  ## path to swath database
59
  db="/nobackupp4/pvotava/DB/export/swath_geo.sql.sqlite3.db"
60
  ## specify working directory
61
  setwd("/nobackupp1/awilso10/mod06")
62
  gisBase="/u/armichae/pr/grass-6.4.2/"
63
  ## path to MOD11A1 file for this tile to align grid/extent
64
  gridfile=list.files("/nobackupp4/datapool/modis/MOD11A1.005/2006.01.27",pattern=paste(tile,".*[.]hdf$",sep=""),recursive=T,full=T)[1]
65
  td=readGDAL(paste("HDF4_EOS:EOS_GRID:\"",gridfile,"\":MODIS_Grid_Daily_1km_LST:Night_view_angl",sep=""))
66
  projection(td)="+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs +datum=WGS84 +ellps=WGS84 "
67
}
68

  
69
if(platform=="litoria"){  #if running on local server, use different paths
70
  ## specify working directory
71
  setwd("~/acrobates/projects/interp")
72
  gisBase="/usr/lib/grass64"
73
   ## location of MOD06 files
74
  datadir="~/acrobates/projects/interp/data/modis/mod35"
75
  ## path to some executables
76
  ncopath=""
77
  swtifpath="sudo MRTDATADIR=\"/usr/local/heg/data\" PGSHOME=/usr/local/heg/TOOLKIT_MTD PWD=/home/adamw /usr/local/heg/bin/swtif"
78
  ## path to swath database
79
  db="/home/adamw/acrobates/projects/interp/data/modis/mod06/swath_geo.sql.sqlite3.db"
80
  ## get grid file
81
  td=raster(paste("~/acrobates/projects/interp/data/modis/mod06/summary/MOD06_",tile,".nc",sep=""),varname="CER")
82
  projection(td)="+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs +datum=WGS84 +ellps=WGS84 "
83
}
84

  
85

  
86
### print some status messages
87
if(verbose) print(paste("Processing tile",tile," for date",date))
88

  
89
## load tile information and get bounding box
90
load(file="modlandTiles.Rdata")
91
tile_bb=tb[tb$tile==tile,] ## identify tile of interest
92
upleft=paste(tile_bb$lat_max,tile_bb$lon_min) #northwest corner
93
lowright=paste(tile_bb$lat_min,tile_bb$lon_max) #southeast corner
94

  
95

  
96
## vars to process
97
vars=as.data.frame(matrix(c(
98
  "Cloud_Mask",              "CM",
99
  "Quality_Assurance",       "QA"),
100
  byrow=T,ncol=2,dimnames=list(1:2,c("variable","varid"))),stringsAsFactors=F)
101

  
102
## vector of variables expected to be in final netcdf file.  If these are not present, the file will be deleted at the end.
103
finalvars=c("CLD","CLD2")
104

  
105

  
106
#####################################################
107
##find swaths in region from sqlite database for the specified date/tile
108
if(verbose) print("Accessing swath ID's from database")
109
con=dbConnect("SQLite", dbname = db)
110
fs=dbGetQuery(con,paste("SELECT * from swath_geo
111
            WHERE east>=",tile_bb$lon_min," AND
112
                  west<=",tile_bb$lon_max," AND
113
                  north>=",tile_bb$lat_min," AND
114
                  south<=",tile_bb$lat_max," AND
115
                  year==",format(as.Date(date,"%Y%m%d"),"%Y")," AND
116
                  day==",as.numeric(format(as.Date(date,"%Y%m%d"),"%j"))
117
  ))
118
con=dbDisconnect(con)
119
fs$id=substr(fs$id,7,19)
120
## find the swaths on disk (using datadir)
121
swaths=list.files(datadir,pattern=paste(fs$id,collapse="|"),recursive=T,full=T)
122

  
123
if(verbose) print(paste(nrow(fs)," swath IDs recieved from database and ",length(swaths)," found on disk"))
124

  
125

  
126
############################################################################
127
############################################################################
128
### Use the HEG tool to grid all available swath data for this date-tile
129
for(file in swaths){
130
  ## Function to generate hegtool parameter file for multi-band HDF-EOS file
131
  print(paste("Starting file",basename(file)))
132
  outfile=paste(tempdir(),"/",basename(file),sep="")
133
  ## First write the parameter file (careful, heg is very finicky!)
134
  hdr=paste("NUM_RUNS = ",length(vars$varid),"|MULTI_BAND_HDFEOS:",length(vars$varid),sep="")
135
  grp=paste("
136
BEGIN
137
INPUT_FILENAME=",file,"
138
OBJECT_NAME=mod35
139
FIELD_NAME=",vars$variable,"|
140
BAND_NUMBER = 1
141
OUTPUT_PIXEL_SIZE_X=1000
142
OUTPUT_PIXEL_SIZE_Y=1000
143
SPATIAL_SUBSET_UL_CORNER = ( ",upleft," )
144
SPATIAL_SUBSET_LR_CORNER = ( ",lowright," )
145
#RESAMPLING_TYPE =",ifelse(grepl("Flag|Mask|Quality",vars),"NN","CUBIC"),"
146
RESAMPLING_TYPE =NN
147
OUTPUT_PROJECTION_TYPE = SIN
148
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 )
149
# projection parameters from http://landweb.nascom.nasa.gov/cgi-bin/QA_WWW/newPage.cgi?fileName=sn_gctp
150
ELLIPSOID_CODE = WGS84
151
OUTPUT_TYPE = HDFEOS
152
OUTPUT_FILENAME = ",outfile,"
153
END
154

  
155
",sep="")
156

  
157
  ## if any remnants from previous runs remain, delete them
158
  if(length(list.files(tempdir(),pattern=basename(file)))>0)
159
    file.remove(list.files(tempdir(),pattern=basename(file),full=T))
160
  ## write it to a file
161
  cat(c(hdr,grp)    , file=paste(tempdir(),"/",basename(file),"_MODparms.txt",sep=""))
162
  ## now run the swath2grid tool
163
  ## write the gridded file
164
  log=system(paste("(",swtifpath," -p ",tempdir(),"/",basename(file),"_MODparms.txt -d ; echo $$)",sep=""),intern=T,ignore.stderr=T)
165
  ## clean up temporary files in working directory
166
#  file.remove(list.files(pattern=
167
#              paste("filetable.temp_",
168
#              as.numeric(log[length(log)]):(as.numeric(log[length(log)])+3),sep="",collapse="|")))  #Look for files with PID within 3 of parent process
169
  if(verbose) print(log)
170
  print(paste("Finished gridding ", file))
171
}  #end looping over swaths
172

  
173
########################
174
## confirm at least one file for this date is present.  If not, quit.
175
outfiles=paste(tempdir(),"/",basename(swaths),sep="")
176
if(!any(file.exists(outfiles))) {
177
  print(paste("########################################   No gridded files for region exist for tile",tile," on date",date))
178
  q("no",status=0)
179
}
180

  
181
#####################################################
182
## Process the gridded files to align exactly with MODLAND tile and produce a daily summary of multiple swaths
183
  
184
## Identify output file
185
  ncfile=paste(outdir,"/MOD06_",tile,"_",date,".nc",sep="")  #this is the 'final' daily output file
186

  
187
## function to convert binary to decimal to assist in identifying correct values
188
## this is helpful when defining QA handling below, but isn't used in processing
189
## b2d=function(x) sum(x * 2^(rev(seq_along(x)) - 1)) #http://tolstoy.newcastle.edu.au/R/e2/help/07/02/10596.html
190
## for example:
191
## b2d(c(T,T))
192

  
193
  ## set Grass to overwrite
194
  Sys.setenv(GRASS_OVERWRITE=1)
195
  Sys.setenv(DEBUG=1)
196
  Sys.setenv(GRASS_GUI="txt")
197

  
198
### Extract various SDSs from a single gridded HDF file and use QA data to throw out 'bad' observations
199
## make temporary working directory
200
  tf=paste(tempdir(),"/grass", Sys.getpid(),"/", sep="")  #temporar
201
  if(!file.exists(tf)) dir.create(tf)
202
  ## create output directory if needed
203
  if(!file.exists(dirname(ncfile))) dir.create(dirname(ncfile),recursive=T)
204
  
205
  ## set up temporary grass instance for this PID
206
  if(verbose) print(paste("Set up temporary grass session in",tf))
207
  initGRASS(gisBase=gisBase,gisDbase=tf,SG=as(td,"SpatialGridDataFrame"),override=T,location="mod06",mapset="PERMANENT",home=tf,pid=Sys.getpid())
208
  system(paste("g.proj -c proj4=\"",projection(td),"\"",sep=""),ignore.stdout=T,ignore.stderr=T)
209

  
210
  ## Define region by importing one MOD11A1 raster.
211
  print("Import one MOD11A1 raster to define grid")
212
  if(platform=="pleiades") execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",gridfile,"\":MODIS_Grid_Daily_1km_LST:Night_view_angl",sep=""),
213
            output="modisgrid",flags=c("quiet","overwrite","o"))
214
   if(platform=="litoria")
215
  execGRASS("r.in.gdal",input=paste("NETCDF:\"/home/adamw/acrobates/projects/interp/data/modis/mod06/summary/MOD06_",tile,".nc\":CER",sep=""),output="modisgrid",flags=c("overwrite","o"))
216
  
217
system("g.region rast=modisgrid save=roi --overwrite",ignore.stdout=F,ignore.stderr=F)
218
system("g.region rast=modisgrid.1 save=roi --overwrite",ignore.stdout=F,ignore.stderr=F)
219

  
220
## Identify which files to process
221
tfs=basename(swaths)
222
## drop swaths that did not produce an output file (typically due to not overlapping the ROI)
223
tfs=tfs[tfs%in%list.files(tempdir())]
224
#tfs=list.files(tempdir(),pattern="temp.*hdf")
225
nfs=length(tfs)
226
if(verbose) print(paste(nfs,"swaths available for processing"))
227

  
228
## loop through scenes and process QA flags
229
  for(i in 1:nfs){
230
     file=paste(tempdir(),"/",tfs[i],sep="")
231
     ## Cloud Mask
232
     execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod35:Cloud_Mask_0",sep=""),
233
              output=paste("CM1_",i,sep=""),flags=c("overwrite","o")) ; print("")
234
    ## extract cloudy and 'probably/confidently clear' pixels
235
    system(paste("r.mapcalc <<EOF
236
                CM_cloud_",i," =  ((CM1_",i," / 2^0) % 2) == 1  &&  ((CM1_",i," / 2^1) % 2^2) == 0 
237
                CM_clear_",i," =  ((CM1_",i," / 2^0) % 2) == 1  &&  ((CM1_",i," / 2^1) % 2^2) >=  2 
238
                CM_path_",i," =   ((CM1_",i," / 2^6) % 2^2) 
239
                CM_cloud2_",i," = ((CM1_",i," / 2^1) % 2^2) 
240
EOF",sep=""))
241

  
242
    ## QA
243
     execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod35:Quality_Assurance_0",sep=""),
244
             output=paste("QA_",i,sep=""),flags=c("overwrite","o")) ; print("")
245
   ## QA_CER
246
#   system(paste("r.mapcalc <<EOF
247
#                 QA_COT_",i,"=   ((QA_",i," / 2^0) % 2^1 )==1
248
#                 QA_COT2_",i,"=  ((QA_",i," / 2^1) % 2^2 )>=2
249
#                 QA_COT3_",i,"=  ((QA_",i," / 2^3) % 2^2 )==0
250
#                 QA_CER_",i,"=   ((QA_",i," / 2^5) % 2^1 )==1
251
#                 QA_CER2_",i,"=  ((QA_",i," / 2^6) % 2^2 )>=2
252
#EOF",sep="")) 
253
#                 QA_CWP_",i,"=   ((QA_",i," / 2^8) % 2^1 )==1
254
#                 QA_CWP2_",i,"=  ((QA_",i," / 2^9) % 2^2 )==3
255

  
256
     
257
 } #end loop through sub daily files
258

  
259
#### Now generate daily averages (or maximum in case of cloud flag)
260
  
261
  system(paste("r.mapcalc <<EOF
262
         CLD_daily=int((max(",paste("if(isnull(CM_cloud_",1:nfs,"),-9999,CM_cloud_",1:nfs,")",sep="",collapse=","),"))) 
263
         CLD2_daily=int((min(",paste("if(isnull(CM_cloud2_",1:nfs,"),-9999,CM_cloud2_",1:nfs,")",sep="",collapse=","),"))) 
264
EOF",sep=""))
265

  
266
execGRASS("r.null",map="CLD_daily",setnull="-9999")
267
execGRASS("r.null",map="CLD2_daily",setnull="-9999")
268

  
269

  
270
  ### Write the files to a netcdf file
271
  ## create image group to facilitate export as multiband netcdf
272
    execGRASS("i.group",group="mod06",input=c("CLD_daily","CLD2_daily")) ; print("")
273
   
274
  if(file.exists(ncfile)) file.remove(ncfile)  #if it exists already, delete it
275
  execGRASS("r.out.gdal",input="mod06",output=ncfile,type="Int16",nodata=-32768,flags=c("quiet"),
276
#      createopt=c("FORMAT=NC4","ZLEVEL=5","COMPRESS=DEFLATE","WRITE_GDAL_TAGS=YES","WRITE_LONLAT=NO"),format="netCDF")  #for compressed netcdf
277
      createopt=c("FORMAT=NC","WRITE_GDAL_TAGS=YES","WRITE_LONLAT=NO"),format="netCDF")
278

  
279
  system(paste(ncopath,"ncecat -O -u time ",ncfile," ",ncfile,sep=""))
280
## create temporary nc file with time information to append to MOD06 data
281
  cat(paste("
282
    netcdf time {
283
      dimensions:
284
        time = 1 ;
285
      variables:
286
        int time(time) ;
287
      time:units = \"days since 2000-01-01 00:00:00\" ;
288
      time:calendar = \"gregorian\";
289
      time:long_name = \"time of observation\"; 
290
    data:
291
      time=",as.integer(as.Date(date,"%Y%m%d")-as.Date("2000-01-01")),";
292
    }"),file=paste(tempdir(),"/time.cdl",sep=""))
293
system(paste("ncgen -o ",tempdir(),"/time.nc ",tempdir(),"/time.cdl",sep=""))
294
system(paste(ncopath,"ncks -A ",tempdir(),"/time.nc ",ncfile,sep=""))
295
## add other attributes
296
  system(paste(ncopath,"ncrename -v Band1,CLD -v Band2,CLD2 ",ncfile,sep=""))
297
  system(paste(ncopath,"ncatted -a scale_factor,CLD,o,d,1 -a units,CLD,o,c,\"none\" -a missing_value,CLD,o,d,-32768 -a long_name,CLD,o,c,\"Cloud Mask\" ",ncfile,sep=""))
298
system(paste(ncopath,"ncatted -a scale_factor,CLD2,o,d,1 -a units,CLD2,o,c,\"none\" -a missing_value,CLD2,o,d,-32768 -a long_name,CLD2,o,c,\"Cloud Mask Flag\" ",ncfile,sep=""))
299

  
300
                                        #  system(paste(ncopath,"ncatted -a sourcecode,global,o,c,",script," ",ncfile,sep=""))
301
   
302
  
303
### delete the temporary files 
304
  unlink_.gislock()
305
  system(paste("rm -frR ",tf,sep=""))
306

  
307

  
308
## Confirm that the file has the correct attributes, otherwise delete it
309
ntime=as.numeric(system(paste("cdo -s ntime ",ncfile),intern=T))
310
## confirm it has all 'final variables as specified above"
311
fvar=all(finalvars%in%strsplit(system(paste("cdo -s showvar ",ncfile),intern=T)," ")[[1]])
312

  
313
  if(ntime!=1|!fvar) {
314
      print(paste("FILE ERROR:  tile ",tile," and date ",date," was not outputted correctly, deleting... "))
315
      file.remove(ncfile)
316
    }
317
    
318
  ## print out some info
319
print(paste(" ###################################################################               Finished ",date,
320
"################################################################"))
321

  
322
## delete old files
323
system("cleartemp")
324

  
325
## quit
326
q("no",status=0)

Also available in: Unified diff