Revision 979e2d4c
Added by Adam M. Wilson over 12 years ago
| climate/procedures/MOD06_L2_data_compile.r | ||
|---|---|---|
| 49 | 49 |
fs$year=format(fs$date,"%Y") |
| 50 | 50 |
fs$time=substr(fs$file,19,22) |
| 51 | 51 |
fs$datetime=as.POSIXct(strptime(paste(substr(fs$file,11,17),substr(fs$file,19,22)), '%Y%j %H%M')) |
| 52 |
fs$dateid=format(fs$date,"%Y%m%d") |
|
| 52 | 53 |
fs$path=as.character(fs$path) |
| 53 | 54 |
fs$file=as.character(fs$file) |
| 54 | 55 |
|
| ... | ... | |
| 61 | 62 |
|
| 62 | 63 |
## vars |
| 63 | 64 |
vars=as.data.frame(matrix(c( |
| 64 |
"Cloud_Effective_Radius","CER", |
|
| 65 |
"Cloud_Effective_Radius_Uncertainty","CERU", |
|
| 66 |
"Cloud_Optical_Thickness","COT", |
|
| 67 |
"Cloud_Optical_Thickness_Uncertainty","COTU", |
|
| 68 |
"Cloud_Water_Path","CWP", |
|
| 69 |
"Cloud_Water_Path_Uncertainty","CWPU", |
|
| 70 |
"Cloud_Phase_Optical_Properties","CPOP", |
|
| 71 |
"Cloud_Multi_Layer_Flag","CMLF", |
|
| 72 |
"Cloud_Mask_1km","CM1", |
|
| 73 |
"Quality_Assurance_1km","QA"), |
|
| 65 |
"Cloud_Effective_Radius", "CER",
|
|
| 66 |
"Cloud_Effective_Radius_Uncertainty", "CERU",
|
|
| 67 |
"Cloud_Optical_Thickness", "COT",
|
|
| 68 |
"Cloud_Optical_Thickness_Uncertainty", "COTU",
|
|
| 69 |
"Cloud_Water_Path", "CWP",
|
|
| 70 |
"Cloud_Water_Path_Uncertainty", "CWPU",
|
|
| 71 |
"Cloud_Phase_Optical_Properties", "CPOP",
|
|
| 72 |
"Cloud_Multi_Layer_Flag", "CMLF",
|
|
| 73 |
"Cloud_Mask_1km", "CM1",
|
|
| 74 |
"Quality_Assurance_1km", "QA"),
|
|
| 74 | 75 |
byrow=T,ncol=2,dimnames=list(1:10,c("variable","varid"))))
|
| 75 | 76 |
|
| 76 | 77 |
|
| ... | ... | |
| 99 | 100 |
OUTPUT_PIXEL_SIZE_Y=1000 |
| 100 | 101 |
SPATIAL_SUBSET_UL_CORNER = ( ",upleft," ) |
| 101 | 102 |
SPATIAL_SUBSET_LR_CORNER = ( ",lowright," ) |
| 102 |
#RESAMPLING_TYPE =",ifelse(grepl("Flag|Mask|Quality",vars$variable),"NN","BICUBIC"),"
|
|
| 103 |
#RESAMPLING_TYPE =",ifelse(grepl("Flag|Mask|Quality",vars$variable),"NN","CUBIC"),"
|
|
| 103 | 104 |
RESAMPLING_TYPE =NN |
| 104 | 105 |
OUTPUT_PROJECTION_TYPE = SIN |
| 105 | 106 |
ELLIPSOID_CODE = WGS84 |
| 106 |
OUTPUT_PROJECTION_PARAMETERS = ( 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 0.0 )
|
|
| 107 |
OUTPUT_PROJECTION_PARAMETERS = ( 6371007.181 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 86400.0 0.0 0.0 )
|
|
| 107 | 108 |
OUTPUT_TYPE = HDFEOS |
| 108 | 109 |
OUTPUT_FILENAME = ",outfile," |
| 109 | 110 |
END |
| ... | ... | |
| 142 | 143 |
############################################################## |
| 143 | 144 |
### Import to GRASS for processing |
| 144 | 145 |
|
| 145 |
fs$grass=paste(fs$month,fs$year,fs$file,sep="_") |
|
| 146 |
#fs$grass=paste(fs$month,fs$year,fs$file,sep="_")
|
|
| 146 | 147 |
td=readGDAL(paste("HDF4_EOS:EOS_GRID:\"",outdir,"/",fs$file[1],"\":mod06:Cloud_Mask_1km_0",sep=""))
|
| 148 |
projection(td)="+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +ellps=WGS84 +units=m +no_defs" |
|
| 147 | 149 |
|
| 148 | 150 |
## fucntion to convert binary to decimal to assist in identifying correct values |
| 149 | 151 |
b2d=function(x) sum(x * 2^(rev(seq_along(x)) - 1)) #http://tolstoy.newcastle.edu.au/R/e2/help/07/02/10596.html |
| ... | ... | |
| 157 | 159 |
|
| 158 | 160 |
initGRASS(gisBase="/usr/lib/grass64",gisDbase=gisDbase,SG=td,override=T, |
| 159 | 161 |
location=gisLocation,mapset=gisMapset) |
| 162 |
getLocationProj() |
|
| 160 | 163 |
|
| 161 |
|
|
| 162 |
## update projection (for some reason the datum doesn't get identified from the HDF files |
|
| 163 |
#cat("name: Sinusoidal (Sanson-Flamsteed)
|
|
| 164 |
#proj: sinu |
|
| 165 |
#datum: wgs84 |
|
| 166 |
#ellps: wgs84 |
|
| 167 |
#lon_0: 0 |
|
| 168 |
#x_0: 0 |
|
| 169 |
#y_0: 0 |
|
| 170 |
#towgs84: 0,0,0,0,0,0,0 |
|
| 171 |
#no_defs: defined",file=paste(gisDbase,"/",gisLocation,"/PERMANENT/PROJ_INFO",sep="")) |
|
| 172 |
## update PROJ_UNITS |
|
| 173 |
#cat("unit: Meter
|
|
| 174 |
#units: Meters |
|
| 175 |
#meters: 1",file=paste(gisDbase,"/",gisLocation,"/PERMANENT/PROJ_UNITS",sep="")) |
|
| 176 |
#system("g.proj -d")
|
|
| 164 |
system("g.mapset mapset=PERMANENT")
|
|
| 165 |
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",outdir,"/",fs$file[1],"\":mod06:Cloud_Mask_1km_0",sep=""),
|
|
| 166 |
output="modisgrid",flags=c("quiet","overwrite","o"))
|
|
| 167 |
system("g.region rast=modisgrid save=roi --overwrite")
|
|
| 168 |
system(paste("g.proj -c proj4=\"+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +ellps=WGS84 +units=m +no_defs\"",sep=""))
|
|
| 169 |
system("g.region roi")
|
|
| 170 |
system("g.region -p")
|
|
| 177 | 171 |
|
| 178 | 172 |
i=1 |
| 179 | 173 |
file=paste(outdir,"/",fs$file[1],sep="") |
| 180 |
|
|
| 181 |
|
|
| 182 |
loadcloud<-function(file){
|
|
| 183 |
## Cloud Mask |
|
| 184 |
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod06:Cloud_Mask_1km_0",sep=""),
|
|
| 185 |
output="CM1",flags=c("overwrite","o"))
|
|
| 186 |
system("g.region rast=CM1")
|
|
| 187 |
|
|
| 188 |
## extract cloudy and 'confidently clear' pixels |
|
| 189 |
system(paste("r.mapcalc <<EOF
|
|
| 190 |
CM_cloud= ((CM1 / 2^0) % 2) == 1 && ((CM1 / 2^1) % 2^2) == 0 |
|
| 191 |
CM_clear= ((CM1 / 2^0) % 2) == 1 && ((CM1 / 2^1) % 2^2) == 3 |
|
| 192 |
EOF")) |
|
| 193 |
|
|
| 194 |
## QA |
|
| 174 |
date=as.Date("2000-03-02")
|
|
| 175 |
|
|
| 176 |
loadcloud<-function(date,fs){
|
|
| 177 |
sink(file=paste("output/",format(date,"%Y%m%d"),".txt",sep=""),type="output")
|
|
| 178 |
## Identify which files to process |
|
| 179 |
tfs=fs$file[fs$date==date] |
|
| 180 |
nfs=length(tfs) |
|
| 181 |
|
|
| 182 |
## create new mapset to hold all data for this day |
|
| 183 |
system(paste("g.mapset --quiet -c mapset=",gisMapset,"_",format(date,"%Y%m%d"),sep=""))
|
|
| 184 |
# file.copy(paste(gisDbase,"/",gisLocation,"/PERMANENT/DEFAULT_WIND",sep=""),paste(gisDbase,"/",gisLocation,"/",gisMapset,"_",format(date,"%Y%m%d"),"/WIND",sep="")) |
|
| 185 |
system("g.region roi@PERMANENT")
|
|
| 186 |
|
|
| 187 |
## loop through scenes and process QA flags |
|
| 188 |
for(i in 1:nfs){
|
|
| 189 |
file=paste(outdir,"/",tfs[i],sep="") |
|
| 190 |
## Cloud Mask |
|
| 191 |
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod06:Cloud_Mask_1km_0",sep=""),
|
|
| 192 |
output=paste("CM1_",i,sep=""),flags=c("overwrite","o","quiet")) ; print("")
|
|
| 193 |
## extract cloudy and 'confidently clear' pixels |
|
| 194 |
system(paste("r.mapcalc <<EOF
|
|
| 195 |
CM_cloud_",i," = ((CM1_",i," / 2^0) % 2) == 1 && ((CM1_",i," / 2^1) % 2^2) == 0 |
|
| 196 |
CM_clear_",i," = ((CM1_",i," / 2^0) % 2) == 1 && ((CM1_",i," / 2^1) % 2^2) == 3 |
|
| 197 |
EOF",sep="")) |
|
| 198 |
|
|
| 199 |
|
|
| 200 |
## QA |
|
| 195 | 201 |
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod06:Quality_Assurance_1km_0",sep=""),
|
| 196 |
output="QA",flags=c("overwrite","o"))
|
|
| 202 |
output=paste("QA_",i,sep=""),flags=c("overwrite","o","quiet")) ; print("")
|
|
| 197 | 203 |
## QA_CER |
| 198 | 204 |
system(paste("r.mapcalc <<EOF
|
| 199 |
QA_COT= ((QA / 2^0) % 2^1 )==1
|
|
| 200 |
QA_COT2= ((QA / 2^1) % 2^2 )==3
|
|
| 201 |
QA_COT3= ((QA / 2^3) % 2^2 )==0
|
|
| 202 |
QA_CER= ((QA / 2^5) % 2^1 )==1
|
|
| 203 |
QA_CER2= ((QA / 2^6) % 2^2 )==3
|
|
| 204 |
EOF")) |
|
| 205 |
# QA_CWP= ((QA / 2^8) % 2^1 )==1
|
|
| 206 |
# QA_CWP2= ((QA / 2^9) % 2^2 )==3
|
|
| 205 |
QA_COT_",i,"= ((QA_",i," / 2^0) % 2^1 )==1
|
|
| 206 |
QA_COT2_",i,"= ((QA_",i," / 2^1) % 2^2 )==3
|
|
| 207 |
QA_COT3_",i,"= ((QA_",i," / 2^3) % 2^2 )==0
|
|
| 208 |
QA_CER_",i,"= ((QA_",i," / 2^5) % 2^1 )==1
|
|
| 209 |
QA_CER2_",i,"= ((QA_",i," / 2^6) % 2^2 )==3
|
|
| 210 |
EOF",sep=""))
|
|
| 211 |
# QA_CWP_",i,"= ((QA_",i," / 2^8) % 2^1 )==1
|
|
| 212 |
# QA_CWP2_",i,"= ((QA_",i," / 2^9) % 2^2 )==3
|
|
| 207 | 213 |
|
| 208 | 214 |
## Optical Thickness |
| 209 | 215 |
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod06:Cloud_Optical_Thickness",sep=""),
|
| 210 |
output="COT",
|
|
| 216 |
output=paste("COT_",i,sep=""),
|
|
| 211 | 217 |
title="cloud_effective_radius", |
| 212 | 218 |
flags=c("overwrite","o")) ; print("")
|
| 213 |
execGRASS("r.null",map="COT",setnull="-9999")
|
|
| 219 |
execGRASS("r.null",map=paste("COT_",i,sep=""),setnull="-9999")
|
|
| 214 | 220 |
## keep only positive COT values where quality is 'useful' and 'very good' & scale to real units |
| 215 |
system(paste("r.mapcalc \"COT=if(QA_COT&&QA_COT2&&QA_COT3&&COT>=0,COT*0.009999999776482582,null())\""))
|
|
| 221 |
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=""))
|
|
| 216 | 222 |
## set COT to 0 in clear-sky pixels |
| 217 |
system(paste("r.mapcalc \"COT2=if(CM_clear,COT,0)\""))
|
|
| 223 |
system(paste("r.mapcalc \"COT2_",i,"=if(CM_clear_",i,"==0,COT_",i,",0)\"",sep=""))
|
|
| 218 | 224 |
|
| 219 |
## Effective radius |
|
| 225 |
## Effective radius ##
|
|
| 220 | 226 |
execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod06:Cloud_Effective_Radius",sep=""),
|
| 221 |
output="CER",
|
|
| 227 |
output=paste("CER_",i,sep=""),
|
|
| 222 | 228 |
title="cloud_effective_radius", |
| 223 | 229 |
flags=c("overwrite","o")) ; print("")
|
| 224 |
execGRASS("r.null",map="CER",setnull="-9999")
|
|
| 230 |
execGRASS("r.null",map=paste("CER_",i,sep=""),setnull="-9999")
|
|
| 225 | 231 |
## keep only positive CER values where quality is 'useful' and 'very good' & scale to real units |
| 226 |
system(paste("r.mapcalc \"CER=if(QA_CER&&QA_CER2&&CER>=0,CER*0.009999999776482582,null())\""))
|
|
| 232 |
system(paste("r.mapcalc \"CER_",i,"=if(QA_CER_",i,"&&QA_CER2_",i,"&&CER_",i,">=0,CER_",i,"*0.009999999776482582,null())\"",sep=""))
|
|
| 233 |
## set CER to 0 in clear-sky pixels |
|
| 234 |
system(paste("r.mapcalc \"CER2_",i,"=if(CM_clear_",i,"==0,CER_",i,",0)\"",sep=""))
|
|
| 227 | 235 |
|
| 228 | 236 |
## Cloud Water Path |
| 229 | 237 |
# execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod06:Cloud_Water_Path",sep=""),
|
| ... | ... | |
| 232 | 240 |
# execGRASS("r.null",map="CWP",setnull="-9999")
|
| 233 | 241 |
# ## keep only positive CWP values where quality is 'useful' and 'very good' & scale to real units |
| 234 | 242 |
# system(paste("r.mapcalc \"CWP=if(QA_CWP&&QA_CWP2,CWP,null())\""))
|
| 235 |
|
|
| 236 | 243 |
|
| 237 |
} ; |
|
| 244 |
} #end loop through sub daily files |
|
| 245 |
|
|
| 246 |
#### Now generate daily averages |
|
| 247 |
|
|
| 248 |
system(paste("r.mapcalc <<EOF
|
|
| 249 |
COT_denom=",paste("!isnull(COT2_",1:nfs,")",sep="",collapse="+"),"
|
|
| 250 |
COT_numer=",paste("if(isnull(COT2_",1:nfs,"),0,COT2_",1:nfs,")",sep="",collapse="+"),"
|
|
| 251 |
COT_daily=COT_numer/COT_denom |
|
| 252 |
CER_denom=",paste("!isnull(CER2_",1:nfs,")",sep="",collapse="+"),"
|
|
| 253 |
CER_numer=",paste("if(isnull(CER2_",1:nfs,"),0,CER2_",1:nfs,")",sep="",collapse="+"),"
|
|
| 254 |
CER_daily=CER_numer/CER_denom |
|
| 255 |
EOF",sep="")) |
|
| 256 |
|
|
| 257 |
#### Write the file to a geotiff |
|
| 258 |
execGRASS("r.out.gdal",input="CER_daily",output=paste("data/modis/MOD06_L2_tif/CER_",format(date,"%Y%m%d"),".tif",sep=""),nodata=-999)
|
|
| 259 |
execGRASS("r.out.gdal",input="COT_daily",output=paste("data/modis/MOD06_L2_tif/COT_",format(date,"%Y%m%d"),".tif",sep=""),nodata=-999)
|
|
| 260 |
|
|
| 261 |
} |
|
| 262 |
|
|
| 263 |
|
|
| 264 |
########################################### |
|
| 265 |
### Now run it |
|
| 266 |
|
|
| 267 |
tdates=sort(unique(fs$date)) |
|
| 268 |
|
|
| 269 |
|
|
| 270 |
mclapply(tdates[1:10],loadcloud,fs=fs,mc.cores=1) |
|
| 238 | 271 |
|
| 272 |
|
|
| 273 |
|
|
| 274 |
|
|
| 275 |
# copy all datasets back to master mapset for summaries |
|
| 276 |
|
|
| 277 |
execGRASS("g.copy",rast=paste("\"COT_daily\",\"COT_",format(date,"%Y%m%d"),"@mod06\"",sep=""))
|
|
| 278 |
|
|
| 279 |
|
|
| 239 | 280 |
## unlock the grass database |
| 240 | 281 |
unlink_.gislock() |
| 241 | 282 |
|
Also available in: Unified diff
MOD06 processing runs well in sequence, but not parallel