Revision 807fa48c
Added by Adam Wilson about 12 years ago
| climate/procedures/MOD06_L2_data_compile_Pleiades.r | ||
|---|---|---|
| 121 | 121 |
|
| 122 | 122 |
### Function to extract various SDSs from a single gridded HDF file and use QA data to throw out 'bad' observations |
| 123 | 123 |
loadcloud<-function(date,fs){
|
| 124 |
### set up unique grass session |
|
| 125 | 124 |
tf=paste(tempdir(),"/grass", Sys.getpid(),"/", sep="") |
| 125 |
print(paste("Set up temporary grass session in",tf))
|
|
| 126 | 126 |
|
| 127 | 127 |
## set up temporary grass instance for this PID |
| 128 | 128 |
initGRASS(gisBase="/nobackupp1/awilso10/software/grass-6.4.3svn",gisDbase=tf,SG=td,override=T,location="mod06",mapset="PERMANENT",home=tf,pid=Sys.getpid()) |
| 129 | 129 |
system(paste("g.proj -c proj4=\"",projection(td),"\"",sep=""))
|
| 130 | 130 |
|
| 131 | 131 |
## Define region by importing one MOD11A1 raster. |
| 132 |
print("Import one MOD11A1 raster to define grid")
|
|
| 132 | 133 |
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",gridfile,"\":MODIS_Grid_Daily_1km_LST:Night_view_angl",sep=""),
|
| 133 | 134 |
output="modisgrid",flags=c("quiet","overwrite","o"))
|
| 134 | 135 |
system("g.region rast=modisgrid save=roi --overwrite")
|
| climate/procedures/MOD06_L2_data_compile_run.r | ||
|---|---|---|
| 1 |
|
|
| 2 |
setwd("/nobackupp1/awilso10/mod06")
|
|
| 3 |
library(multicore) |
|
| 4 |
### get list of files to process |
|
| 5 |
datadir="/nobackupp4/datapool/modis/MOD06_L2.005/" |
|
| 6 |
|
|
| 7 |
fs=data.frame( |
|
| 8 |
path=list.files(datadir,full=T,recursive=T,pattern="hdf"), |
|
| 9 |
file=basename(list.files(datadir,full=F,recursive=T,pattern="hdf"))) |
|
| 10 |
fs$date=as.Date(substr(fs$file,11,17),"%Y%j") |
|
| 11 |
fs$month=format(fs$date,"%m") |
|
| 12 |
fs$year=format(fs$date,"%Y") |
|
| 13 |
fs$time=substr(fs$file,19,22) |
|
| 14 |
fs$datetime=as.POSIXct(strptime(paste(substr(fs$file,11,17),substr(fs$file,19,22)), '%Y%j %H%M')) |
|
| 15 |
fs$dateid=format(fs$date,"%Y%m%d") |
|
| 16 |
fs$path=as.character(fs$path) |
|
| 17 |
fs$file=as.character(fs$file) |
|
| 18 |
|
|
| 19 |
## get all unique dates |
|
| 20 |
alldates=unique(fs$dateid) |
|
| 21 |
|
|
| 22 |
## load tile information |
|
| 23 |
load(file="modlandTiles.Rdata") |
|
| 24 |
### use MODIS tile as ROI |
|
| 25 |
#modt=readOGR("modgrid","modis_sinusoidal_grid_world",)
|
|
| 26 |
#modt@data[,colnames(tb)[3:6]]=tb[match(paste(modt$h,modt$v),paste(tb$ih,tb$iv)),3:6] |
|
| 27 |
#write.csv(modt@data,file="modistile.csv") |
|
| 28 |
|
|
| 29 |
|
|
| 30 |
## write it out |
|
| 31 |
save(fs,tb,file="allfiles.Rdata") |
|
| 32 |
#save(alldates,file="alldates.Rdata") |
|
| 33 |
|
|
| 34 |
## identify which have been completed |
|
| 35 |
outdir="2_daily" |
|
| 36 |
done=alldates%in%substr(list.files(outdir),5,12) |
|
| 37 |
table(done) |
|
| 38 |
notdone=alldates[!done] |
|
| 39 |
|
|
| 40 |
#notdone=alldates[1:4] |
|
| 41 |
|
|
| 42 |
save(notdone,file="notdone.Rdata") |
|
| 43 |
|
|
| 44 |
#write.table(paste("i=",notdone[1:10],sep=""),file="notdone.txt",row.names=F)
|
|
| 45 |
|
|
| 46 |
## vars |
|
| 47 |
vars=as.data.frame(matrix(c( |
|
| 48 |
"Cloud_Effective_Radius", "CER", |
|
| 49 |
"Cloud_Effective_Radius_Uncertainty", "CERU", |
|
| 50 |
"Cloud_Optical_Thickness", "COT", |
|
| 51 |
"Cloud_Optical_Thickness_Uncertainty", "COTU", |
|
| 52 |
"Cloud_Water_Path", "CWP", |
|
| 53 |
"Cloud_Water_Path_Uncertainty", "CWPU", |
|
| 54 |
"Cloud_Phase_Optical_Properties", "CPOP", |
|
| 55 |
"Cloud_Multi_Layer_Flag", "CMLF", |
|
| 56 |
"Cloud_Mask_1km", "CM1", |
|
| 57 |
"Quality_Assurance_1km", "QA"), |
|
| 58 |
byrow=T,ncol=2,dimnames=list(1:10,c("variable","varid"))),stringsAsFactors=F)
|
|
| 59 |
save(vars,file="vars.Rdata") |
|
| 60 |
|
|
| 61 |
library(multicore) |
|
| 62 |
mclapply(1:length(notdone),function(i) system(paste("Rscript --verbose --vanilla /u/awilso10/environmental-layers/climate/procedures/MOD06_L2_data_compile_Pleiades.r i=",i,sep="")))
|
|
| 63 |
|
|
| 64 |
|
|
| 65 |
## finish up and quit R |
|
| 66 |
q("no")
|
|
| climate/procedures/MOD06_L2_process.r | ||
|---|---|---|
| 1 |
################################################################################### |
|
| 2 |
### R code to aquire and process MOD06_L2 cloud data from the MODIS platform |
|
| 3 |
|
|
| 4 |
|
|
| 5 |
## load libraries |
|
| 6 |
require(sp) |
|
| 7 |
require(rgdal) |
|
| 8 |
require(reshape) |
|
| 9 |
require(ncdf4) |
|
| 10 |
require(geosphere) |
|
| 11 |
require(raster) |
|
| 12 |
require(spgrass6) |
|
| 13 |
library(multicore) |
|
| 14 |
|
|
| 15 |
## specify some working directories |
|
| 16 |
setwd("/nobackupp1/awilso10/mod06")
|
|
| 17 |
|
|
| 18 |
outdir="2_daily" |
|
| 19 |
|
|
| 20 |
print(paste("tempdir()=",tempdir()))
|
|
| 21 |
print(paste("TMPDIR=",Sys.getenv("TMPDIR")))
|
|
| 22 |
|
|
| 23 |
### get list of files to process |
|
| 24 |
datadir="/nobackupp4/datapool/modis/MOD06_L2.005/" |
|
| 25 |
|
|
| 26 |
fs=data.frame( |
|
| 27 |
path=list.files(datadir,full=T,recursive=T,pattern="hdf"), |
|
| 28 |
file=basename(list.files(datadir,full=F,recursive=T,pattern="hdf"))) |
|
| 29 |
fs$date=as.Date(substr(fs$file,11,17),"%Y%j") |
|
| 30 |
fs$month=format(fs$date,"%m") |
|
| 31 |
fs$year=format(fs$date,"%Y") |
|
| 32 |
fs$time=substr(fs$file,19,22) |
|
| 33 |
fs$datetime=as.POSIXct(strptime(paste(substr(fs$file,11,17),substr(fs$file,19,22)), '%Y%j %H%M')) |
|
| 34 |
fs$dateid=format(fs$date,"%Y%m%d") |
|
| 35 |
fs$path=as.character(fs$path) |
|
| 36 |
fs$file=as.character(fs$file) |
|
| 37 |
|
|
| 38 |
## get all unique dates |
|
| 39 |
alldates=unique(fs$dateid) |
|
| 40 |
|
|
| 41 |
## load tile information |
|
| 42 |
load(file="modlandTiles.Rdata") |
|
| 43 |
### use MODIS tile as ROI |
|
| 44 |
#modt=readOGR("modgrid","modis_sinusoidal_grid_world",)
|
|
| 45 |
#modt@data[,colnames(tb)[3:6]]=tb[match(paste(modt$h,modt$v),paste(tb$ih,tb$iv)),3:6] |
|
| 46 |
#write.csv(modt@data,file="modistile.csv") |
|
| 47 |
tile="h11v08" #can move this to submit script if needed |
|
| 48 |
|
|
| 49 |
|
|
| 50 |
## identify which have been completed |
|
| 51 |
done=alldates%in%substr(list.files(outdir),5,12) |
|
| 52 |
table(done) |
|
| 53 |
notdone=alldates[!done] #these are the dates that still need to be processed |
|
| 54 |
|
|
| 55 |
## vars to process |
|
| 56 |
vars=as.data.frame(matrix(c( |
|
| 57 |
"Cloud_Effective_Radius", "CER", |
|
| 58 |
"Cloud_Effective_Radius_Uncertainty", "CERU", |
|
| 59 |
"Cloud_Optical_Thickness", "COT", |
|
| 60 |
"Cloud_Optical_Thickness_Uncertainty", "COTU", |
|
| 61 |
"Cloud_Water_Path", "CWP", |
|
| 62 |
"Cloud_Water_Path_Uncertainty", "CWPU", |
|
| 63 |
"Cloud_Phase_Optical_Properties", "CPOP", |
|
| 64 |
"Cloud_Multi_Layer_Flag", "CMLF", |
|
| 65 |
"Cloud_Mask_1km", "CM1", |
|
| 66 |
"Quality_Assurance_1km", "QA"), |
|
| 67 |
byrow=T,ncol=2,dimnames=list(1:10,c("variable","varid"))),stringsAsFactors=F)
|
|
| 68 |
|
|
| 69 |
|
|
| 70 |
############################################################################ |
|
| 71 |
############################################################################ |
|
| 72 |
### Define functions to process a particular date-tile |
|
| 73 |
|
|
| 74 |
swath2grid=function(file,vars,upleft,lowright){
|
|
| 75 |
## Function to generate hegtool parameter file for multi-band HDF-EOS file |
|
| 76 |
print(paste("Starting file",basename(file)))
|
|
| 77 |
outfile=paste(tempdir(),"/",basename(file),sep="") |
|
| 78 |
## First write the parameter file (careful, heg is very finicky!) |
|
| 79 |
hdr=paste("NUM_RUNS = ",length(vars$varid),"|MULTI_BAND_HDFEOS:",length(vars$varid),sep="")
|
|
| 80 |
grp=paste("
|
|
| 81 |
BEGIN |
|
| 82 |
INPUT_FILENAME=",file," |
|
| 83 |
OBJECT_NAME=mod06 |
|
| 84 |
FIELD_NAME=",vars$variable,"| |
|
| 85 |
BAND_NUMBER = 1 |
|
| 86 |
OUTPUT_PIXEL_SIZE_X=1000 |
|
| 87 |
OUTPUT_PIXEL_SIZE_Y=1000 |
|
| 88 |
SPATIAL_SUBSET_UL_CORNER = ( ",upleft," ) |
|
| 89 |
SPATIAL_SUBSET_LR_CORNER = ( ",lowright," ) |
|
| 90 |
#RESAMPLING_TYPE =",ifelse(grepl("Flag|Mask|Quality",vars),"NN","CUBIC"),"
|
|
| 91 |
RESAMPLING_TYPE =NN |
|
| 92 |
OUTPUT_PROJECTION_TYPE = SIN |
|
| 93 |
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 ) |
|
| 94 |
# projection parameters from http://landweb.nascom.nasa.gov/cgi-bin/QA_WWW/newPage.cgi?fileName=sn_gctp |
|
| 95 |
ELLIPSOID_CODE = WGS84 |
|
| 96 |
OUTPUT_TYPE = HDFEOS |
|
| 97 |
OUTPUT_FILENAME = ",outfile," |
|
| 98 |
END |
|
| 99 |
|
|
| 100 |
",sep="") |
|
| 101 |
|
|
| 102 |
## if any remnants from previous runs remain, delete them |
|
| 103 |
if(length(list.files(tempdir(),pattern=basename(file)))>0) |
|
| 104 |
file.remove(list.files(tempdir(),pattern=basename(file),full=T)) |
|
| 105 |
## write it to a file |
|
| 106 |
cat(c(hdr,grp) , file=paste(tempdir(),"/",basename(file),"_MODparms.txt",sep="")) |
|
| 107 |
## now run the swath2grid tool |
|
| 108 |
## write the tiff file |
|
| 109 |
log=system(paste("/nobackupp4/pvotava/software/heg/bin/swtif -p ",paste(tempdir(),"/",basename(file),"_MODparms.txt -d",sep=""),sep=""),intern=T)
|
|
| 110 |
## clean up temporary files in working directory |
|
| 111 |
# file.remove(paste("filetable.temp_",pid,sep=""))
|
|
| 112 |
print(log) |
|
| 113 |
## confirm file is present |
|
| 114 |
print(paste("Confirming output file (",outfile,") is present and readable by GDAL"))
|
|
| 115 |
gi=GDALinfo(outfile); print(gi) |
|
| 116 |
print(paste("Finished ", file))
|
|
| 117 |
} |
|
| 118 |
|
|
| 119 |
|
|
| 120 |
############################################################## |
|
| 121 |
### Import to GRASS for processing |
|
| 122 |
|
|
| 123 |
## function to convert binary to decimal to assist in identifying correct values |
|
| 124 |
## this is helpful when defining QA handling below, but isn't used in processing |
|
| 125 |
## b2d=function(x) sum(x * 2^(rev(seq_along(x)) - 1)) #http://tolstoy.newcastle.edu.au/R/e2/help/07/02/10596.html |
|
| 126 |
## for example: |
|
| 127 |
## b2d(c(T,T)) |
|
| 128 |
|
|
| 129 |
## set Grass to overwrite |
|
| 130 |
Sys.setenv(GRASS_OVERWRITE=1) |
|
| 131 |
Sys.setenv(DEBUG=1) |
|
| 132 |
Sys.setenv(GRASS_GUI="txt") |
|
| 133 |
|
|
| 134 |
### Function to extract various SDSs from a single gridded HDF file and use QA data to throw out 'bad' observations |
|
| 135 |
loadcloud<-function(date,fs){
|
|
| 136 |
tf=paste(tempdir(),"/grass", Sys.getpid(),"/", sep="") |
|
| 137 |
print(paste("Set up temporary grass session in",tf))
|
|
| 138 |
|
|
| 139 |
## load a MOD11A1 file to define grid |
|
| 140 |
gridfile=list.files("/nobackupp4/datapool/modis/MOD11A1.005/2006.01.27/",pattern=paste(tile,".*hdf$",sep=""),full=T)[1]
|
|
| 141 |
td=readGDAL(paste("HDF4_EOS:EOS_GRID:\"",gridfile,"\":MODIS_Grid_Daily_1km_LST:Night_view_angl",sep=""))
|
|
| 142 |
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 " |
|
| 143 |
|
|
| 144 |
## set up temporary grass instance for this PID |
|
| 145 |
initGRASS(gisBase="/nobackupp1/awilso10/software/grass-6.4.3svn",gisDbase=tf,SG=td,override=T,location="mod06",mapset="PERMANENT",home=tf,pid=Sys.getpid()) |
|
| 146 |
system(paste("g.proj -c proj4=\"",projection(td),"\"",sep=""))
|
|
| 147 |
|
|
| 148 |
## Define region by importing one MOD11A1 raster. |
|
| 149 |
print("Import one MOD11A1 raster to define grid")
|
|
| 150 |
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",gridfile,"\":MODIS_Grid_Daily_1km_LST:Night_view_angl",sep=""),
|
|
| 151 |
output="modisgrid",flags=c("quiet","overwrite","o"))
|
|
| 152 |
system("g.region rast=modisgrid save=roi --overwrite")
|
|
| 153 |
|
|
| 154 |
## Identify which files to process |
|
| 155 |
tfs=fs$file[fs$dateid==date] |
|
| 156 |
nfs=length(tfs) |
|
| 157 |
|
|
| 158 |
## loop through scenes and process QA flags |
|
| 159 |
for(i in 1:nfs){
|
|
| 160 |
file=paste(tempdir(),"/",tfs[i],sep="") |
|
| 161 |
## Cloud Mask |
|
| 162 |
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod06:Cloud_Mask_1km_0",sep=""),
|
|
| 163 |
output=paste("CM1_",i,sep=""),flags=c("overwrite","o")) ; print("")
|
|
| 164 |
## extract cloudy and 'confidently clear' pixels |
|
| 165 |
system(paste("r.mapcalc <<EOF
|
|
| 166 |
CM_cloud_",i," = ((CM1_",i," / 2^0) % 2) == 1 && ((CM1_",i," / 2^1) % 2^2) == 0 |
|
| 167 |
CM_clear_",i," = ((CM1_",i," / 2^0) % 2) == 1 && ((CM1_",i," / 2^1) % 2^2) == 3 |
|
| 168 |
EOF",sep="")) |
|
| 169 |
|
|
| 170 |
## QA |
|
| 171 |
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod06:Quality_Assurance_1km_0",sep=""),
|
|
| 172 |
output=paste("QA_",i,sep=""),flags=c("overwrite","o")) ; print("")
|
|
| 173 |
## QA_CER |
|
| 174 |
system(paste("r.mapcalc <<EOF
|
|
| 175 |
QA_COT_",i,"= ((QA_",i," / 2^0) % 2^1 )==1 |
|
| 176 |
QA_COT2_",i,"= ((QA_",i," / 2^1) % 2^2 )==3 |
|
| 177 |
QA_COT3_",i,"= ((QA_",i," / 2^3) % 2^2 )==0 |
|
| 178 |
QA_CER_",i,"= ((QA_",i," / 2^5) % 2^1 )==1 |
|
| 179 |
QA_CER2_",i,"= ((QA_",i," / 2^6) % 2^2 )==3 |
|
| 180 |
EOF",sep="")) |
|
| 181 |
# QA_CWP_",i,"= ((QA_",i," / 2^8) % 2^1 )==1 |
|
| 182 |
# QA_CWP2_",i,"= ((QA_",i," / 2^9) % 2^2 )==3 |
|
| 183 |
|
|
| 184 |
## Optical Thickness |
|
| 185 |
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod06:Cloud_Optical_Thickness",sep=""),
|
|
| 186 |
output=paste("COT_",i,sep=""),
|
|
| 187 |
title="cloud_effective_radius", |
|
| 188 |
flags=c("overwrite","o")) ; print("")
|
|
| 189 |
execGRASS("r.null",map=paste("COT_",i,sep=""),setnull="-9999")
|
|
| 190 |
## keep only positive COT values where quality is 'useful' and 'very good' & scale to real units |
|
| 191 |
system(paste("r.mapcalc \"COT_",i,"=if(QA_COT_",i,"&&QA_COT2_",i,"&&QA_COT3_",i,"&&COT_",i,">=0,COT_",i,"*0.009999999776482582,null())\"",sep=""))
|
|
| 192 |
## set COT to 0 in clear-sky pixels |
|
| 193 |
system(paste("r.mapcalc \"COT2_",i,"=if(CM_clear_",i,"==0,COT_",i,",0)\"",sep=""))
|
|
| 194 |
|
|
| 195 |
## Effective radius ## |
|
| 196 |
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod06:Cloud_Effective_Radius",sep=""),
|
|
| 197 |
output=paste("CER_",i,sep=""),
|
|
| 198 |
title="cloud_effective_radius", |
|
| 199 |
flags=c("overwrite","o")) ; print("")
|
|
| 200 |
execGRASS("r.null",map=paste("CER_",i,sep=""),setnull="-9999")
|
|
| 201 |
## keep only positive CER values where quality is 'useful' and 'very good' & scale to real units |
|
| 202 |
system(paste("r.mapcalc \"CER_",i,"=if(QA_CER_",i,"&&QA_CER2_",i,"&&CER_",i,">=0,CER_",i,"*0.009999999776482582,null())\"",sep=""))
|
|
| 203 |
## set CER to 0 in clear-sky pixels |
|
| 204 |
system(paste("r.mapcalc \"CER2_",i,"=if(CM_clear_",i,"==0,CER_",i,",0)\"",sep=""))
|
|
| 205 |
|
|
| 206 |
## Cloud Water Path |
|
| 207 |
# execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod06:Cloud_Water_Path",sep=""),
|
|
| 208 |
# output=paste("CWP_",i,sep=""),title="cloud_water_path",
|
|
| 209 |
# flags=c("overwrite","o")) ; print("")
|
|
| 210 |
# execGRASS("r.null",map=paste("CWP_",i,sep=""),setnull="-9999")
|
|
| 211 |
## keep only positive CER values where quality is 'useful' and 'very good' & scale to real units |
|
| 212 |
# system(paste("r.mapcalc \"CWP_",i,"=if(QA_CWP_",i,"&&QA_CWP2_",i,"&&CWP_",i,">=0,CWP_",i,"*0.009999999776482582,null())\"",sep=""))
|
|
| 213 |
## set CER to 0 in clear-sky pixels |
|
| 214 |
# system(paste("r.mapcalc \"CWP2_",i,"=if(CM_clear_",i,"==0,CWP_",i,",0)\"",sep=""))
|
|
| 215 |
|
|
| 216 |
} #end loop through sub daily files |
|
| 217 |
|
|
| 218 |
#### Now generate daily averages (or maximum in case of cloud flag) |
|
| 219 |
|
|
| 220 |
system(paste("r.mapcalc <<EOF
|
|
| 221 |
COT_denom=",paste("!isnull(COT2_",1:nfs,")",sep="",collapse="+"),"
|
|
| 222 |
COT_numer=",paste("if(isnull(COT2_",1:nfs,"),0,COT2_",1:nfs,")",sep="",collapse="+"),"
|
|
| 223 |
COT_daily=COT_numer/COT_denom |
|
| 224 |
CER_denom=",paste("!isnull(CER2_",1:nfs,")",sep="",collapse="+"),"
|
|
| 225 |
CER_numer=",paste("if(isnull(CER2_",1:nfs,"),0,CER2_",1:nfs,")",sep="",collapse="+"),"
|
|
| 226 |
CER_daily=CER_numer/CER_denom |
|
| 227 |
CLD_daily=max(",paste("if(isnull(CM_cloud_",1:nfs,"),0,CM_cloud_",1:nfs,")",sep="",collapse=","),")
|
|
| 228 |
EOF",sep="")) |
|
| 229 |
|
|
| 230 |
#### Write the files to a geotiff |
|
| 231 |
execGRASS("r.out.gdal",input="CER_daily",output=paste(outdir,"/CER_",date,".tif",sep=""),nodata=-999,flags=c("quiet"))
|
|
| 232 |
execGRASS("r.out.gdal",input="COT_daily",output=paste(outdir,"/COT_",date,".tif",sep=""),nodata=-999,flags=c("quiet"))
|
|
| 233 |
execGRASS("r.out.gdal",input="CLD_daily",output=paste(outdir,"/CLD_",date,".tif",sep=""),nodata=99,flags=c("quiet"))
|
|
| 234 |
|
|
| 235 |
### delete the temporary files |
|
| 236 |
unlink_.gislock() |
|
| 237 |
system("/nobackupp1/awilso10/software/grass-6.4.3svn/etc/clean_temp")
|
|
| 238 |
system(paste("rm -R ",tf,sep=""))
|
|
| 239 |
} |
|
| 240 |
|
|
| 241 |
|
|
| 242 |
########################################### |
|
| 243 |
### Define a wrapper function that will call the two functions above (gridding and QA-handling) for a single tile-date |
|
| 244 |
|
|
| 245 |
mod06<-function(date,tile){
|
|
| 246 |
print(paste("Processing date ",date," for tile",tile))
|
|
| 247 |
##################################################### |
|
| 248 |
## Run the gridding procedure |
|
| 249 |
tile_bb=tb[tb$tile==tile,] ## identify tile of interest |
|
| 250 |
lapply(fs$path[fs$dateid==date],swath2grid,vars=vars,upleft=paste(tile_bb$lat_max,tile_bb$lon_min),lowright=paste(tile_bb$lat_min,tile_bb$lon_max)) |
|
| 251 |
|
|
| 252 |
##################################################### |
|
| 253 |
## Process the gridded files |
|
| 254 |
|
|
| 255 |
## temporary objects to test function below |
|
| 256 |
# i=1 |
|
| 257 |
#file=paste(outdir,"/",fs$file[1],sep="") |
|
| 258 |
#date=as.Date("2000-05-23")
|
|
| 259 |
|
|
| 260 |
## run themod06 processing for this date |
|
| 261 |
loadcloud(date,fs=fs) |
|
| 262 |
## print out some info |
|
| 263 |
print(paste(" ################################################################### Finished ",date,"
|
|
| 264 |
################################################################")) |
|
| 265 |
} |
|
| 266 |
|
|
| 267 |
## test it |
|
| 268 |
##date=notdone[1] |
|
| 269 |
##mod06(date,tile) |
|
| 270 |
|
|
| 271 |
## run it for all dates |
|
| 272 |
mclapply(notdone[1:100],mod06,tile) |
|
| 273 |
|
|
| 274 |
|
|
| 275 |
## quit R |
|
| 276 |
q("no")
|
|
| 277 |
|
|
| 278 |
|
|
| climate/procedures/Pleiades.R | ||
|---|---|---|
| 1 | 1 |
#### Script to facilitate processing of MOD06 data |
| 2 | 2 |
|
| 3 |
|
|
| 4 | 3 |
setwd("/nobackupp1/awilso10/mod06")
|
| 5 | 4 |
|
| 6 |
### get list of files to process |
|
| 7 |
datadir="/nobackupp4/datapool/modis/MOD06_L2.005/" |
|
| 8 |
|
|
| 9 |
fs=data.frame( |
|
| 10 |
path=list.files(datadir,full=T,recursive=T,pattern="hdf"), |
|
| 11 |
file=basename(list.files(datadir,full=F,recursive=T,pattern="hdf"))) |
|
| 12 |
fs$date=as.Date(substr(fs$file,11,17),"%Y%j") |
|
| 13 |
fs$month=format(fs$date,"%m") |
|
| 14 |
fs$year=format(fs$date,"%Y") |
|
| 15 |
fs$time=substr(fs$file,19,22) |
|
| 16 |
fs$datetime=as.POSIXct(strptime(paste(substr(fs$file,11,17),substr(fs$file,19,22)), '%Y%j %H%M')) |
|
| 17 |
fs$dateid=format(fs$date,"%Y%m%d") |
|
| 18 |
fs$path=as.character(fs$path) |
|
| 19 |
fs$file=as.character(fs$file) |
|
| 20 |
|
|
| 21 |
## get all unique dates |
|
| 22 |
alldates=unique(fs$dateid) |
|
| 23 |
|
|
| 24 |
|
|
| 25 | 5 |
## get MODLAND tile information |
| 26 | 6 |
tb=read.table("http://landweb.nascom.nasa.gov/developers/sn_tiles/sn_bound_10deg.txt",skip=6,nrows=648,header=T)
|
| 27 | 7 |
tb$tile=paste("h",sprintf("%02d",tb$ih),"v",sprintf("%02d",tb$iv),sep="")
|
| 28 |
### use MODIS tile as ROI |
|
| 29 |
#modt=readOGR("modgrid","modis_sinusoidal_grid_world",)
|
|
| 30 |
#modt@data[,colnames(tb)[3:6]]=tb[match(paste(modt$h,modt$v),paste(tb$ih,tb$iv)),3:6] |
|
| 31 |
#write.csv(modt@data,file="modistile.csv") |
|
| 32 |
|
|
| 33 |
|
|
| 34 |
## write it out |
|
| 35 |
save(fs,tb,file="allfiles.Rdata") |
|
| 36 |
#save(alldates,file="alldates.Rdata") |
|
| 37 |
|
|
| 38 |
## identify which have been completed |
|
| 39 |
outdir="2_daily" |
|
| 40 |
done=alldates%in%substr(list.files(outdir),5,12) |
|
| 41 |
table(done) |
|
| 42 |
notdone=alldates[!done] |
|
| 43 |
|
|
| 44 |
#notdone=alldates[1:4] |
|
| 45 |
|
|
| 46 |
save(notdone,file="notdone.Rdata") |
|
| 47 |
|
|
| 48 |
|
|
| 49 |
## vars |
|
| 50 |
vars=as.data.frame(matrix(c( |
|
| 51 |
"Cloud_Effective_Radius", "CER", |
|
| 52 |
"Cloud_Effective_Radius_Uncertainty", "CERU", |
|
| 53 |
"Cloud_Optical_Thickness", "COT", |
|
| 54 |
"Cloud_Optical_Thickness_Uncertainty", "COTU", |
|
| 55 |
"Cloud_Water_Path", "CWP", |
|
| 56 |
"Cloud_Water_Path_Uncertainty", "CWPU", |
|
| 57 |
"Cloud_Phase_Optical_Properties", "CPOP", |
|
| 58 |
"Cloud_Multi_Layer_Flag", "CMLF", |
|
| 59 |
"Cloud_Mask_1km", "CM1", |
|
| 60 |
"Quality_Assurance_1km", "QA"), |
|
| 61 |
byrow=T,ncol=2,dimnames=list(1:10,c("variable","varid"))),stringsAsFactors=F)
|
|
| 62 |
save(vars,file="vars.Rdata") |
|
| 63 |
|
|
| 8 |
save(tb,file="modlandTiles.Rdata") |
|
| 64 | 9 |
|
| 65 | 10 |
### Submission script |
| 66 | 11 |
|
| 67 | 12 |
cat(paste("
|
| 68 |
#PBS -S /bin/sh |
|
| 69 |
#PBS -J 700-899 |
|
| 70 |
###PBS -J 1-",length(notdone)," |
|
| 13 |
#PBS -S /bin/bash |
|
| 14 |
#PBS -l select=1:ncpus=4:mpiprocs=4:model=wes |
|
| 15 |
####old PBS -l select=64:ncpus=4:mpiprocs=4:model=wes |
|
| 16 |
####### old: select=48:ncpus=8:mpiprocs=8:model=neh |
|
| 71 | 17 |
#PBS -l walltime=0:10:00 |
| 72 |
#PBS -l ncpus=100 |
|
| 73 | 18 |
#PBS -j oe |
| 74 |
#PBS -o log/log_^array_index^ |
|
| 75 | 19 |
#PBS -m e |
| 20 |
#PBS -V |
|
| 21 |
####PBS -W group_list=s1007 |
|
| 22 |
#PBS -q devel |
|
| 23 |
###PBS -o log/log_^array_index^ |
|
| 24 |
#PBS -o log/log_DataCompile |
|
| 76 | 25 |
#PBS -M adam.wilson@yale.edu |
| 77 | 26 |
#PBS -N MOD06 |
| 78 | 27 |
|
| 28 |
source /usr/share/modules/init/bash |
|
| 29 |
|
|
| 79 | 30 |
## cd to working directory |
| 80 | 31 |
cd /nobackupp1/awilso10/mod06 |
| 81 | 32 |
|
| 82 | 33 |
## set some memory limits |
| 83 | 34 |
# ulimit -d 1500000 -m 1500000 -v 1500000 #limit memory usage |
| 84 | 35 |
source /usr/local/lib/global.profile |
| 36 |
source /u/awilso10/.bashrc |
|
| 85 | 37 |
## export a few important variables |
| 86 | 38 |
export PATH=$PATH:/nobackupp1/awilso10/bin:/nobackupp1/awilso10/software/bin |
| 87 | 39 |
export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/nobackupp1/awilso10/software/lib |
| 88 | 40 |
export MRTDATADIR=/nobackupp1/awilso10/software/heg/data |
| 89 | 41 |
export PGSHOME=/nobackupp1/awilso10/software/heg |
| 90 |
export MRTBINDIR=/nobackup1/awilso10/software/TOOLKIT_MTD |
|
| 42 |
export MRTBINDIR=/nobackupp1/awilso10/software/TOOLKIT_MTD
|
|
| 91 | 43 |
export R_LIBS=\"/u/awilso10/R/x86_64-unknown-linux-gnu-library/2.15/\" |
| 44 |
export TMPDIR=/nobackupp1/awilso10/mod06/tmp |
|
| 92 | 45 |
## load modules |
| 93 | 46 |
module load gcc mpi-sgi/mpt.2.06r6 hdf4 udunits R |
| 94 | 47 |
## Run the script! |
| 95 |
Rscript --verbose --vanilla /u/awilso10/environmental-layers/climate/procedures/MOD06_L2_data_compile_Pleiades.r i=${PBS_ARRAY_INDEX}
|
|
| 96 |
rm -r $TMPDIR |
|
| 48 |
TMPDIR=$TMPDIR Rscript --verbose --vanilla /u/awilso10/environmental-layers/climate/procedures/MOD06_L2_process.r |
|
| 97 | 49 |
exit 0 |
| 98 | 50 |
",sep=""),file="MOD06_process") |
| 99 | 51 |
|
| 100 | 52 |
### Check the file |
| 101 | 53 |
system("cat MOD06_process")
|
| 102 |
#system("chmod +x MOD06_process")
|
|
| 54 |
#system("cat ~/environmental-layers/climate/procedures/MOD06_L2_process.r")
|
|
| 55 |
|
|
| 56 |
## Submit it (and keep the pid)! |
|
| 57 |
pid=system("qsub -q devel MOD06_process",intern=T); pid; pid=strsplit(pid,split="[.]")[[1]][1]
|
|
| 58 |
|
|
| 59 |
#system("qsub MOD06_process")
|
|
| 103 | 60 |
|
| 104 |
## Submit it! |
|
| 105 |
#system("qsub -q devel MOD06_process")
|
|
| 106 |
system("qsub MOD06_process")
|
|
| 61 |
## work in interactive mode |
|
| 62 |
#system("qsub -I -lselect=1:ncpus=2:model=wes -q devel")
|
|
| 107 | 63 |
|
| 108 | 64 |
## check progress |
| 109 | 65 |
system("qstat -u awilso10")
|
| 110 |
system("qstat -t 391843[]")
|
|
| 111 |
system("qstat -f 391843[2]")
|
|
| 66 |
system(paste("qstat -t -x",pid))
|
|
| 112 | 67 |
|
| 113 |
#system("qstat devel ")
|
|
| 68 |
system("qstat devel ")
|
|
| 114 | 69 |
#system("qstat | grep awilso10")
|
| 115 | 70 |
|
| 116 |
#print(paste(max(0,length(system("qstat",intern=T))-2)," processes running"))
|
|
| 117 |
# system("ssh c0-8.farm.caes.ucdavis.edu")
|
|
| 118 |
# system("qalter -p +1024 25964") #decrease priority of job to run extraction below.
|
|
| 119 |
system("cat log/InterpScript.o55934.2")
|
|
| 120 | 71 |
|
| 121 |
## check log |
|
| 122 |
system(paste("cat",list.files("log",pattern="InterpScript",full=T)[100]))
|
|
| 123 |
#system(paste("cat",list.files("log",pattern="InterpScript",full=T)[13]," | grep \"Temporary Directory\""))
|
|
| 72 |
### copy the files back to Yale |
|
| 73 |
system("scp 2_daily/* adamw@acrobates.eeb.yale.edu:/data/personal/adamw/projects/interp/")
|
|
Also available in: Unified diff
Script is successfully running and producing the summary files, though the output looks strange. Maybe a problem with sinusoidal output of HEG? Use Pleiades.R to drive the submission and MOD06_L2_process as the processing script