Project

General

Profile

« Previous | Next » 

Revision f27d277d

Added by Adam Wilson over 11 years ago

Simplified MOD35 script and converted output to single band p(Clear) rather than multiple bands

View differences:

climate/procedures/MOD35_L2_process.r
35 35

  
36 36

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

  
69 69
if(platform=="litoria"){  #if running on local server, use different paths
70 70
  ## specify working directory
71
  setwd("~/acrobates/projects/interp")
71
  setwd("~/acrobates/adamw/projects/interp")
72 72
  gisBase="/usr/lib/grass64"
73 73
   ## location of MOD06 files
74
  datadir="~/acrobates/projects/interp/data/modis/mod35"
74
  datadir="~/acrobates/adamw/projects/interp/data/modis/mod35"
75 75
  ## path to some executables
76 76
  ncopath=""
77 77
  swtifpath="sudo MRTDATADIR=\"/usr/local/heg/data\" PGSHOME=/usr/local/heg/TOOLKIT_MTD PWD=/home/adamw /usr/local/heg/bin/swtif"
78 78
  ## path to swath database
79
  db="/home/adamw/acrobates/projects/interp/data/modis/mod06/swath_geo.sql.sqlite3.db"
79
  db="~/acrobates/adamw/projects/interp/data/modis/mod06/swath_geo.sql.sqlite3.db"
80 80
  ## get grid file
81
  td=raster(paste("~/acrobates/projects/interp/data/modis/mod06/summary/MOD06_",tile,".nc",sep=""),varname="CER")
81
  td=raster(paste("~/acrobates/adamw/projects/interp/data/modis/mod06/summary/MOD06_",tile,".nc",sep=""),varname="CER")
82 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 83
}
84 84

  
......
100 100
  byrow=T,ncol=2,dimnames=list(1:2,c("variable","varid"))),stringsAsFactors=F)
101 101

  
102 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")
103
finalvars=c("PClear")
104 104

  
105 105

  
106 106
#####################################################
107 107
##find swaths in region from sqlite database for the specified date/tile
108 108
if(verbose) print("Accessing swath ID's from database")
109 109
con=dbConnect("SQLite", dbname = db)
110
fs=dbGetQuery(con,paste("SELECT * from swath_geo
110
fs=dbGetQuery(con,paste("SELECT * from swath_geo6
111 111
            WHERE east>=",tile_bb$lon_min," AND
112 112
                  west<=",tile_bb$lon_max," AND
113 113
                  north>=",tile_bb$lat_min," AND
......
209 209

  
210 210
  ## Define region by importing one MOD11A1 raster.
211 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

  
212
  if(platform=="pleiades") {
213
    execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",gridfile,"\":MODIS_Grid_Daily_1km_LST:Night_view_angl",sep=""),
214
              output="modisgrid",flags=c("quiet","overwrite","o"))
215
    system("g.region rast=modisgrid save=roi --overwrite",ignore.stdout=F,ignore.stderr=F)
216
  }
217

  
218
if(platform=="litoria"){
219
  execGRASS("r.in.gdal",input=paste("NETCDF:\"/home/adamw/acrobates/adamw/projects/interp/data/modis/mod06/summary/MOD06_",tile,".nc\":CER",sep=""),
220
            output="modisgrid",flags=c("overwrite","o"))
221
  system("g.region rast=modisgrid.1 save=roi --overwrite",ignore.stdout=F,ignore.stderr=F)
222
}
220 223
## Identify which files to process
221 224
tfs=basename(swaths)
222 225
## drop swaths that did not produce an output file (typically due to not overlapping the ROI)
......
231 234
     ## Cloud Mask
232 235
     execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod35:Cloud_Mask_0",sep=""),
233 236
              output=paste("CM1_",i,sep=""),flags=c("overwrite","o")) ; print("")
234
    ## extract cloudy and 'probably/confidently clear' pixels
237
    ## QA      ## extract first bit to keep only "useful" values of cloud mask
238
     execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod35:Quality_Assurance_0",sep=""),
239
             output=paste("QA_",i,sep=""),flags=c("overwrite","o")) ; print("")
