| 1 |
d2dbbe0c
|
Adam M. Wilson
|
###################################################################################
|
| 2 |
|
|
### R code to aquire and process MOD06_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="20030102"
|
| 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 |
|
|
|
| 49 |
|
|
|
| 50 |
|
|
## get year and doy from date
|
| 51 |
|
|
year=format(as.Date(date,"%Y%m%d"),"%Y")
|
| 52 |
|
|
doy=format(as.Date(date,"%Y%m%d"),"%j")
|
| 53 |
|
|
|
| 54 |
|
|
if(platform=="pleiades"){
|
| 55 |
|
|
## location of MOD06 files
|
| 56 |
|
|
datadir=paste("/nobackupp4/datapool/modis/MOD06_L2.005/",year,"/",doy,"/",sep="")
|
| 57 |
|
|
## path to some executables
|
| 58 |
|
|
ncopath="/nasa/sles11/nco/4.0.8/gcc/mpt/bin/"
|
| 59 |
|
|
swtifpath="/nobackupp1/awilso10/software/heg/bin/swtif"
|
| 60 |
|
|
## path to swath database
|
| 61 |
|
|
db="/nobackupp4/pvotava/DB/export/swath_geo.sql.sqlite3.db"
|
| 62 |
|
|
## specify working directory
|
| 63 |
|
|
setwd("/nobackupp1/awilso10/mod06")
|
| 64 |
|
|
gisBase="/u/armichae/pr/grass-6.4.2/"
|
| 65 |
|
|
## path to MOD11A1 file for this tile to align grid/extent
|
| 66 |
|
|
gridfile=list.files("/nobackupp4/datapool/modis/MOD11A1.005/2006.01.27",pattern=paste(tile,".*[.]hdf$",sep=""),recursive=T,full=T)[1]
|
| 67 |
|
|
td=readGDAL(paste("HDF4_EOS:EOS_GRID:\"",gridfile,"\":MODIS_Grid_Daily_1km_LST:Night_view_angl",sep=""))
|
| 68 |
|
|
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 "
|
| 69 |
|
|
}
|
| 70 |
|
|
|
| 71 |
|
|
if(platform=="litoria"){ #if running on local server, use different paths
|
| 72 |
|
|
## specify working directory
|
| 73 |
|
|
setwd("~/acrobates/projects/interp")
|
| 74 |
|
|
gisBase="/usr/lib/grass64"
|
| 75 |
|
|
## location of MOD06 files
|
| 76 |
|
|
datadir="~/acrobates/projects/interp/data/modis/Venezuela/MOD06"
|
| 77 |
|
|
## path to some executables
|
| 78 |
|
|
ncopath=""
|
| 79 |
|
|
swtifpath="sudo MRTDATADIR=\"/usr/local/heg/data\" PGSHOME=/usr/local/heg/TOOLKIT_MTD PWD=/home/adamw /usr/local/heg/bin/swtif"
|
| 80 |
|
|
## path to swath database
|
| 81 |
|
|
db="/home/adamw/acrobates/projects/interp/data/modis/mod06/swath_geo.sql.sqlite3.db"
|
| 82 |
|
|
## get grid file
|
| 83 |
|
|
td=raster(paste("~/acrobates/projects/interp/data/modis/mod06/summary/MOD06_",tile,".nc",sep=""),varname="CER")
|
| 84 |
|
|
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 "
|
| 85 |
|
|
}
|
| 86 |
|
|
|
| 87 |
|
|
|
| 88 |
|
|
### print some status messages
|
| 89 |
|
|
if(verbose) print(paste("Processing tile",tile," for date",date))
|
| 90 |
|
|
|
| 91 |
|
|
## load tile information and get bounding box
|
| 92 |
|
|
load(file="modlandTiles.Rdata")
|
| 93 |
|
|
tile_bb=tb[tb$tile==tile,] ## identify tile of interest
|
| 94 |
|
|
upleft=paste(tile_bb$lat_max,tile_bb$lon_min) #northwest corner
|
| 95 |
|
|
lowright=paste(tile_bb$lat_min,tile_bb$lon_max) #southeast corner
|
| 96 |
|
|
|
| 97 |
|
|
|
| 98 |
|
|
## vars to process
|
| 99 |
|
|
vars=as.data.frame(matrix(c(
|
| 100 |
|
|
"Cloud_Effective_Radius", "CER",
|
| 101 |
|
|
"Cloud_Effective_Radius_Uncertainty", "CERU",
|
| 102 |
|
|
"Cloud_Optical_Thickness", "COT",
|
| 103 |
|
|
"Cloud_Optical_Thickness_Uncertainty", "COTU",
|
| 104 |
|
|
# "Cloud_Water_Path", "CWP",
|
| 105 |
|
|
# "Cloud_Water_Path_Uncertainty", "CWPU",
|
| 106 |
|
|
# "Cloud_Phase_Optical_Properties", "CPOP",
|
| 107 |
|
|
# "Cloud_Multi_Layer_Flag", "CMLF",
|
| 108 |
|
|
"Cloud_Mask_1km", "CM1",
|
| 109 |
|
|
"Quality_Assurance_1km", "QA"),
|
| 110 |
|
|
byrow=T,ncol=2,dimnames=list(1:6,c("variable","varid"))),stringsAsFactors=F)
|
| 111 |
|
|
|
| 112 |
|
|
## vector of variables expected to be in final netcdf file. If these are not present, the file will be deleted at the end.
|
| 113 |
|
|
finalvars=c("CER","COT","CLD")
|
| 114 |
|
|
|
| 115 |
|
|
|
| 116 |
|
|
#####################################################
|
| 117 |
|
|
##find swaths in region from sqlite database for the specified date/tile
|
| 118 |
|
|
if(verbose) print("Accessing swath ID's from database")
|
| 119 |
|
|
con=dbConnect("SQLite", dbname = db)
|
| 120 |
|
|
fs=dbGetQuery(con,paste("SELECT * from swath_geo
|
| 121 |
|
|
WHERE east>=",tile_bb$lon_min," AND
|
| 122 |
|
|
west<=",tile_bb$lon_max," AND
|
| 123 |
|
|
north>=",tile_bb$lat_min," AND
|
| 124 |
|
|
south<=",tile_bb$lat_max," AND
|
| 125 |
|
|
year==",format(as.Date(date,"%Y%m%d"),"%Y")," AND
|
| 126 |
|
|
day==",as.numeric(format(as.Date(date,"%Y%m%d"),"%j"))
|
| 127 |
|
|
))
|
| 128 |
|
|
con=dbDisconnect(con)
|
| 129 |
|
|
fs$id=substr(fs$id,7,19)
|
| 130 |
|
|
## find the swaths on disk (using datadir)
|
| 131 |
|
|
swaths=list.files(datadir,pattern=paste(fs$id,collapse="|"),recursive=T,full=T)
|
| 132 |
|
|
|
| 133 |
|
|
if(verbose) print(paste(nrow(fs)," swath IDs recieved from database and ",length(swaths)," found on disk"))
|
| 134 |
|
|
|
| 135 |
|
|
|
| 136 |
|
|
############################################################################
|
| 137 |
|
|
############################################################################
|
| 138 |
|
|
### Use the HEG tool to grid all available swath data for this date-tile
|
| 139 |
|
|
for(file in swaths){
|
| 140 |
|
|
## Function to generate hegtool parameter file for multi-band HDF-EOS file
|
| 141 |
|
|
print(paste("Starting file",basename(file)))
|
| 142 |
|
|
outfile=paste(tempdir(),"/",basename(file),sep="")
|
| 143 |
|
|
## First write the parameter file (careful, heg is very finicky!)
|
| 144 |
|
|
hdr=paste("NUM_RUNS = ",length(vars$varid),"|MULTI_BAND_HDFEOS:",length(vars$varid),sep="")
|
| 145 |
|
|
grp=paste("
|
| 146 |
|
|
BEGIN
|
| 147 |
|
|
INPUT_FILENAME=",file,"
|
| 148 |
|
|
OBJECT_NAME=mod06
|
| 149 |
|
|
FIELD_NAME=",vars$variable,"|
|
| 150 |
|
|
BAND_NUMBER = 1
|
| 151 |
|
|
OUTPUT_PIXEL_SIZE_X=1000
|
| 152 |
|
|
OUTPUT_PIXEL_SIZE_Y=1000
|
| 153 |
|
|
SPATIAL_SUBSET_UL_CORNER = ( ",upleft," )
|
| 154 |
|
|
SPATIAL_SUBSET_LR_CORNER = ( ",lowright," )
|
| 155 |
|
|
#RESAMPLING_TYPE =",ifelse(grepl("Flag|Mask|Quality",vars),"NN","CUBIC"),"
|
| 156 |
|
|
RESAMPLING_TYPE =NN
|
| 157 |
|
|
OUTPUT_PROJECTION_TYPE = SIN
|
| 158 |
|
|
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 )
|
| 159 |
|
|
# projection parameters from http://landweb.nascom.nasa.gov/cgi-bin/QA_WWW/newPage.cgi?fileName=sn_gctp
|
| 160 |
|
|
ELLIPSOID_CODE = WGS84
|
| 161 |
|
|
OUTPUT_TYPE = HDFEOS
|
| 162 |
|
|
OUTPUT_FILENAME = ",outfile,"
|
| 163 |
|
|
END
|
| 164 |
|
|
|
| 165 |
|
|
",sep="")
|
| 166 |
|
|
|
| 167 |
|
|
## if any remnants from previous runs remain, delete them
|
| 168 |
|
|
if(length(list.files(tempdir(),pattern=basename(file)))>0)
|
| 169 |
|
|
file.remove(list.files(tempdir(),pattern=basename(file),full=T))
|
| 170 |
|
|
## write it to a file
|
| 171 |
|
|
cat(c(hdr,grp) , file=paste(tempdir(),"/",basename(file),"_MODparms.txt",sep=""))
|
| 172 |
|
|
## now run the swath2grid tool
|
| 173 |
|
|
## write the gridded file
|
| 174 |
|
|
log=system(paste("(",swtifpath," -p ",tempdir(),"/",basename(file),"_MODparms.txt -d ; echo $$)",sep=""),intern=T,ignore.stderr=T)
|
| 175 |
|
|
## clean up temporary files in working directory
|
| 176 |
|
|
# file.remove(list.files(pattern=
|
| 177 |
|
|
# paste("filetable.temp_",
|
| 178 |
|
|
# as.numeric(log[length(log)]):(as.numeric(log[length(log)])+3),sep="",collapse="|"))) #Look for files with PID within 3 of parent process
|
| 179 |
|
|
if(verbose) print(log)
|
| 180 |
|
|
print(paste("Finished gridding ", file))
|
| 181 |
|
|
} #end looping over swaths
|
| 182 |
|
|
|
| 183 |
|
|
########################
|
| 184 |
|
|
## confirm at least one file for this date is present. If not, quit.
|
| 185 |
|
|
outfiles=paste(tempdir(),"/",basename(swaths),sep="")
|
| 186 |
|
|
if(!any(file.exists(outfiles))) {
|
| 187 |
|
|
print(paste("######################################## No gridded files for region exist for tile",tile," on date",date))
|
| 188 |
|
|
q("no",status=0)
|
| 189 |
|
|
}
|
| 190 |
|
|
|
| 191 |
|
|
#####################################################
|
| 192 |
|
|
## Process the gridded files to align exactly with MODLAND tile and produce a daily summary of multiple swaths
|
| 193 |
|
|
|
| 194 |
|
|
## Identify output file
|
| 195 |
|
|
ncfile=paste(outdir,"/MOD06_",tile,"_",date,".nc",sep="") #this is the 'final' daily output file
|
| 196 |
|
|
|
| 197 |
|
|
## function to convert binary to decimal to assist in identifying correct values
|
| 198 |
|
|
## this is helpful when defining QA handling below, but isn't used in processing
|
| 199 |
|
|
## b2d=function(x) sum(x * 2^(rev(seq_along(x)) - 1)) #http://tolstoy.newcastle.edu.au/R/e2/help/07/02/10596.html
|
| 200 |
|
|
## for example:
|
| 201 |
|
|
## b2d(c(T,T))
|
| 202 |
|
|
|
| 203 |
|
|
## set Grass to overwrite
|
| 204 |
|
|
Sys.setenv(GRASS_OVERWRITE=1)
|
| 205 |
|
|
Sys.setenv(DEBUG=1)
|
| 206 |
|
|
Sys.setenv(GRASS_GUI="txt")
|
| 207 |
|
|
|
| 208 |
|
|
### Extract various SDSs from a single gridded HDF file and use QA data to throw out 'bad' observations
|
| 209 |
|
|
## make temporary working directory
|
| 210 |
|
|
tf=paste(tempdir(),"/grass", Sys.getpid(),"/", sep="") #temporar
|
| 211 |
|
|
if(!file.exists(tf)) dir.create(tf)
|
| 212 |
|
|
## create output directory if needed
|
| 213 |
|
|
if(!file.exists(dirname(ncfile))) dir.create(dirname(ncfile),recursive=T)
|
| 214 |
|
|
|
| 215 |
|
|
## set up temporary grass instance for this PID
|
| 216 |
|
|
if(verbose) print(paste("Set up temporary grass session in",tf))
|
| 217 |
|
|
initGRASS(gisBase=gisBase,gisDbase=tf,SG=as(td,"SpatialGridDataFrame"),override=T,location="mod06",mapset="PERMANENT",home=tf,pid=Sys.getpid())
|
| 218 |
|
|
system(paste("g.proj -c proj4=\"",projection(td),"\"",sep=""),ignore.stdout=T,ignore.stderr=T)
|
| 219 |
|
|
|
| 220 |
|
|
## Define region by importing one MOD11A1 raster.
|
| 221 |
|
|
print("Import one MOD11A1 raster to define grid")
|
| 222 |
|
|
if(platform=="pleiades") execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",gridfile,"\":MODIS_Grid_Daily_1km_LST:Night_view_angl",sep=""),
|
| 223 |
|
|
output="modisgrid",flags=c("quiet","overwrite","o"))
|
| 224 |
|
|
if(platform=="litoria") 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"))
|
| 225 |
|
|
|
| 226 |
|
|
system("g.region rast=modisgrid save=roi --overwrite",ignore.stdout=F,ignore.stderr=F)
|
| 227 |
|
|
system("g.region rast=modisgrid.1 save=roi --overwrite",ignore.stdout=F,ignore.stderr=F)
|
| 228 |
|
|
|
| 229 |
|
|
## Identify which files to process
|
| 230 |
|
|
tfs=basename(swaths)
|
| 231 |
|
|
## drop swaths that did not produce an output file (typically due to not overlapping the ROI)
|
| 232 |
|
|
tfs=tfs[tfs%in%list.files(tempdir())]
|
| 233 |
|
|
nfs=length(tfs)
|
| 234 |
|
|
if(verbose) print(paste(nfs,"swaths available for processing"))
|
| 235 |
|
|
|
| 236 |
|
|
## loop through scenes and process QA flags
|
| 237 |
|
|
for(i in 1:nfs){
|
| 238 |
|
|
file=paste(tempdir(),"/",tfs[i],sep="")
|
| 239 |
|
|
## Cloud Mask
|
| 240 |
|
|
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod06:Cloud_Mask_1km_0",sep=""),
|
| 241 |
|
|
output=paste("CM1_",i,sep=""),flags=c("overwrite","o")) ; print("")
|
| 242 |
|
|
## extract cloudy and 'probably/confidently clear' pixels
|
| 243 |
|
|
system(paste("r.mapcalc <<EOF
|
| 244 |
|
|
CM_cloud_",i," = ((CM1_",i," / 2^0) % 2) == 1 && ((CM1_",i," / 2^1) % 2^2) == 0
|
| 245 |
|
|
CM_clear_",i," = ((CM1_",i," / 2^0) % 2) == 1 && ((CM1_",i," / 2^1) % 2^2) > 2
|
| 246 |
|
|
CM_path_",i," = ((CM1_",i," / 2^6) % 2^2)
|
| 247 |
|
|
CM_cloud2_",i," = ((CM1_",i," / 2^1) % 2^2)
|
| 248 |
|
|
EOF",sep=""))
|
| 249 |
|
|
## Set CM_cloud2 to null if it is "01" (uncertain)
|
| 250 |
|
|
# execGRASS("r.null",map=paste("CM_cloud2_",i,sep=""),setnull="-9999,1")
|
| 251 |
|
|
|
| 252 |
|
|
## QA
|
| 253 |
|
|
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod06:Quality_Assurance_1km_0",sep=""),
|
| 254 |
|
|
output=paste("QA_",i,sep=""),flags=c("overwrite","o")) ; print("")
|
| 255 |
|
|
## QA_CER
|
| 256 |
|
|
system(paste("r.mapcalc <<EOF
|
| 257 |
|
|
QA_COT_",i,"= ((QA_",i," / 2^0) % 2^1 )==1
|
| 258 |
|
|
QA_COT2_",i,"= ((QA_",i," / 2^1) % 2^2 )>=2
|
| 259 |
|
|
QA_COT3_",i,"= ((QA_",i," / 2^3) % 2^2 )==0
|
| 260 |
|
|
QA_CER_",i,"= ((QA_",i," / 2^5) % 2^1 )==1
|
| 261 |
|
|
QA_CER2_",i,"= ((QA_",i," / 2^6) % 2^2 )>=2
|
| 262 |
|
|
EOF",sep=""))
|
| 263 |
|
|
# QA_CWP_",i,"= ((QA_",i," / 2^8) % 2^1 )==1
|
| 264 |
|
|
# QA_CWP2_",i,"= ((QA_",i," / 2^9) % 2^2 )==3
|
| 265 |
|
|
|
| 266 |
|
|
## Optical Thickness
|
| 267 |
|
|
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod06:Cloud_Optical_Thickness",sep=""),
|
| 268 |
|
|
output=paste("COT_",i,sep=""),
|
| 269 |
|
|
title="cloud_effective_radius",
|
| 270 |
|
|
flags=c("overwrite","o")) ; print("")
|
| 271 |
|
|
execGRASS("r.null",map=paste("COT_",i,sep=""),setnull="-9999")
|
| 272 |
|
|
## keep only positive COT values where quality is 'useful' and '>= good' & scale to real units
|
| 273 |
|
|
system(paste("r.mapcalc \"COT_",i,"=if(QA_COT_",i,"&&QA_COT2_",i,"&&QA_COT3_",i,"&&COT_",i,">=0,COT_",i,",null())\"",sep=""))
|
| 274 |
|
|
## set null COT to 0 in clear-sky pixels
|
| 275 |
|
|
system(paste("r.mapcalc \"COT2_",i,"=if(isnull(COT_",i,")&&CM_clear_",i,"==1,0,COT_",i,")\"",sep=""))
|
| 276 |
|
|
|
| 277 |
|
|
## Effective radius ##
|
| 278 |
|
|
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod06:Cloud_Effective_Radius",sep=""),
|
| 279 |
|
|
output=paste("CER_",i,sep=""),
|
| 280 |
|
|
title="cloud_effective_radius",
|
| 281 |
|
|
flags=c("overwrite","o")) ; print("")
|
| 282 |
|
|
execGRASS("r.null",map=paste("CER_",i,sep=""),setnull="-9999")
|
| 283 |
|
|
## keep only positive CER values where quality is 'useful' and '>= good' & scale to real units
|
| 284 |
|
|
system(paste("r.mapcalc \"CER_",i,"=if(QA_CER_",i,"&&QA_CER2_",i,"&&CER_",i,">=0,CER_",i,",null())\"",sep=""))
|
| 285 |
|
|
## set null CER to 0 in clear-sky pixels
|
| 286 |
|
|
system(paste("r.mapcalc \"CER2_",i,"=if(isnull(CER_",i,")&&CM_clear_",i,"==1,0,CER_",i,")\"",sep=""))
|
| 287 |
|
|
|
| 288 |
|
|
## Cloud Water Path
|
| 289 |
|
|
# execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod06:Cloud_Water_Path",sep=""),
|
| 290 |
|
|
# output=paste("CWP_",i,sep=""),title="cloud_water_path",
|
| 291 |
|
|
# flags=c("overwrite","o")) ; print("")
|
| 292 |
|
|
# execGRASS("r.null",map=paste("CWP_",i,sep=""),setnull="-9999")
|
| 293 |
|
|
# system(paste("r.mapcalc \"CWP_",i,"=if(CWP_",i,"<0,null(),CWP_",i,")\"",sep=""))
|
| 294 |
|
|
## keep only positive CER values where quality is 'useful' and 'very good' & scale to real units
|
| 295 |
|
|
# system(paste("r.mapcalc \"CWP_",i,"=if(QA_CWP_",i,"&&QA_CWP2_",i,"&&CWP_",i,">=0,CWP_",i,",null())\"",sep=""))
|
| 296 |
|
|
## set CER to 0 in clear-sky pixels
|
| 297 |
|
|
# system(paste("r.mapcalc \"CWP2_",i,"=if(CM_clear_",i,"==0,CWP_",i,",0)\"",sep=""))
|
| 298 |
|
|
|
| 299 |
|
|
} #end loop through sub daily files
|
| 300 |
|
|
|
| 301 |
|
|
#### Now generate daily averages (or maximum in case of cloud flag)
|
| 302 |
|
|
|
| 303 |
|
|
system(paste("r.mapcalc <<EOF
|
| 304 |
|
|
COT_numer=",paste("if(isnull(COT_",1:nfs,"),0,COT_",1:nfs,")",sep="",collapse="+"),"
|
| 305 |
|
|
COT_denom=",paste("!isnull(COT_",1:nfs,")",sep="",collapse="+"),"
|
| 306 |
|
|
COT_daily=int(COT_numer/COT_denom)
|
| 307 |
|
|
COT2_numer=",paste("if(isnull(COT2_",1:nfs,"),0,COT2_",1:nfs,")",sep="",collapse="+"),"
|
| 308 |
|
|
COT2_denom=",paste("!isnull(COT2_",1:nfs,")",sep="",collapse="+"),"
|
| 309 |
|
|
COT2_daily=int(COT2_numer/COT2_denom)
|
| 310 |
|
|
CER_numer=",paste("if(isnull(CER_",1:nfs,"),0,CER_",1:nfs,")",sep="",collapse="+"),"
|
| 311 |
|
|
CER_denom=",paste("!isnull(CER_",1:nfs,")",sep="",collapse="+"),"
|
| 312 |
|
|
CER_daily=int(CER_numer/CER_denom)
|
| 313 |
|
|
CER2_numer=",paste("if(isnull(CER2_",1:nfs,"),0,CER2_",1:nfs,")",sep="",collapse="+"),"
|
| 314 |
|
|
CER2_denom=",paste("!isnull(CER2_",1:nfs,")",sep="",collapse="+"),"
|
| 315 |
|
|
CER2_daily=int(CER2_numer/CER2_denom)
|
| 316 |
|
|
CLD_daily=int((max(",paste("if(isnull(CM_cloud_",1:nfs,"),-9999,CM_cloud_",1:nfs,")",sep="",collapse=","),")))
|
| 317 |
|
|
CLD2_daily=int((min(",paste("if(isnull(CM_cloud2_",1:nfs,"),-9999,CM_cloud2_",1:nfs,")",sep="",collapse=","),")))
|
| 318 |
|
|
EOF",sep=""))
|
| 319 |
|
|
|
| 320 |
|
|
execGRASS("r.null",map="CLD_daily",setnull="-9999")
|
| 321 |
|
|
execGRASS("r.null",map="CLD2_daily",setnull="-9999")
|
| 322 |
|
|
|
| 323 |
|
|
|
| 324 |
|
|
### Write the files to a netcdf file
|
| 325 |
|
|
## create image group to facilitate export as multiband netcdf
|
| 326 |
|
|
execGRASS("i.group",group="mod06",input=c("CER_daily","CER2_daily","COT_daily","COT2_daily","CLD_daily","CLD2_daily")) ; print("")
|
| 327 |
|
|
|
| 328 |
|
|
if(file.exists(ncfile)) file.remove(ncfile) #if it exists already, delete it
|
| 329 |
|
|
execGRASS("r.out.gdal",input="mod06",output=ncfile,type="Int16",nodata=-32768,flags=c("quiet"),
|
| 330 |
|
|
# createopt=c("FORMAT=NC4","ZLEVEL=5","COMPRESS=DEFLATE","WRITE_GDAL_TAGS=YES","WRITE_LONLAT=NO"),format="netCDF") #for compressed netcdf
|
| 331 |
|
|
createopt=c("FORMAT=NC","WRITE_GDAL_TAGS=YES","WRITE_LONLAT=NO"),format="netCDF")
|
| 332 |
|
|
|
| 333 |
|
|
system(paste(ncopath,"ncecat -O -u time ",ncfile," ",ncfile,sep=""))
|
| 334 |
|
|
## create temporary nc file with time information to append to MOD06 data
|
| 335 |
|
|
cat(paste("
|
| 336 |
|
|
netcdf time {
|
| 337 |
|
|
dimensions:
|
| 338 |
|
|
time = 1 ;
|
| 339 |
|
|
variables:
|
| 340 |
|
|
int time(time) ;
|
| 341 |
|
|
time:units = \"days since 2000-01-01 00:00:00\" ;
|
| 342 |
|
|
time:calendar = \"gregorian\";
|
| 343 |
|
|
time:long_name = \"time of observation\";
|
| 344 |
|
|
data:
|
| 345 |
|
|
time=",as.integer(as.Date(date,"%Y%m%d")-as.Date("2000-01-01")),";
|
| 346 |
|
|
}"),file=paste(tempdir(),"/time.cdl",sep=""))
|
| 347 |
|
|
system(paste("ncgen -o ",tempdir(),"/time.nc ",tempdir(),"/time.cdl",sep=""))
|
| 348 |
|
|
system(paste(ncopath,"ncks -A ",tempdir(),"/time.nc ",ncfile,sep=""))
|
| 349 |
|
|
## add other attributes
|
| 350 |
|
|
system(paste(ncopath,"ncrename -v Band1,CER -v Band2,CER2 -v Band3,COT -v Band4,COT2 -v Band5,CLD -v Band6,CLD2 ",ncfile,sep=""))
|
| 351 |
|
|
system(paste(ncopath,"ncatted -a scale_factor,CER,o,d,0.01 -a units,CER,o,c,\"micron\" -a missing_value,CER,o,d,-32768 -a long_name,CER,o,c,\"Cloud Particle Effective Radius\" ",ncfile,sep=""))
|
| 352 |
|
|
system(paste(ncopath,"ncatted -a scale_factor,CER2,o,d,0.01 -a units,CER2,o,c,\"micron\" -a missing_value,CER2,o,d,-32768 -a long_name,CER2,o,c,\"Cloud Particle Effective Radius with clear sky set to zero\" ",ncfile,sep=""))
|
| 353 |
|
|
|
| 354 |
|
|
system(paste(ncopath,"ncatted -a scale_factor,COT,o,d,0.01 -a units,COT,o,c,\"none\" -a missing_value,COT,o,d,-32768 -a long_name,COT,o,c,\"Cloud Optical Thickness\" ",ncfile,sep=""))
|
| 355 |
|
|
system(paste(ncopath,"ncatted -a scale_factor,COT2,o,d,0.01 -a units,COT2,o,c,\"none\" -a missing_value,COT2,o,d,-32768 -a long_name,COT2,o,c,\"Cloud Optical Thickness with clear sky set to zero\" ",ncfile,sep=""))
|
| 356 |
|
|
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=""))
|
| 357 |
|
|
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=""))
|
| 358 |
|
|
|
| 359 |
|
|
# system(paste(ncopath,"ncatted -a sourcecode,global,o,c,",script," ",ncfile,sep=""))
|
| 360 |
|
|
|
| 361 |
|
|
|
| 362 |
|
|
### delete the temporary files
|
| 363 |
|
|
unlink_.gislock()
|
| 364 |
|
|
system(paste("rm -frR ",tf,sep=""))
|
| 365 |
|
|
|
| 366 |
|
|
|
| 367 |
|
|
## Confirm that the file has the correct attributes, otherwise delete it
|
| 368 |
|
|
ntime=as.numeric(system(paste("cdo -s ntime ",ncfile),intern=T))
|
| 369 |
|
|
## confirm it has all 'final variables as specified above"
|
| 370 |
|
|
fvar=all(finalvars%in%strsplit(system(paste("cdo -s showvar ",ncfile),intern=T)," ")[[1]])
|
| 371 |
|
|
|
| 372 |
|
|
if(ntime!=1|!fvar) {
|
| 373 |
|
|
print(paste("FILE ERROR: tile ",tile," and date ",date," was not outputted correctly, deleting... "))
|
| 374 |
|
|
file.remove(ncfile)
|
| 375 |
|
|
}
|
| 376 |
|
|
|
| 377 |
|
|
## print out some info
|
| 378 |
|
|
print(paste(" ################################################################### Finished ",date,
|
| 379 |
|
|
"################################################################"))
|
| 380 |
|
|
|
| 381 |
|
|
## delete old files
|
| 382 |
|
|
system("cleartemp")
|
| 383 |
|
|
|
| 384 |
|
|
## quit
|
| 385 |
|
|
q("no",status=0)
|
