1 |
3133718a
|
Adam M. Wilson
|
###################################################################################
|
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 |
848e34fd
|
Adam M. Wilson @ pfe
|
'profile','p',0,'logical',
|
28 |
3133718a
|
Adam M. Wilson
|
'help', 'h', 0, 'logical'
|
29 |
|
|
), ncol=4, byrow=TRUE))
|
30 |
|
|
if ( !is.null(opta$help) )
|
31 |
|
|
{
|
32 |
|
|
prg <- commandArgs()[1];
|
33 |
|
|
cat(paste("Usage: ", prg, " --date | -d <file> :: The date to process\n", sep=""));
|
34 |
|
|
q(status=1);
|
35 |
|
|
}
|
36 |
|
|
|
37 |
e4e30b86
|
Adam M. Wilson @ pfe
|
testing=F
|
38 |
171a4e16
|
Adam M. Wilson @ pfe
|
platform="pleiades"
|
39 |
3133718a
|
Adam M. Wilson
|
|
40 |
848e34fd
|
Adam M. Wilson @ pfe
|
## record profiling information if requested
|
41 |
|
|
if(opta$profile) Rprof("/nobackupp1/awilso10/mod35/log/profile.out")
|
42 |
|
|
|
43 |
3133718a
|
Adam M. Wilson
|
## default date and tile to play with (will be overwritten below when running in batch)
|
44 |
171a4e16
|
Adam M. Wilson @ pfe
|
if(testing){
|
45 |
|
|
date="20090129"
|
46 |
e4e30b86
|
Adam M. Wilson @ pfe
|
tile="h11v08"
|
47 |
171a4e16
|
Adam M. Wilson @ pfe
|
tile="h17v00"
|
48 |
|
|
verbose=T
|
49 |
|
|
}
|
50 |
3133718a
|
Adam M. Wilson
|
|
51 |
|
|
## now update using options if given
|
52 |
171a4e16
|
Adam M. Wilson @ pfe
|
if(!testing){
|
53 |
|
|
date=opta$date
|
54 |
|
|
tile=opta$tile
|
55 |
|
|
verbose=opta$verbose #print out extensive information for debugging?
|
56 |
|
|
}
|
57 |
3133718a
|
Adam M. Wilson
|
## get year and doy from date
|
58 |
|
|
year=format(as.Date(date,"%Y%m%d"),"%Y")
|
59 |
|
|
doy=format(as.Date(date,"%Y%m%d"),"%j")
|
60 |
|
|
|
61 |
|
|
if(platform=="pleiades"){
|
62 |
aba23d60
|
Adam M. Wilson @ pfe
|
## location of MOD35 files
|
63 |
35d59dc1
|
Adam M. Wilson @ pfe
|
datadir=paste("/nobackupp4/datapool/modis/MOD35_L2.006/",year,"/",doy,"/",sep="")
|
64 |
3133718a
|
Adam M. Wilson
|
## path to some executables
|
65 |
|
|
ncopath="/nasa/sles11/nco/4.0.8/gcc/mpt/bin/"
|
66 |
3bb88310
|
Adam M. Wilson @ pfe
|
swtifpath="/nobackupp1/awilso10/software/heg/bin/swtif_2.12"
|
67 |
|
|
# swtifpath="/nobackupp4/pvotava/software/heg/2.12/bin/swtif"
|
68 |
3133718a
|
Adam M. Wilson
|
## path to swath database
|
69 |
|
|
db="/nobackupp4/pvotava/DB/export/swath_geo.sql.sqlite3.db"
|
70 |
|
|
## specify working directory
|
71 |
5af36cdd
|
Adam M. Wilson @ pfe
|
outdir=paste("/nobackupp1/awilso10/mod35/daily/",tile,"/",sep="") #directory for separate daily files
|
72 |
|
|
basedir="/nobackupp1/awilso10/mod35/" #directory to hold files temporarily before transferring to lou
|
73 |
aba23d60
|
Adam M. Wilson @ pfe
|
setwd(tempdir())
|
74 |
|
|
## grass database
|
75 |
3133718a
|
Adam M. Wilson
|
gisBase="/u/armichae/pr/grass-6.4.2/"
|
76 |
|
|
## path to MOD11A1 file for this tile to align grid/extent
|
77 |
e4e30b86
|
Adam M. Wilson @ pfe
|
gridfile=list.files("/nobackupp1/awilso10/mod35/MODTILES/",pattern=tile,full=T)[1]
|
78 |
|
|
td=raster(gridfile)
|
79 |
3133718a
|
Adam M. Wilson
|
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 "
|
80 |
|
|
}
|
81 |
|
|
|
82 |
|
|
if(platform=="litoria"){ #if running on local server, use different paths
|
83 |
|
|
## specify working directory
|
84 |
f27d277d
|
Adam M. Wilson
|
setwd("~/acrobates/adamw/projects/interp")
|
85 |
aba23d60
|
Adam M. Wilson @ pfe
|
outdir=paste("daily/",tile,"/",sep="") #directory for separate daily files
|
86 |
|
|
basedir=outdir
|
87 |
3133718a
|
Adam M. Wilson
|
gisBase="/usr/lib/grass64"
|
88 |
|
|
## location of MOD06 files
|
89 |
f27d277d
|
Adam M. Wilson
|
datadir="~/acrobates/adamw/projects/interp/data/modis/mod35"
|
90 |
3133718a
|
Adam M. Wilson
|
## path to some executables
|
91 |
|
|
ncopath=""
|
92 |
|
|
swtifpath="sudo MRTDATADIR=\"/usr/local/heg/data\" PGSHOME=/usr/local/heg/TOOLKIT_MTD PWD=/home/adamw /usr/local/heg/bin/swtif"
|
93 |
|
|
## path to swath database
|
94 |
f27d277d
|
Adam M. Wilson
|
db="~/acrobates/adamw/projects/interp/data/modis/mod06/swath_geo.sql.sqlite3.db"
|
95 |
3133718a
|
Adam M. Wilson
|
## get grid file
|
96 |
f27d277d
|
Adam M. Wilson
|
td=raster(paste("~/acrobates/adamw/projects/interp/data/modis/mod06/summary/MOD06_",tile,".nc",sep=""),varname="CER")
|
97 |
3133718a
|
Adam M. Wilson
|
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 "
|
98 |
|
|
}
|
99 |
|
|
|
100 |
|
|
|
101 |
|
|
### print some status messages
|
102 |
848e34fd
|
Adam M. Wilson @ pfe
|
if(verbose) writeLines(paste("STATUS: Beginning ",tile,date))
|
103 |
3133718a
|
Adam M. Wilson
|
|
104 |
|
|
## load tile information and get bounding box
|
105 |
aba23d60
|
Adam M. Wilson @ pfe
|
load(file="/nobackupp1/awilso10/mod35/modlandTiles.Rdata")
|
106 |
3133718a
|
Adam M. Wilson
|
tile_bb=tb[tb$tile==tile,] ## identify tile of interest
|
107 |
be64daa8
|
Adam M. Wilson @ pfe
|
|
108 |
|
|
## get bounds of swath to keep and feed into grass when generating tile
|
109 |
|
|
## expand a little (0.5 deg) to ensure that there is no clipping of pixels on the edges
|
110 |
|
|
## tile will later be aligned with MODLAND tile so the extra will eventually be trimmed
|
111 |
5af36cdd
|
Adam M. Wilson @ pfe
|
upleft=paste(min(90,tile_bb$lat_max+0.5),max(-180,tile_bb$lon_min-0.5)) #northwest corner
|
112 |
be64daa8
|
Adam M. Wilson @ pfe
|
lowright=paste(max(-90,tile_bb$lat_min-0.5),min(180,tile_bb$lon_max+0.5)) #southeast corner
|
113 |
3133718a
|
Adam M. Wilson
|
|
114 |
|
|
## vector of variables expected to be in final netcdf file. If these are not present, the file will be deleted at the end.
|
115 |
5af36cdd
|
Adam M. Wilson @ pfe
|
finalvars=c("CMday","CMnight")
|
116 |
3133718a
|
Adam M. Wilson
|
|
117 |
|
|
|
118 |
|
|
#####################################################
|
119 |
|
|
##find swaths in region from sqlite database for the specified date/tile
|
120 |
|
|
if(verbose) print("Accessing swath ID's from database")
|
121 |
|
|
con=dbConnect("SQLite", dbname = db)
|
122 |
f27d277d
|
Adam M. Wilson
|
fs=dbGetQuery(con,paste("SELECT * from swath_geo6
|
123 |
3133718a
|
Adam M. Wilson
|
WHERE east>=",tile_bb$lon_min," AND
|
124 |
|
|
west<=",tile_bb$lon_max," AND
|
125 |
|
|
north>=",tile_bb$lat_min," AND
|
126 |
|
|
south<=",tile_bb$lat_max," AND
|
127 |
|
|
year==",format(as.Date(date,"%Y%m%d"),"%Y")," AND
|
128 |
|
|
day==",as.numeric(format(as.Date(date,"%Y%m%d"),"%j"))
|
129 |
|
|
))
|
130 |
|
|
con=dbDisconnect(con)
|
131 |
|
|
fs$id=substr(fs$id,7,19)
|
132 |
|
|
## find the swaths on disk (using datadir)
|
133 |
|
|
swaths=list.files(datadir,pattern=paste(fs$id,collapse="|"),recursive=T,full=T)
|
134 |
|
|
|
135 |
848e34fd
|
Adam M. Wilson @ pfe
|
### print some status messages
|
136 |
|
|
if(verbose) writeLines(paste("STATUS:swaths tile:",tile,"date:",date,"swathIDs:",nrow(fs)," swathsOnDisk:",length(swaths)))
|
137 |
3133718a
|
Adam M. Wilson
|
|
138 |
|
|
|
139 |
171a4e16
|
Adam M. Wilson @ pfe
|
## define function that grids swaths
|
140 |
|
|
swtif<-function(file,var){
|
141 |
|
|
outfile=paste(tempdir(),"/",var$varid,"_",basename(file),sep="") #gridded path
|
142 |
|
|
## First write the parameter file (careful, heg is very finicky!)
|
143 |
|
|
hdr=paste("NUM_RUNS = 1")
|
144 |
|
|
grp=paste("
|
145 |
3133718a
|
Adam M. Wilson
|
BEGIN
|
146 |
|
|
INPUT_FILENAME=",file,"
|
147 |
|
|
OBJECT_NAME=mod35
|
148 |
171a4e16
|
Adam M. Wilson @ pfe
|
FIELD_NAME=",var$variable,"|
|
149 |
|
|
BAND_NUMBER = ",var$band,"
|
150 |
5af36cdd
|
Adam M. Wilson @ pfe
|
OUTPUT_PIXEL_SIZE_X=926.6
|
151 |
|
|
OUTPUT_PIXEL_SIZE_Y=926.6
|
152 |
|
|
# MODIS 1km Resolution
|
153 |
3133718a
|
Adam M. Wilson
|
SPATIAL_SUBSET_UL_CORNER = ( ",upleft," )
|
154 |
|
|
SPATIAL_SUBSET_LR_CORNER = ( ",lowright," )
|
155 |
171a4e16
|
Adam M. Wilson @ pfe
|
RESAMPLING_TYPE =",var$method,"
|
156 |
3133718a
|
Adam M. Wilson
|
OUTPUT_PROJECTION_TYPE = SIN
|
157 |
|
|
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 )
|
158 |
|
|
# projection parameters from http://landweb.nascom.nasa.gov/cgi-bin/QA_WWW/newPage.cgi?fileName=sn_gctp
|
159 |
|
|
ELLIPSOID_CODE = WGS84
|
160 |
|
|
OUTPUT_TYPE = HDFEOS
|
161 |
171a4e16
|
Adam M. Wilson @ pfe
|
OUTPUT_FILENAME = ",outfile,"
|
162 |
3133718a
|
Adam M. Wilson
|
END
|
163 |
|
|
",sep="")
|
164 |
|
|
## write it to a file
|
165 |
|
|
cat(c(hdr,grp) , file=paste(tempdir(),"/",basename(file),"_MODparms.txt",sep=""))
|
166 |
|
|
## now run the swath2grid tool
|
167 |
|
|
## write the gridded file
|
168 |
e0f4e7ec
|
Adam M. Wilson @ pfe
|
system(paste(swtifpath," -p ",tempdir(),"/",basename(file),"_MODparms.txt -d -tmpLatLondir ",tempdir(),sep=""),intern=F,ignore.stderr=F)
|
169 |
171a4e16
|
Adam M. Wilson @ pfe
|
print(paste("Finished processing variable",var$variable,"from ",basename(file),"to",outfile))
|
170 |
|
|
}
|
171 |
a981dca4
|
Adam M. Wilson @ pfe
|
|
172 |
171a4e16
|
Adam M. Wilson @ pfe
|
## vars to grid
|
173 |
|
|
vars=as.data.frame(matrix(c(
|
174 |
|
|
"Cloud_Mask", "CM", "NN", 1,
|
175 |
|
|
"Quality_Assurance", "QA", "NN", 1,
|
176 |
|
|
"Solar_Zenith", "SolZen", "NN", 1,
|
177 |
|
|
"Sensor_Zenith", "SenZen", "CUBIC", 1
|
178 |
|
|
),
|
179 |
|
|
byrow=T,ncol=4,dimnames=list(1:4,c("variable","varid","method","band"))),stringsAsFactors=F)
|
180 |
a981dca4
|
Adam M. Wilson @ pfe
|
|
181 |
|
|
|
182 |
171a4e16
|
Adam M. Wilson @ pfe
|
############################################################################
|
183 |
|
|
############################################################################
|
184 |
|
|
### Use the HEG tool to grid all available swath data for this date-tile
|
185 |
|
|
for(file in swaths){
|
186 |
|
|
print(paste("Starting file",basename(file)))
|
187 |
|
|
## run swtif for each band
|
188 |
|
|
lapply(1:nrow(vars),function(i) swtif(file,vars[i,]))
|
189 |
3133718a
|
Adam M. Wilson
|
} #end looping over swaths
|
190 |
|
|
|
191 |
171a4e16
|
Adam M. Wilson @ pfe
|
|
192 |
|
|
#############################################################################
|
193 |
848e34fd
|
Adam M. Wilson @ pfe
|
## check for zero dimension in HDFs
|
194 |
|
|
## occasionlly swtif will output a hdf with a resolution of 0. Not sure why, but drop them here.
|
195 |
|
|
CMcheck=list.files(pattern="CM_.*hdf$") #list of files to check
|
196 |
|
|
CM_0=do.call(c,lapply(CMcheck, function(f) any(res(raster(f))==0)))
|
197 |
|
|
keep=sub("CM_","",CMcheck[!CM_0])
|
198 |
|
|
if(length(keep)<length(CMcheck)){writeLines(paste("Warning (Resolution of zero): ",paste(sub("CM_","",CMcheck)[!sub("CM_","",CMcheck)%in%keep],collapse=",")," from ",tile," for ",date))}
|
199 |
|
|
outfiles=list.files(tempdir(),full=T,pattern=paste(keep,"$",sep="",collapse="|"))
|
200 |
|
|
if(length(outfiles)==0) {
|
201 |
3133718a
|
Adam M. Wilson
|
print(paste("######################################## No gridded files for region exist for tile",tile," on date",date))
|
202 |
|
|
q("no",status=0)
|
203 |
|
|
}
|
204 |
|
|
|
205 |
848e34fd
|
Adam M. Wilson @ pfe
|
## confirm at least one file for this date is present. If not, quit.
|
206 |
|
|
#outfiles=list.files(tempdir(),full=T,pattern=paste(basename(swaths),"$",sep="",collapse="|"))
|
207 |
|
|
#if(!any(file.exists(outfiles))) {
|
208 |
|
|
# print(paste("######################################## No gridded files for region exist for tile",tile," on date",date))
|
209 |
|
|
# q("no",status=0)
|
210 |
|
|
#}
|
211 |
|
|
|
212 |
|
|
|
213 |
e0f4e7ec
|
Adam M. Wilson @ pfe
|
plot=F
|
214 |
|
|
if(plot){
|
215 |
|
|
i=1
|
216 |
95354b03
|
Adam M. Wilson @ pfe
|
system(paste("gdalinfo ",outfiles[19]))
|
217 |
|
|
d=lapply(outfiles,function(r) raster(r))
|
218 |
e5c2e69b
|
Adam M. Wilson @ pfe
|
summary(d[[6]])
|
219 |
e0f4e7ec
|
Adam M. Wilson @ pfe
|
}
|
220 |
a981dca4
|
Adam M. Wilson @ pfe
|
#system(paste("scp ",outfiles[1]," adamw@acrobates.eeb.yale.edu:/data/personal/adamw/projects/interp/tmp/",sep=""))
|
221 |
e0f4e7ec
|
Adam M. Wilson @ pfe
|
|
222 |
3133718a
|
Adam M. Wilson
|
#####################################################
|
223 |
|
|
## Process the gridded files to align exactly with MODLAND tile and produce a daily summary of multiple swaths
|
224 |
|
|
|
225 |
|
|
## function to convert binary to decimal to assist in identifying correct values
|
226 |
|
|
## this is helpful when defining QA handling below, but isn't used in processing
|
227 |
|
|
## b2d=function(x) sum(x * 2^(rev(seq_along(x)) - 1)) #http://tolstoy.newcastle.edu.au/R/e2/help/07/02/10596.html
|
228 |
|
|
## for example:
|
229 |
|
|
## b2d(c(T,T))
|
230 |
|
|
|
231 |
|
|
## set Grass to overwrite
|
232 |
|
|
Sys.setenv(GRASS_OVERWRITE=1)
|
233 |
|
|
Sys.setenv(DEBUG=1)
|
234 |
|
|
Sys.setenv(GRASS_GUI="txt")
|
235 |
|
|
|
236 |
|
|
### Extract various SDSs from a single gridded HDF file and use QA data to throw out 'bad' observations
|
237 |
|
|
## make temporary working directory
|
238 |
|
|
tf=paste(tempdir(),"/grass", Sys.getpid(),"/", sep="") #temporar
|
239 |
|
|
if(!file.exists(tf)) dir.create(tf)
|
240 |
|
|
## create output directory if needed
|
241 |
aba23d60
|
Adam M. Wilson @ pfe
|
## Identify output file
|
242 |
5af36cdd
|
Adam M. Wilson @ pfe
|
ncfile=paste(outdir,"MOD35_",tile,"_",date,".nc",sep="") #this is the 'final' daily output file
|
243 |
|
|
if(!file.exists(dirname(ncfile))) dir.create(dirname(ncfile),recursive=T)
|
244 |
aba23d60
|
Adam M. Wilson @ pfe
|
|
245 |
3133718a
|
Adam M. Wilson
|
## set up temporary grass instance for this PID
|
246 |
|
|
if(verbose) print(paste("Set up temporary grass session in",tf))
|
247 |
a981dca4
|
Adam M. Wilson @ pfe
|
initGRASS(gisBase=gisBase,gisDbase=tf,SG=as(td,"SpatialGridDataFrame"),override=T,location="mod35",mapset="PERMANENT",home=tf,pid=Sys.getpid())
|
248 |
3133718a
|
Adam M. Wilson
|
system(paste("g.proj -c proj4=\"",projection(td),"\"",sep=""),ignore.stdout=T,ignore.stderr=T)
|
249 |
|
|
|
250 |
|
|
## Define region by importing one MOD11A1 raster.
|
251 |
|
|
print("Import one MOD11A1 raster to define grid")
|
252 |
f27d277d
|
Adam M. Wilson
|
if(platform=="pleiades") {
|
253 |
e4e30b86
|
Adam M. Wilson @ pfe
|
execGRASS("r.in.gdal",input=td@file@name,output="modisgrid",flags=c("quiet","overwrite","o"))
|
254 |
f27d277d
|
Adam M. Wilson
|
system("g.region rast=modisgrid save=roi --overwrite",ignore.stdout=F,ignore.stderr=F)
|
255 |
|
|
}
|
256 |
|
|
|
257 |
|
|
if(platform=="litoria"){
|
258 |
|
|
execGRASS("r.in.gdal",input=paste("NETCDF:\"/home/adamw/acrobates/adamw/projects/interp/data/modis/mod06/summary/MOD06_",tile,".nc\":CER",sep=""),
|
259 |
|
|
output="modisgrid",flags=c("overwrite","o"))
|
260 |
|
|
system("g.region rast=modisgrid.1 save=roi --overwrite",ignore.stdout=F,ignore.stderr=F)
|
261 |
|
|
}
|
262 |
a981dca4
|
Adam M. Wilson @ pfe
|
|
263 |
3133718a
|
Adam M. Wilson
|
## Identify which files to process
|
264 |
171a4e16
|
Adam M. Wilson @ pfe
|
tfs=unique(sub("CM_|QA_|SenZen_|SolZen_","",basename(outfiles)))
|
265 |
3133718a
|
Adam M. Wilson
|
#tfs=list.files(tempdir(),pattern="temp.*hdf")
|
266 |
|
|
nfs=length(tfs)
|
267 |
|
|
if(verbose) print(paste(nfs,"swaths available for processing"))
|
268 |
|
|
|
269 |
|
|
## loop through scenes and process QA flags
|
270 |
|
|
for(i in 1:nfs){
|
271 |
848e34fd
|
Adam M. Wilson @ pfe
|
bfile=tfs[i]
|
272 |
171a4e16
|
Adam M. Wilson @ pfe
|
## Read in the data from the HDFs
|
273 |
3133718a
|
Adam M. Wilson
|
## Cloud Mask
|
274 |
848e34fd
|
Adam M. Wilson @ pfe
|
GDALinfo(paste("CM_",bfile,sep=""),returnStats=F,silent=T)
|
275 |
171a4e16
|
Adam M. Wilson @ pfe
|
execGRASS("r.in.gdal",input=paste("CM_",bfile,sep=""),
|
276 |
3133718a
|
Adam M. Wilson
|
output=paste("CM1_",i,sep=""),flags=c("overwrite","o")) ; print("")
|
277 |
171a4e16
|
Adam M. Wilson @ pfe
|
## QA ## extract first bit to keep only "useful" values of cloud mask
|
278 |
|
|
execGRASS("r.in.gdal",input=paste("QA_",bfile,sep=""),
|
279 |
f27d277d
|
Adam M. Wilson
|
output=paste("QA_",i,sep=""),flags=c("overwrite","o")) ; print("")
|
280 |
171a4e16
|
Adam M. Wilson @ pfe
|
## Sensor Zenith ## extract first bit to keep only "low angle" observations
|
281 |
|
|
execGRASS("r.in.gdal",input=paste("SenZen_",bfile,sep=""),
|
282 |
a981dca4
|
Adam M. Wilson @ pfe
|
output=paste("SZ_",i,sep=""),flags=c("overwrite","o")) ; print("")
|
283 |
171a4e16
|
Adam M. Wilson @ pfe
|
## Solar Zenith ## extract first bit to keep only "low angle" observations
|
284 |
|
|
execGRASS("r.in.gdal",input=paste("SolZen_",bfile,sep=""),
|
285 |
a981dca4
|
Adam M. Wilson @ pfe
|
output=paste("SoZ_",i,sep=""),flags=c("overwrite","o")) ; print("")
|
286 |
171a4e16
|
Adam M. Wilson @ pfe
|
## produce the summaries
|
287 |
a981dca4
|
Adam M. Wilson @ pfe
|
system(paste("r.mapcalc <<EOF
|
288 |
|
|
CM_fill_",i," = if(isnull(CM1_",i,"),1,0)
|
289 |
f27d277d
|
Adam M. Wilson
|
QA_useful_",i," = if((QA_",i," / 2^0) % 2==1,1,0)
|
290 |
95354b03
|
Adam M. Wilson @ pfe
|
SZ_low_",i," = if(SZ_",i,"<6000,1,0)
|
291 |
a981dca4
|
Adam M. Wilson @ pfe
|
SoZ_low_",i," = if(SoZ_",i,"<8500,1,0)
|
292 |
5af36cdd
|
Adam M. Wilson @ pfe
|
CM_dayflag_",i," = if((CM1_",i," / 2^3) % 2==1,1,0)
|
293 |
171a4e16
|
Adam M. Wilson @ pfe
|
CM_cloud_",i," = if((CM1_",i," / 2^0) % 2==1,(CM1_",i," / 2^1) % 2^2,null())
|
294 |
95354b03
|
Adam M. Wilson @ pfe
|
SZday_",i," = if(CM_dayflag_",i,"==1,SZ_",i,",null())
|
295 |
|
|
SZnight_",i," = if(CM_dayflag_",i,"==0,SZ_",i,",null())
|
296 |
171a4e16
|
Adam M. Wilson @ pfe
|
CMday_",i," = if(SoZ_low_",i,"==1&SZ_low_",i,"==1&QA_useful_",i,"==1&CM_dayflag_",i,"==1,CM_cloud_",i,",null())
|
297 |
|
|
CMnight_",i," = if(SZ_low_",i,"==1&QA_useful_",i,"==1&CM_dayflag_",i,"==0,CM_cloud_",i,",null())
|
298 |
3133718a
|
Adam M. Wilson
|
EOF",sep=""))
|
299 |
|
|
|
300 |
a981dca4
|
Adam M. Wilson @ pfe
|
# CM_dayflag_",i," = if((CM1_",i," / 2^3) % 2==1,1,0)
|
301 |
171a4e16
|
Adam M. Wilson @ pfe
|
# CM_dscore_",i," = if((CM_dayflag_",i,"==0|isnull(CM1_",i,")),0,if(QA_useful_",i,"==0,1,if(SZ_",i,">=6000,2,if(SoZ_",i,">=8500,3,4))))
|
302 |
|
|
# CM_nscore_",i," = if((CM_dayflag_",i,"==1|isnull(CM1_",i,")),0,if(QA_useful_",i,"==0,1,if(SZ_",i,">=6000,2,4)))
|
303 |
a981dca4
|
Adam M. Wilson @ pfe
|
|
304 |
|
|
drawplot=F
|
305 |
|
|
if(drawplot){
|
306 |
|
|
d2=stack(
|
307 |
171a4e16
|
Adam M. Wilson @ pfe
|
# raster(readRAST6(paste("QA_useful_",i,sep=""))),
|
308 |
e4e30b86
|
Adam M. Wilson @ pfe
|
raster(readRAST6(paste("CM1_",i,sep=""))),
|
309 |
95354b03
|
Adam M. Wilson @ pfe
|
raster(readRAST6(paste("CM_cloud_",i,sep=""))),
|
310 |
|
|
raster(readRAST6(paste("CM_dayflag_",i,sep=""))),
|
311 |
|
|
raster(readRAST6(paste("CMday_",i,sep=""))),
|
312 |
|
|
raster(readRAST6(paste("CMnight_",i,sep=""))),
|
313 |
171a4e16
|
Adam M. Wilson @ pfe
|
# raster(readRAST6(paste("CM_fill_",i,sep=""))),
|
314 |
|
|
# raster(readRAST6(paste("SoZ_",i,sep=""))),
|
315 |
95354b03
|
Adam M. Wilson @ pfe
|
raster(readRAST6(paste("SZ_",i,sep="")))
|
316 |
a981dca4
|
Adam M. Wilson @ pfe
|
)
|
317 |
171a4e16
|
Adam M. Wilson @ pfe
|
plot(d2,add=F)
|
318 |
a981dca4
|
Adam M. Wilson @ pfe
|
}
|
319 |
|
|
|
320 |
|
|
|
321 |
3133718a
|
Adam M. Wilson
|
} #end loop through sub daily files
|
322 |
|
|
|
323 |
a981dca4
|
Adam M. Wilson @ pfe
|
## select lowest view angle
|
324 |
|
|
## use r.series to find minimum
|
325 |
171a4e16
|
Adam M. Wilson @ pfe
|
system(paste("r.series input=",paste("SZnight_",1:nfs,sep="",collapse=",")," output=SZnight_min method=min_raster",sep=""))
|
326 |
|
|
system(paste("r.series input=",paste("SZday_",1:nfs,sep="",collapse=",")," output=SZday_min method=min_raster",sep=""))
|
327 |
a981dca4
|
Adam M. Wilson @ pfe
|
## select cloud observation with lowest sensor zenith for day and night
|
328 |
|
|
system(
|
329 |
|
|
paste("r.mapcalc <<EOF
|
330 |
171a4e16
|
Adam M. Wilson @ pfe
|
CMday_daily=",paste(paste("if((SZday_min+1)==",1:nfs,",CMday_",1:nfs,",",sep="",collapse=" "),"null()",paste(rep(")",times=nfs),sep="",collapse="")),"
|
331 |
|
|
CMnight_daily=",paste(paste("if((SZnight_min+1)==",1:nfs,",CMnight_",1:nfs,",",sep="",collapse=" "),"null()",paste(rep(")",times=nfs),sep="",collapse=""))
|
332 |
a981dca4
|
Adam M. Wilson @ pfe
|
))
|
333 |
|
|
|
334 |
e5c2e69b
|
Adam M. Wilson @ pfe
|
execGRASS("r.null",map="CMday_daily",setnull="255") ; print("")
|
335 |
|
|
execGRASS("r.null",map="CMnight_daily",setnull="255") ; print("")
|
336 |
|
|
|
337 |
171a4e16
|
Adam M. Wilson @ pfe
|
if(plot){
|
338 |
e4e30b86
|
Adam M. Wilson @ pfe
|
ps=1:nfs
|
339 |
e5c2e69b
|
Adam M. Wilson @ pfe
|
ps=c(10,11,13,14)
|
340 |
95354b03
|
Adam M. Wilson @ pfe
|
sz1=brick(lapply(ps,function(i) raster(readRAST6(paste("SZday_",i,sep="")))))
|
341 |
|
|
d=brick(lapply(ps,function(i) raster(readRAST6(paste("CMday_",i,sep="")))))
|
342 |
171a4e16
|
Adam M. Wilson @ pfe
|
d2=brick(list(raster(readRAST6("SZday_min")),raster(readRAST6("SZnight_min")),raster(readRAST6("CMday_daily")),raster(readRAST6("CMnight_daily"))))
|
343 |
|
|
library(rasterVis)
|
344 |
|
|
levelplot(sz1,col.regions=rainbow(100),at=seq(min(sz1@data@min),max(sz1@data@max),len=100))
|
345 |
|
|
levelplot(d)
|
346 |
|
|
levelplot(d2)
|
347 |
|
|
}
|
348 |
f27d277d
|
Adam M. Wilson
|
|
349 |
3133718a
|
Adam M. Wilson
|
|
350 |
|
|
### Write the files to a netcdf file
|
351 |
|
|
## create image group to facilitate export as multiband netcdf
|
352 |
e5c2e69b
|
Adam M. Wilson @ pfe
|
execGRASS("i.group",group="mod35",input=c("CMday_daily","CMnight_daily"),flags=c("quiet")) ; print("")
|
353 |
e4e30b86
|
Adam M. Wilson @ pfe
|
|
354 |
|
|
if(file.exists(ncfile)) file.remove(ncfile) #if it exists already, delete it
|
355 |
e5c2e69b
|
Adam M. Wilson @ pfe
|
execGRASS("r.out.gdal",input="mod35",output=ncfile,type="Byte",nodata=255,flags=c("verbose"),
|
356 |
3133718a
|
Adam M. Wilson
|
# createopt=c("FORMAT=NC4","ZLEVEL=5","COMPRESS=DEFLATE","WRITE_GDAL_TAGS=YES","WRITE_LONLAT=NO"),format="netCDF") #for compressed netcdf
|
357 |
|
|
createopt=c("FORMAT=NC","WRITE_GDAL_TAGS=YES","WRITE_LONLAT=NO"),format="netCDF")
|
358 |
|
|
|
359 |
|
|
system(paste(ncopath,"ncecat -O -u time ",ncfile," ",ncfile,sep=""))
|
360 |
|
|
## create temporary nc file with time information to append to MOD06 data
|
361 |
|
|
cat(paste("
|
362 |
|
|
netcdf time {
|
363 |
|
|
dimensions:
|
364 |
|
|
time = 1 ;
|
365 |
|
|
variables:
|
366 |
|
|
int time(time) ;
|
367 |
|
|
time:units = \"days since 2000-01-01 00:00:00\" ;
|
368 |
|
|
time:calendar = \"gregorian\";
|
369 |
|
|
time:long_name = \"time of observation\";
|
370 |
|
|
data:
|
371 |
|
|
time=",as.integer(as.Date(date,"%Y%m%d")-as.Date("2000-01-01")),";
|
372 |
|
|
}"),file=paste(tempdir(),"/time.cdl",sep=""))
|
373 |
|
|
system(paste("ncgen -o ",tempdir(),"/time.nc ",tempdir(),"/time.cdl",sep=""))
|
374 |
|
|
system(paste(ncopath,"ncks -A ",tempdir(),"/time.nc ",ncfile,sep=""))
|
375 |
|
|
## add other attributes
|
376 |
e5c2e69b
|
Adam M. Wilson @ pfe
|
## need to delete _FillValue becuase r.out.gdal incorrectly calls zero values missing if there are no other missing values in the raster.
|
377 |
|
|
## so need to delete then re-add. If you just change the value, ncatted will change the values in the raster in addition to the attribute.
|
378 |
5af36cdd
|
Adam M. Wilson @ pfe
|
system(paste(ncopath,"ncrename -v Band1,CMday -v Band2,CMnight ",ncfile,sep=""))
|
379 |
c24e32a8
|
Adam M. Wilson @ pfe
|
system(paste(ncopath,"ncatted ",
|
380 |
5af36cdd
|
Adam M. Wilson @ pfe
|
" -a units,CMday,o,c,\"Cloud Flag (0-3)\" ",
|
381 |
|
|
" -a missing_value,CMday,o,b,255 ",
|
382 |
e5c2e69b
|
Adam M. Wilson @ pfe
|
" -a _FillValue,CMday,d,, ",
|
383 |
5af36cdd
|
Adam M. Wilson @ pfe
|
" -a valid_range,CMday,o,b,\"0,3\" ",
|
384 |
e4e30b86
|
Adam M. Wilson @ pfe
|
" -a long_name,CMday,o,c,\"Cloud Flag from day pixels\" ",
|
385 |
5af36cdd
|
Adam M. Wilson @ pfe
|
" -a units,CMnight,o,c,\"Cloud Flag (0-3)\" ",
|
386 |
|
|
" -a missing_value,CMnight,o,b,255 ",
|
387 |
e5c2e69b
|
Adam M. Wilson @ pfe
|
" -a _FillValue,CMnight,d,, ",
|
388 |
5af36cdd
|
Adam M. Wilson @ pfe
|
" -a valid_range,CMnight,o,b,\"0,3\" ",
|
389 |
e4e30b86
|
Adam M. Wilson @ pfe
|
" -a long_name,CMnight,o,c,\"Cloud Flag from night pixels\" ",
|
390 |
c24e32a8
|
Adam M. Wilson @ pfe
|
ncfile,sep=""))
|
391 |
e5c2e69b
|
Adam M. Wilson @ pfe
|
## add the fillvalue attribute back (without changing the actual values)
|
392 |
|
|
system(paste(ncopath,"ncatted -a _FillValue,CMday,o,b,255 ",ncfile,sep=""))
|
393 |
|
|
system(paste(ncopath,"ncatted -a _FillValue,CMnight,o,b,255 ",ncfile,sep=""))
|
394 |
3133718a
|
Adam M. Wilson
|
|
395 |
|
|
|
396 |
|
|
## Confirm that the file has the correct attributes, otherwise delete it
|
397 |
|
|
ntime=as.numeric(system(paste("cdo -s ntime ",ncfile),intern=T))
|
398 |
|
|
## confirm it has all 'final variables as specified above"
|
399 |
|
|
fvar=all(finalvars%in%strsplit(system(paste("cdo -s showvar ",ncfile),intern=T)," ")[[1]])
|
400 |
|
|
|
401 |
|
|
if(ntime!=1|!fvar) {
|
402 |
|
|
print(paste("FILE ERROR: tile ",tile," and date ",date," was not outputted correctly, deleting... "))
|
403 |
|
|
file.remove(ncfile)
|
404 |
|
|
}
|
405 |
aba23d60
|
Adam M. Wilson @ pfe
|
############ copy files to lou
|
406 |
b3344197
|
Adam M. Wilson @ pfe
|
#if(platform=="pleiades"){
|
407 |
|
|
# archivedir=paste("MOD35/",outdir,"/",sep="") #directory to create on lou
|
408 |
|
|
# system(paste("ssh -q bridge2 \"ssh -q lou mkdir -p ",archivedir,"\"",sep=""))
|
409 |
|
|
# system(paste("ssh -q bridge2 \"scp -q ",ncfile," lou:",archivedir,"\"",sep=""))
|
410 |
|
|
# file.remove(ncfile)
|
411 |
|
|
# file.remove(paste(ncfile,".aux.xml",sep=""))
|
412 |
|
|
#}
|
413 |
aba23d60
|
Adam M. Wilson @ pfe
|
|
414 |
|
|
|
415 |
|
|
### delete the temporary files
|
416 |
b3344197
|
Adam M. Wilson @ pfe
|
# unlink_.gislock()
|
417 |
|
|
# system(paste("rm -frR ",tempdir(),sep=""))
|
418 |
aba23d60
|
Adam M. Wilson @ pfe
|
|
419 |
|
|
|
420 |
848e34fd
|
Adam M. Wilson @ pfe
|
### print some status messages
|
421 |
|
|
if(verbose) writeLines(paste("STATUS:end tile:",tile,"date:",date,"swathIDs:",nrow(fs)," swathsOnDisk:",length(swaths),"fileExists:",file.exists(ncfile)))
|
422 |
|
|
|
423 |
|
|
## turn off the profiler
|
424 |
|
|
if(opta$profile) Rprof(NULL)
|
425 |
3133718a
|
Adam M. Wilson
|
|
426 |
|
|
|
427 |
|
|
## quit
|
428 |
|
|
q("no",status=0)
|