240
     ## extract first two bits of cloud mask
235 241
    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 
242
                QA_useful_",i," =  if((QA_",i," / 2^0) % 2==1,1,0)
243
                CM_cloud_",i," =  if((CM1_",i," / 2^0) % 2==1,(CM1_",i," / 2^1) % 2^2,-9999)
244
                Pclear1_",i," = if(CM_cloud_",i,"==0,0,if(CM_cloud_",i,"==1,66,if(CM_cloud_",i,"==2,95,if(CM_cloud_",i,"==3,99,-9999))))
245
                Pclear_",i," = if(QA_useful_",i,"==1,Pclear1_",i,",-9999)
238 246
                CM_path_",i," =   ((CM1_",i," / 2^6) % 2^2) 
239
                CM_cloud2_",i," = ((CM1_",i," / 2^1) % 2^2) 
240 247
EOF",sep=""))
241 248

  
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
     
249
     execGRASS("r.null",map=paste("Pclear_",i,sep=""),setnull="-9999")
250
     execGRASS("r.null",map=paste("CM_cloud_",i,sep=""),setnull="-9999")
251
   
252
    
257 253
 } #end loop through sub daily files
258 254

  
259 255
#### Now generate daily averages (or maximum in case of cloud flag)
260 256
  
261 257
  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=","),"))) 
258
         Pclear_daily=int((min(",paste("if(isnull(Pclear_",1:nfs,"),9999,Pclear_",1:nfs,")",sep="",collapse=","),"))) 
264 259
EOF",sep=""))
265 260

  
266
execGRASS("r.null",map="CLD_daily",setnull="-9999")
267
execGRASS("r.null",map="CLD2_daily",setnull="-9999")
261
#         CLD_daily=int((min(",paste("if(isnull(CM_cloud_",1:nfs,"),9999,CM_cloud_",1:nfs,")",sep="",collapse=","),"))) 
262

  
263
## reset null values
264
#execGRASS("r.null",map="CLD_daily",setnull="9999")
265
execGRASS("r.null",map="Pclear_daily",setnull="9999")
268 266

  
269 267

  
270 268
  ### Write the files to a netcdf file
271 269
  ## create image group to facilitate export as multiband netcdf
272
    execGRASS("i.group",group="mod06",input=c("CLD_daily","CLD2_daily")) ; print("")
270
    execGRASS("i.group",group="mod35",input=c("Pclear_daily")) ; print("")
273 271
   
274 272
  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"),
273
  execGRASS("r.out.gdal",input="mod35",output=ncfile,type="Byte",nodata=255,flags=c("quiet"),
276 274
#      createopt=c("FORMAT=NC4","ZLEVEL=5","COMPRESS=DEFLATE","WRITE_GDAL_TAGS=YES","WRITE_LONLAT=NO"),format="netCDF")  #for compressed netcdf
277 275
      createopt=c("FORMAT=NC","WRITE_GDAL_TAGS=YES","WRITE_LONLAT=NO"),format="netCDF")
278 276

  
......
293 291
system(paste("ncgen -o ",tempdir(),"/time.nc ",tempdir(),"/time.cdl",sep=""))
294 292
system(paste(ncopath,"ncks -A ",tempdir(),"/time.nc ",ncfile,sep=""))
295 293
## 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=""))
294
  system(paste(ncopath,"ncrename -v Band1,PClear ",ncfile,sep=""))
295
  system(paste(ncopath,"ncatted -a scale_factor,PClear,o,d,1 -a units,PClear,o,c,\"Probability (%)\" -a missing_value,PClear,o,d,255 -a long_name,PClear,o,c,\"Probability of Clear Sky\" ",ncfile,sep=""))
299 296

  
300 297
                                        #  system(paste(ncopath,"ncatted -a sourcecode,global,o,c,",script," ",ncfile,sep=""))
301 298
   
......
320 317
"################################################################"))
321 318

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

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

Also available in: Unified diff