Revision 165c04fb
Added by Adam Wilson over 11 years ago
| climate/procedures/MOD06_L2_process.r | ||
|---|---|---|
| 26 | 26 |
} |
| 27 | 27 |
|
| 28 | 28 |
|
| 29 |
date=opta$date #date="20030301"
|
|
| 29 |
date=opta$date #date="20101225"
|
|
| 30 | 30 |
tile=opta$tile #tile="h11v08" |
| 31 | 31 |
verbose=opta$verbose #print out extensive information for debugging? |
| 32 | 32 |
outdir=paste("daily/",tile,"/",sep="") #directory for separate daily files
|
| ... | ... | |
| 42 | 42 |
require(raster) |
| 43 | 43 |
library(rgdal) |
| 44 | 44 |
require(spgrass6) |
| 45 |
require(RSQLite) |
|
| 45 |
#require(RSQLite)
|
|
| 46 | 46 |
|
| 47 | 47 |
|
| 48 | 48 |
## specify some working directories |
| ... | ... | |
| 115 | 115 |
## write the gridded file and save the log including the pid of the parent process |
| 116 | 116 |
# log=system(paste("( /nobackupp1/awilso10/software/heg/bin/swtif -p ",tempdir(),"/",basename(file),"_MODparms.txt -d ; echo $$)",sep=""),intern=T)
|
| 117 | 117 |
log=system(paste("( /nobackupp1/awilso10/software/heg/bin/swtif -p ",tempdir(),"/",basename(file),"_MODparms.txt -d ; echo $$)",sep=""),intern=T,ignore.stderr=T)
|
| 118 |
log=system(paste("(sudo MRTDATADIR=/usr/local/heg/data PGSHOME=/usr/local/heg/TOOLKIT_MTD PWD=/home/adamw /usr/local/heg/bin/swtif -p ",tempdir(),"/",basename(file),"_MODparms.txt -d )",sep=""))
|
|
| 119 |
|
|
| 118 | 120 |
## clean up temporary files in working directory |
| 119 | 121 |
file.remove(list.files(pattern= |
| 120 | 122 |
paste("filetable.temp_",
|
| ... | ... | |
| 147 | 149 |
|
| 148 | 150 |
## set up temporary grass instance for this PID |
| 149 | 151 |
initGRASS(gisBase="/u/armichae/pr/grass-6.4.2/",gisDbase=tf,SG=td,override=T,location="mod06",mapset="PERMANENT",home=tf,pid=Sys.getpid()) |
| 152 |
initGRASS(gisBase="/usr/lib/grass64/",SG=td,gisDbase=tf,override=T,location="mod06",mapset="PERMANENT",home=tf,pid=Sys.getpid()) |
|
| 150 | 153 |
system(paste("g.proj -c proj4=\"",projection(td),"\"",sep=""),ignore.stdout=T,ignore.stderr=T)
|
| 151 | 154 |
|
| 152 | 155 |
## Define region by importing one MOD11A1 raster. |
| 153 | 156 |
print("Import one MOD11A1 raster to define grid")
|
| 154 |
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",gridfile,"\":MODIS_Grid_Daily_1km_LST:Night_view_angl",sep=""),
|
|
| 155 |
output="modisgrid",flags=c("quiet","overwrite","o"))
|
|
| 156 |
system("g.region rast=modisgrid save=roi --overwrite",ignore.stdout=T,ignore.stderr=T)
|
|
| 157 |
|
|
| 157 |
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",gridfile,"\":MODIS_Grid_Daily_1km_LST:Night_view_angl",sep=""), output="modisgrid",flags=c("quiet","overwrite","o"))
|
|
| 158 |
execGRASS("r.in.gdal",input=paste("NETCDF:\"",gridfile,"\":CER",sep=""), output="modisgrid",flags=c("o"))
|
|
| 159 |
system("g.region rast=modisgrid.1 save=roi --overwrite",ignore.stdout=T,ignore.stderr=T)
|
|
| 158 | 160 |
## Identify which files to process |
| 159 | 161 |
tfs=fs$file[fs$dateid==date] |
| 160 | 162 |
## drop swaths that did not produce an output file (typically due to not overlapping the ROI) |
| 161 |
tfs=tfs[tfs%in%list.files(tempdir())] |
|
| 163 |
tfs=tfs[tfs%in%list.files(tempdir(),pattern="hdf$")]
|
|
| 162 | 164 |
nfs=length(tfs) |
| 163 | 165 |
|
| 164 | 166 |
## loop through scenes and process QA flags |
| ... | ... | |
| 167 | 169 |
## Cloud Mask |
| 168 | 170 |
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod06:Cloud_Mask_1km_0",sep=""),
|
| 169 | 171 |
output=paste("CM1_",i,sep=""),flags=c("overwrite","o")) ; print("")
|
| 170 |
## extract cloudy and 'confidently clear' pixels |
|
| 171 | 172 |
system(paste("r.mapcalc <<EOF
|
| 172 |
CM_cloud_",i," = ((CM1_",i," / 2^0) % 2) == 1 && ((CM1_",i," / 2^1) % 2^2) == 0 |
|
| 173 |
CM_clear_",i," = ((CM1_",i," / 2^0) % 2) == 1 && ((CM1_",i," / 2^1) % 2^2) == 3 |
|
| 173 |
CM_determined_",i," = ((CM1_",i," / 2^0) % 2) |
|
| 174 |
CM_state_",i," = ((CM1_",i," / 2^1) % 2^2) |
|
| 175 |
EOF",sep="")) |
|
| 176 |
|
|
| 177 |
## extract cloudy and 'confidently clear' pixels |
|
| 178 |
system(paste("r.mapcalc <<EOF
|
|
| 179 |
CM_cloud_",i," = (((CM1_",i," / 2^0) % 2) == 1) && (((CM1_",i," / 2^1) % 2^2) == 0 ) |
|
| 180 |
CM_clear_",i," = ((CM1_",i," / 2^0) % 2) == 1 && ( ((CM1_",i," / 2^1) % 2^2) > 2 ) |
|
| 174 | 181 |
EOF",sep="")) |
| 175 | 182 |
## QA |
| 176 | 183 |
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod06:Quality_Assurance_1km_0",sep=""),
|
| ... | ... | |
| 178 | 185 |
## QA_CER |
| 179 | 186 |
system(paste("r.mapcalc <<EOF
|
| 180 | 187 |
QA_COT_",i,"= ((QA_",i," / 2^0) % 2^1 )==1 |
| 181 |
QA_COT2_",i,"= ((QA_",i," / 2^1) % 2^2 )==3
|
|
| 188 |
QA_COT2_",i,"= ((QA_",i," / 2^1) % 2^2 )>=2
|
|
| 182 | 189 |
QA_COT3_",i,"= ((QA_",i," / 2^3) % 2^2 )==0 |
| 183 | 190 |
QA_CER_",i,"= ((QA_",i," / 2^5) % 2^1 )==1 |
| 184 |
QA_CER2_",i,"= ((QA_",i," / 2^6) % 2^2 )==3
|
|
| 191 |
QA_CER2_",i,"= ((QA_",i," / 2^6) % 2^2 )>=2
|
|
| 185 | 192 |
EOF",sep="")) |
| 186 | 193 |
# QA_CWP_",i,"= ((QA_",i," / 2^8) % 2^1 )==1 |
| 187 | 194 |
# QA_CWP2_",i,"= ((QA_",i," / 2^9) % 2^2 )==3 |
Also available in: Unified diff
Merging with Pleiades copies