| 1 |
1 |
###################################################################################
|
| 2 |
|
### R code to aquire and process MOD06_L2 cloud data from the MODIS platform
|
|
2 |
### R code to aquire and process MOD35_L2 cloud data from the MODIS platform
|
| 3 |
3 |
|
| 4 |
4 |
|
| 5 |
5 |
# Redirect all warnings to stderr()
|
| ... | ... | |
| 17 |
17 |
require(raster)
|
| 18 |
18 |
require(rgdal)
|
| 19 |
19 |
require(spgrass6)
|
| 20 |
|
require(RSQLite)
|
| 21 |
20 |
|
| 22 |
21 |
## get options
|
| 23 |
22 |
opta <- getopt(matrix(c(
|
| ... | ... | |
| 33 |
32 |
q(status=1);
|
| 34 |
33 |
}
|
| 35 |
34 |
|
|
35 |
testing=T
|
|
36 |
platform="pleiades"
|
| 36 |
37 |
|
| 37 |
38 |
## 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
|
|
39 |
if(testing){
|
|
40 |
date="20090129"
|
|
41 |
tile="h11v08"
|
|
42 |
tile="h17v00"
|
|
43 |
verbose=T
|
|
44 |
}
|
| 42 |
45 |
|
| 43 |
46 |
## 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 |
|
|
|
47 |
if(!testing){
|
|
48 |
date=opta$date
|
|
49 |
tile=opta$tile
|
|
50 |
verbose=opta$verbose #print out extensive information for debugging?
|
|
51 |
}
|
| 50 |
52 |
## get year and doy from date
|
| 51 |
53 |
year=format(as.Date(date,"%Y%m%d"),"%Y")
|
| 52 |
54 |
doy=format(as.Date(date,"%Y%m%d"),"%j")
|
| 53 |
55 |
|
| 54 |
56 |
if(platform=="pleiades"){
|
| 55 |
|
## location of MOD06 files
|
| 56 |
|
datadir=paste("/nobackupp4/datapool/modis/MOD06_L2.005/",year,"/",doy,"/",sep="")
|
|
57 |
## location of MOD35 files
|
|
58 |
datadir=paste("/nobackupp4/datapool/modis/MOD09GA.005/",year,"/",doy,"/",sep="")
|
| 57 |
59 |
## path to some executables
|
| 58 |
60 |
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 |
61 |
## specify working directory
|
| 63 |
|
setwd("/nobackupp1/awilso10/mod06")
|
|
62 |
outdir=paste("/nobackupp1/awilso10/mod09/daily/",tile,"/",sep="") #directory for separate daily files
|
|
63 |
basedir="/nobackupp1/awilso10/mod09/" #directory to hold files temporarily before transferring to lou
|
|
64 |
setwd(tempdir())
|
|
65 |
## grass database
|
| 64 |
66 |
gisBase="/u/armichae/pr/grass-6.4.2/"
|
| 65 |
67 |
## 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 |
gridfile=list.files("/nobackupp1/awilso10/mod35/MODTILES/",pattern=tile,full=T)[1]
|
|
69 |
td=raster(gridfile)
|
| 68 |
70 |
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 |
71 |
}
|
| 70 |
72 |
|
| 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 |
73 |
### print some status messages
|
| 89 |
74 |
if(verbose) print(paste("Processing tile",tile," for date",date))
|
| 90 |
75 |
|
| 91 |
76 |
## load tile information and get bounding box
|
| 92 |
|
load(file="modlandTiles.Rdata")
|
|
77 |
load(file="/nobackupp1/awilso10/mod35/modlandTiles.Rdata")
|
| 93 |
78 |
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 |
79 |
|
| 116 |
80 |
#####################################################
|
| 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 |
81 |
|
| 194 |
|
## Identify output file
|
| 195 |
|
ncfile=paste(outdir,"/MOD06_",tile,"_",date,".nc",sep="") #this is the 'final' daily output file
|
| 196 |
|
|
| 197 |
82 |
## function to convert binary to decimal to assist in identifying correct values
|
| 198 |
83 |
## this is helpful when defining QA handling below, but isn't used in processing
|
| 199 |
84 |
## b2d=function(x) sum(x * 2^(rev(seq_along(x)) - 1)) #http://tolstoy.newcastle.edu.au/R/e2/help/07/02/10596.html
|
| ... | ... | |
| 210 |
95 |
tf=paste(tempdir(),"/grass", Sys.getpid(),"/", sep="") #temporar
|
| 211 |
96 |
if(!file.exists(tf)) dir.create(tf)
|
| 212 |
97 |
## create output directory if needed
|
|
98 |
## Identify output file
|
|
99 |
ncfile=paste(outdir,"MOD09cloud_",tile,"_",date,".nc",sep="") #this is the 'final' daily output file
|
| 213 |
100 |
if(!file.exists(dirname(ncfile))) dir.create(dirname(ncfile),recursive=T)
|
| 214 |
|
|
|
101 |
|
| 215 |
102 |
## set up temporary grass instance for this PID
|
| 216 |
103 |
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())
|
|
104 |
initGRASS(gisBase=gisBase,gisDbase=tf,SG=as(td,"SpatialGridDataFrame"),override=T,location="mod09",mapset="PERMANENT",home=tf,pid=Sys.getpid())
|
| 218 |
105 |
system(paste("g.proj -c proj4=\"",projection(td),"\"",sep=""),ignore.stdout=T,ignore.stderr=T)
|
| 219 |
106 |
|
| 220 |
107 |
## Define region by importing one MOD11A1 raster.
|
| 221 |
108 |
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)
|
|
109 |
if(platform=="pleiades") {
|
|
110 |
execGRASS("r.in.gdal",input=td@file@name,output="modisgrid",flags=c("quiet","overwrite","o"))
|
|
111 |
system("g.region rast=modisgrid save=roi --overwrite",ignore.stdout=F,ignore.stderr=F)
|
|
112 |
}
|
|
113 |
|
|
114 |
if(platform=="litoria"){
|
|
115 |
execGRASS("r.in.gdal",input=paste("NETCDF:\"/home/adamw/acrobates/adamw/projects/interp/data/modis/mod06/summary/MOD06_",tile,".nc\":CER",sep=""),
|
|
116 |
output="modisgrid",flags=c("overwrite","o"))
|
|
117 |
system("g.region rast=modisgrid.1 save=roi --overwrite",ignore.stdout=F,ignore.stderr=F)
|
|
118 |
}
|
| 228 |
119 |
|
| 229 |
120 |
## 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())]
|
|
121 |
tfs=unique(sub("CM_|QA_|SenZen_|SolZen_","",basename(outfiles)))
|
|
122 |
#tfs=list.files(tempdir(),pattern="temp.*hdf")
|
| 233 |
123 |
nfs=length(tfs)
|
| 234 |
124 |
if(verbose) print(paste(nfs,"swaths available for processing"))
|
| 235 |
125 |
|
| 236 |
126 |
## loop through scenes and process QA flags
|
| 237 |
127 |
for(i in 1:nfs){
|
| 238 |
|
file=paste(tempdir(),"/",tfs[i],sep="")
|
|
128 |
bfile=tfs[i]
|
|
129 |
## Read in the data from the HDFs
|
| 239 |
130 |
## Cloud Mask
|
| 240 |
|
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod06:Cloud_Mask_1km_0",sep=""),
|
|
131 |
execGRASS("r.in.gdal",input=paste("CM_",bfile,sep=""),
|
| 241 |
132 |
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=""),
|
|
133 |
## QA ## extract first bit to keep only "useful" values of cloud mask
|
|
134 |
execGRASS("r.in.gdal",input=paste("QA_",bfile,sep=""),
|
| 254 |
135 |
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=""))
|
|
136 |
## Sensor Zenith ## extract first bit to keep only "low angle" observations
|
|
137 |
execGRASS("r.in.gdal",input=paste("SenZen_",bfile,sep=""),
|
|
138 |
output=paste("SZ_",i,sep=""),flags=c("overwrite","o")) ; print("")
|
| 276 |
139 |
|
| 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=""))
|
|
140 |
## check for interpolation artefacts
|
|
141 |
# execGRASS("r.stats",input=paste("SZ_",i,sep=""),output=paste("SZ_",i,".txt",sep=""),flags=c("c"))
|
|
142 |
# execGRASS("r.clump",input=paste("SZ_",i,sep=""),output=paste("SZ_",i,"_clump",sep=""))
|
|
143 |
# execGRASS("r.stats",input=paste("SZ_",i,"_clump",sep=""),output="-",flags=c("c"))
|
|
144 |
|
|
145 |
## write out the table of weights for use in the neighborhood analysis to identify bad pixels from swtif
|
|
146 |
p=75 #must be odd
|
|
147 |
mat=matrix(rep(0,p*p),nrow=p)
|
|
148 |
mat[0.5+p/2,]=1
|
|
149 |
cat(mat,file="weights.txt")
|
|
150 |
execGRASS("r.neighbors",input=paste("SZ_",i,sep=""),output=paste("SZ_",i,"clump",sep=""),method="range",size=p,weight="weights.txt") # too slow!
|
|
151 |
system(paste("r.mapcalc \"SZ_",i,"_clump2=SZ_",i,"clump==0\"",sep=""))
|
|
152 |
|
|
153 |
# p=-50:50
|
|
154 |
# system(paste("r.mapcalc \"SZ_",i,"_clump=if(min(",paste("SZ_",i,"[0,",p,"]",sep="",collapse=","),")==max(",paste("SZ_",i,"[0,",p,"]",sep="",collapse=","),"),1,0)\"",sep=""))
|
|
155 |
# system(paste("r.mapcalc \"SZ_",i,"_clump=if(min(",paste("SZ_",i,"[0,",min(p,"]",sep="",collapse=","),")==max(",paste("SZ_",i,"[0,",p,"]",sep="",collapse=","),"),1,0)\"",sep=""))
|
|
156 |
# vals=do.call(rbind.data.frame,strsplit(execGRASS("r.stats",input=paste("SZ_",i,"_clump",sep=""),output="-",flags=c("c"),intern=T),split=" "))
|
|
157 |
# colnames(vals)=c("value","count")
|
|
158 |
# vals$count=as.numeric(as.character(vals$count))
|
|
159 |
# vals$value=as.numeric(as.character(vals$value))
|
|
160 |
# vals=na.omit(vals)
|
|
161 |
# vals$count[vals$value==1&vals$count>10]
|
|
162 |
#
|
|
163 |
#plot(p~value,data=vals)
|
|
164 |
# print(sum(vals$p[vals$p>.1]))
|
| 298 |
165 |
|
| 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=","),")))
|
|
166 |
## Solar Zenith ## extract first bit to keep only "low angle" observations
|
|
167 |
execGRASS("r.in.gdal",input=paste("SolZen_",bfile,sep=""),
|
|
168 |
output=paste("SoZ_",i,sep=""),flags=c("overwrite","o")) ; print("")
|
|
169 |
## produce the summaries
|
|
170 |
system(paste("r.mapcalc <<EOF
|
|
171 |
CM_fill_",i," = if(isnull(CM1_",i,"),1,0)
|
|
172 |
QA_useful_",i," = if((QA_",i," / 2^0) % 2==1,1,0)
|
|
173 |
SZ_low_",i," = if(SZ_",i,"_clump2==0&SZ_",i,"<6000,1,0)
|
|
174 |
SoZ_low_",i," = if(SoZ_",i,"<8500,1,0)
|
|
175 |
CM_dayflag_",i," = if((CM1_",i," / 2^3) % 2==1,1,0)
|
|
176 |
CM_cloud_",i," = if((CM1_",i," / 2^0) % 2==1,(CM1_",i," / 2^1) % 2^2,null())
|
|
177 |
SZday_",i," = if(SZ_",i,"_clump2==0&CM_dayflag_",i,"==1,SZ_",i,",null())
|
|
178 |
SZnight_",i," = if(SZ_",i,"_clump2==0&CM_dayflag_",i,"==0,SZ_",i,",null())
|
|
179 |
CMday_",i," = if(SoZ_low_",i,"==1&SZ_low_",i,"==1&QA_useful_",i,"==1&CM_dayflag_",i,"==1,CM_cloud_",i,",null())
|
|
180 |
CMnight_",i," = if(SZ_low_",i,"==1&QA_useful_",i,"==1&CM_dayflag_",i,"==0,CM_cloud_",i,",null())
|
| 318 |
181 |
EOF",sep=""))
|
| 319 |
182 |
|
| 320 |
|
execGRASS("r.null",map="CLD_daily",setnull="-9999")
|
| 321 |
|
execGRASS("r.null",map="CLD2_daily",setnull="-9999")
|
|
183 |
# CM_dayflag_",i," = if((CM1_",i," / 2^3) % 2==1,1,0)
|
|
184 |
# 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))))
|
|
185 |
# CM_nscore_",i," = if((CM_dayflag_",i,"==1|isnull(CM1_",i,")),0,if(QA_useful_",i,"==0,1,if(SZ_",i,">=6000,2,4)))
|
|
186 |
|
|
187 |
drawplot=F
|
|
188 |
if(drawplot){
|
|
189 |
d2=stack(
|
|
190 |
# raster(readRAST6(paste("QA_useful_",i,sep=""))),
|
|
191 |
raster(readRAST6(paste("CM1_",i,sep=""))),
|
|
192 |
# raster(readRAST6(paste("CM_cloud_",i,sep=""))),
|
|
193 |
# raster(readRAST6(paste("CM_dayflag_",i,sep=""))),
|
|
194 |
# raster(readRAST6(paste("CMday_",i,sep=""))),
|
|
195 |
# raster(readRAST6(paste("CMnight_",i,sep=""))),
|
|
196 |
# raster(readRAST6(paste("CM_fill_",i,sep=""))),
|
|
197 |
# raster(readRAST6(paste("SoZ_",i,sep=""))),
|
|
198 |
raster(readRAST6(paste("SZ_",i,sep=""))),
|
|
199 |
raster(readRAST6(paste("SZ_",i,"_clump",sep=""))),
|
|
200 |
raster(readRAST6(paste("SZ_",i,"_clump2",sep="")))
|
|
201 |
)
|
|
202 |
plot(d2,add=F)
|
|
203 |
}
|
|
204 |
|
|
205 |
|
|
206 |
} #end loop through sub daily files
|
|
207 |
|
|
208 |
## select lowest view angle
|
|
209 |
## use r.series to find minimum
|
|
210 |
system(paste("r.series input=",paste("SZnight_",1:nfs,sep="",collapse=",")," output=SZnight_min method=min_raster",sep=""))
|
|
211 |
system(paste("r.series input=",paste("SZday_",1:nfs,sep="",collapse=",")," output=SZday_min method=min_raster",sep=""))
|
|
212 |
|
|
213 |
## select cloud observation with lowest sensor zenith for day and night
|
|
214 |
system(
|
|
215 |
paste("r.mapcalc <<EOF
|
|
216 |
CMday_daily=",paste(paste("if((SZday_min+1)==",1:nfs,",CMday_",1:nfs,",",sep="",collapse=" "),"null()",paste(rep(")",times=nfs),sep="",collapse="")),"
|
|
217 |
CMnight_daily=",paste(paste("if((SZnight_min+1)==",1:nfs,",CMnight_",1:nfs,",",sep="",collapse=" "),"null()",paste(rep(")",times=nfs),sep="",collapse=""))
|
|
218 |
))
|
|
219 |
|
|
220 |
if(plot){
|
|
221 |
ps=1:nfs
|
|
222 |
ps=c(12,14,17)
|
|
223 |
sz1=brick(lapply(ps,function(i) raster(readRAST6(paste("SZnight_",i,sep="")))))
|
|
224 |
sz_clump=brick(lapply(ps,function(i) raster(readRAST6(paste("SZ_",i,"_clump2",sep="")))))
|
|
225 |
d=brick(lapply(ps,function(i) raster(readRAST6(paste("CMnight_",i,sep="")))))
|
|
226 |
d2=brick(list(raster(readRAST6("SZday_min")),raster(readRAST6("SZnight_min")),raster(readRAST6("CMday_daily")),raster(readRAST6("CMnight_daily"))))
|
|
227 |
library(rasterVis)
|
|
228 |
levelplot(sz1,col.regions=rainbow(100),at=seq(min(sz1@data@min),max(sz1@data@max),len=100))
|
|
229 |
levelplot(sz_clump)
|
|
230 |
levelplot(d)
|
|
231 |
levelplot(d2)
|
|
232 |
}
|
| 322 |
233 |
|
| 323 |
234 |
|
| 324 |
235 |
### Write the files to a netcdf file
|
| 325 |
236 |
## 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"),
|
|
237 |
execGRASS("i.group",group="mod35",input=c("CMday_daily","CMnight_daily")) ; print("")
|
|
238 |
|
|
239 |
if(file.exists(ncfile)) file.remove(ncfile) #if it exists already, delete it
|
|
240 |
execGRASS("r.out.gdal",input="mod35",output=ncfile,type="Byte",nodata=255,flags=c("quiet"),
|
| 330 |
241 |
# createopt=c("FORMAT=NC4","ZLEVEL=5","COMPRESS=DEFLATE","WRITE_GDAL_TAGS=YES","WRITE_LONLAT=NO"),format="netCDF") #for compressed netcdf
|
| 331 |
242 |
createopt=c("FORMAT=NC","WRITE_GDAL_TAGS=YES","WRITE_LONLAT=NO"),format="netCDF")
|
| 332 |
243 |
|
| ... | ... | |
| 347 |
258 |
system(paste("ncgen -o ",tempdir(),"/time.nc ",tempdir(),"/time.cdl",sep=""))
|
| 348 |
259 |
system(paste(ncopath,"ncks -A ",tempdir(),"/time.nc ",ncfile,sep=""))
|
| 349 |
260 |
## 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=""))
|
|
261 |
system(paste(ncopath,"ncrename -v Band1,CMday -v Band2,CMnight ",ncfile,sep=""))
|
|
262 |
system(paste(ncopath,"ncatted ",
|
|
263 |
" -a units,CMday,o,c,\"Cloud Flag (0-3)\" ",
|
|
264 |
" -a missing_value,CMday,o,b,255 ",
|
|
265 |
" -a _FillValue,CMday,o,b,255 ",
|
|
266 |
" -a valid_range,CMday,o,b,\"0,3\" ",
|
|
267 |
" -a long_name,CMday,o,c,\"Cloud Flag from day pixels\" ",
|
|
268 |
" -a units,CMnight,o,c,\"Cloud Flag (0-3)\" ",
|
|
269 |
" -a missing_value,CMnight,o,b,255 ",
|
|
270 |
" -a _FillValue,CMnight,o,b,255 ",
|
|
271 |
" -a valid_range,CMnight,o,b,\"0,3\" ",
|
|
272 |
" -a long_name,CMnight,o,c,\"Cloud Flag from night pixels\" ",
|
|
273 |
ncfile,sep=""))
|
|
274 |
#system(paste(ncopath,"ncatted -a sourcecode,global,o,c,",script," ",ncfile,sep=""))
|
| 360 |
275 |
|
| 361 |
|
|
| 362 |
|
### delete the temporary files
|
| 363 |
|
unlink_.gislock()
|
| 364 |
|
system(paste("rm -frR ",tf,sep=""))
|
| 365 |
|
|
| 366 |
276 |
|
| 367 |
277 |
## Confirm that the file has the correct attributes, otherwise delete it
|
| 368 |
278 |
ntime=as.numeric(system(paste("cdo -s ntime ",ncfile),intern=T))
|
| ... | ... | |
| 373 |
283 |
print(paste("FILE ERROR: tile ",tile," and date ",date," was not outputted correctly, deleting... "))
|
| 374 |
284 |
file.remove(ncfile)
|
| 375 |
285 |
}
|
| 376 |
|
|
|
286 |
############ copy files to lou
|
|
287 |
#if(platform=="pleiades"){
|
|
288 |
# archivedir=paste("MOD35/",outdir,"/",sep="") #directory to create on lou
|
|
289 |
# system(paste("ssh -q bridge2 \"ssh -q lou mkdir -p ",archivedir,"\"",sep=""))
|
|
290 |
# system(paste("ssh -q bridge2 \"scp -q ",ncfile," lou:",archivedir,"\"",sep=""))
|
|
291 |
# file.remove(ncfile)
|
|
292 |
# file.remove(paste(ncfile,".aux.xml",sep=""))
|
|
293 |
#}
|
|
294 |
|
|
295 |
|
|
296 |
### delete the temporary files
|
|
297 |
# unlink_.gislock()
|
|
298 |
# system(paste("rm -frR ",tempdir(),sep=""))
|
|
299 |
|
|
300 |
|
| 377 |
301 |
## print out some info
|
| 378 |
|
print(paste(" ################################################################### Finished ",date,
|
| 379 |
|
"################################################################"))
|
|
302 |
print(paste("####################### Finished ",tile,"-",date, "###################################"))
|
| 380 |
303 |
|
| 381 |
304 |
## delete old files
|
| 382 |
|
system("cleartemp")
|
|
305 |
#system("cleartemp")
|
| 383 |
306 |
|
| 384 |
307 |
## quit
|
| 385 |
308 |
q("no",status=0)
|
updated to use new swtif