| 1 |
35d59dc1
|
Adam M. Wilson @ pfe
|
### Process a folder of daily MOD35 HDF files to produce a climatology
|
| 2 |
|
|
|
| 3 |
aba23d60
|
Adam M. Wilson @ pfe
|
.libPaths("/pleiades/u/awilso10/R/x86_64-unknown-linux-gnu-library/2.15/")
|
| 4 |
|
|
|
| 5 |
35d59dc1
|
Adam M. Wilson @ pfe
|
## import commandline arguments
|
| 6 |
aba23d60
|
Adam M. Wilson @ pfe
|
library(getopt,lib="/pleiades/u/awilso10/R/x86_64-unknown-linux-gnu-library/2.15/")
|
| 7 |
35d59dc1
|
Adam M. Wilson @ pfe
|
## get options
|
| 8 |
|
|
opta <- getopt(matrix(c(
|
| 9 |
|
|
'tile', 't', 1, 'character',
|
| 10 |
|
|
'verbose','v',1,'logical'
|
| 11 |
|
|
), ncol=4, byrow=TRUE))
|
| 12 |
|
|
|
| 13 |
|
|
tile=opta$tile #tile="h11v08"
|
| 14 |
|
|
verbose=opta$verbose #print out extensive information for debugging?
|
| 15 |
|
|
|
| 16 |
aba23d60
|
Adam M. Wilson @ pfe
|
## set working directory
|
| 17 |
|
|
setwd("/u/awilso10/MOD35")
|
| 18 |
|
|
|
| 19 |
35d59dc1
|
Adam M. Wilson @ pfe
|
### directory containing daily files
|
| 20 |
|
|
outdir=paste("daily/",tile,"/",sep="") #directory for separate daily files
|
| 21 |
|
|
|
| 22 |
|
|
### directory to hold climatology
|
| 23 |
|
|
outdir2="summary" #directory for combined daily files and summarized files
|
| 24 |
|
|
if(!file.exists(outdir2)) dir.create(outdir2)
|
| 25 |
|
|
|
| 26 |
|
|
### path to NCO
|
| 27 |
|
|
ncopath="/nasa/sles11/nco/4.0.8/gcc/mpt/bin/"
|
| 28 |
|
|
|
| 29 |
|
|
### Vector of variables that must be in file or they will be deleted.
|
| 30 |
|
|
### Formated as output from system(paste("cdo -s showvar ",fdly$path[i]),intern=T)
|
| 31 |
|
|
#finalvars=" CER COT CLD"
|
| 32 |
|
|
|
| 33 |
|
|
|
| 34 |
|
|
################################################################################
|
| 35 |
|
|
## Get list of all daily files
|
| 36 |
be64daa8
|
Adam M. Wilson @ pfe
|
if(verbose) print(paste("Checking daily output in preparation for generating climatology:",tile))
|
| 37 |
35d59dc1
|
Adam M. Wilson @ pfe
|
|
| 38 |
|
|
fdly=data.frame(path=list.files(outdir,pattern="nc$",full=T),stringsAsFactors=F)
|
| 39 |
|
|
fdly$file=basename(fdly$path)
|
| 40 |
|
|
fdly$dateid=substr(fdly$file,14,21)
|
| 41 |
|
|
fdly$date=as.Date(substr(fdly$file,14,21),"%Y%m%d")
|
| 42 |
|
|
fdly$month=format(fdly$date,"%m")
|
| 43 |
|
|
fdly$year=format(fdly$date,"%Y")
|
| 44 |
|
|
nrow(fdly)
|
| 45 |
|
|
|
| 46 |
|
|
## print some summaries
|
| 47 |
be64daa8
|
Adam M. Wilson @ pfe
|
if(verbose) print(paste("Summary of available daily files:",tile))
|
| 48 |
35d59dc1
|
Adam M. Wilson @ pfe
|
print(table(fdly$year))
|
| 49 |
|
|
print(table(fdly$month))
|
| 50 |
|
|
#print(table(fdly$fvar))
|
| 51 |
|
|
|
| 52 |
|
|
#################################################################################
|
| 53 |
|
|
## Combine the year-by-year files into a single daily file in the summary directory (for archiving)
|
| 54 |
be64daa8
|
Adam M. Wilson @ pfe
|
if(verbose) print(paste("Merging daily files into single file output:",tile))
|
| 55 |
35d59dc1
|
Adam M. Wilson @ pfe
|
|
| 56 |
|
|
## create temporary directory to put intermediate files (will be deleted when R quits)
|
| 57 |
|
|
tsdir=paste(tempdir(),"/summary",sep="")
|
| 58 |
|
|
if(!file.exists(tsdir)) dir.create(tsdir,recursive=T)
|
| 59 |
|
|
|
| 60 |
|
|
## merge all daily files to create a single file with all dates
|
| 61 |
|
|
system(paste(ncopath,"ncrcat -O ",outdir,"/*nc ",outdir2,"/MOD35_",tile,"_daily.nc",sep=""))
|
| 62 |
1523a1e7
|
Adam M. Wilson
|
|
| 63 |
35d59dc1
|
Adam M. Wilson @ pfe
|
## Update attributes
|
| 64 |
|
|
system(paste(ncopath,"ncatted ",
|
| 65 |
c24e32a8
|
Adam M. Wilson @ pfe
|
" -a valid_min,PClear,o,b,0 ",
|
| 66 |
|
|
" -a valid_max,PClear,o,b,100 ",
|
| 67 |
|
|
#" -a valid_range,PClear,o,b,\"0,255\" ",
|
| 68 |
|
|
#" -a missing_value,PClear,o,b,255 ",
|
| 69 |
|
|
#" -a _FillValue,PClear,d,b,255 ",
|
| 70 |
35d59dc1
|
Adam M. Wilson @ pfe
|
" -a units,time,o,c,\"days since 2000-1-1 0:0:0\" ",
|
| 71 |
|
|
" -a title,global,o,c,\"MODIS Cloud Product (MOD35) Daily Timeseries\" ",
|
| 72 |
|
|
" -a institution,global,o,c,\"Yale University\" ",
|
| 73 |
|
|
" -a source,global,o,c,\"MODIS Cloud Mask (MOD35)\" ",
|
| 74 |
|
|
" -a comment,global,o,c,\"Compiled by Adam M. Wilson (adam.wilson@yale.edu)\" ",
|
| 75 |
|
|
outdir2,"/MOD35_",tile,"_daily.nc",sep=""))
|
| 76 |
|
|
|
| 77 |
78f94710
|
Adam M. Wilson
|
## report on daily file:
|
| 78 |
|
|
#ncfile=paste(outdir2,"/MOD35_",tile,"_daily.nc",sep="")
|
| 79 |
|
|
system(paste("ncdump -h ",outdir2,"/MOD35_",tile,"_daily.nc | head -20 ",sep=""))
|
| 80 |
|
|
system(paste("cdo showyear ",outdir2,"/MOD35_",tile,"_daily.nc",sep=""))
|
| 81 |
|
|
system(paste("cdo showmon ",outdir2,"/MOD35_",tile,"_daily.nc",sep=""))
|
| 82 |
|
|
#system(paste("cdo showdate ",outdir2,"/MOD35_",tile,"_daily.nc",sep=""))
|
| 83 |
|
|
|
| 84 |
|
|
|
| 85 |
35d59dc1
|
Adam M. Wilson @ pfe
|
### produce a monthly timeseries?
|
| 86 |
|
|
#system(paste("cdo -O monmean ",outdir2,"/MOD35_",tile,"_daily.nc ",tsdir,"/MOD35_",tile,"_monmean.nc",sep=""))
|
| 87 |
|
|
|
| 88 |
|
|
#############################
|
| 89 |
|
|
## Generate the Climatologies
|
| 90 |
be64daa8
|
Adam M. Wilson @ pfe
|
if(verbose) print(paste("Generate monthly climatologies: ",tile))
|
| 91 |
35d59dc1
|
Adam M. Wilson @ pfe
|
|
| 92 |
|
|
myear=as.integer(max(fdly$year)) #this year will be used in all dates of monthly climatologies (and day will = 15)
|
| 93 |
|
|
|
| 94 |
|
|
## Monthly means
|
| 95 |
be64daa8
|
Adam M. Wilson @ pfe
|
if(verbose) print(paste("Calculating the monthly means:",tile))
|
| 96 |
c24e32a8
|
Adam M. Wilson @ pfe
|
system(paste("cdo -O -b I8 -v sorttimestamp -setyear,",myear," -setday,15 -mulc,-1 -subc,100 -ymonmean ",outdir2,"/MOD35_",tile,"_daily.nc ",tsdir,"/MOD35_",tile,"_ymonmean.nc",sep=""),wait=T)
|
| 97 |
|
|
system(paste(ncopath,"ncrename -v PClear,PCloud ",tsdir,"/MOD35_",tile,"_ymonmean.nc",sep=""))
|
| 98 |
|
|
system(paste(ncopath,"ncatted ",
|
| 99 |
|
|
" -a long_name,PCloud,o,c,\"Mean Probability of Cloud\" ",
|
| 100 |
|
|
" -a missing_value,PCloud,o,b,255 ",
|
| 101 |
|
|
" -a _FillValue,PCloud,d,b,255 ",
|
| 102 |
|
|
tsdir,"/MOD35_",tile,"_ymonmean.nc",sep=""))
|
| 103 |
b83b17ae
|
Adam M. Wilson
|
|
| 104 |
|
|
|
| 105 |
c24e32a8
|
Adam M. Wilson @ pfe
|
## Monthly Mean
|
| 106 |
|
|
#months=c("01","02","03","04","05","06","07","08","09","10","11","12")
|
| 107 |
|
|
# month="02"
|
| 108 |
|
|
|
| 109 |
aba23d60
|
Adam M. Wilson @ pfe
|
#system(paste("cdo -O sorttimestamp -setyear,",myear," -mulc,-1 -subc,100 -ydrunmean,30 ",outdir2,"/MOD35_",tile,"_daily.nc ",outdir2,"/MOD35_",tile,"_ydrunmean30.nc &",sep=""),wait=T)
|
| 110 |
c24e32a8
|
Adam M. Wilson @ pfe
|
#system(paste("scp summary/MOD35_",tile,".nc adamw@acrobates.eeb.yale.edu:/data/personal/adamw/projects/interp/data/modis/mod35/",sep=""))
|
| 111 |
|
|
#system(paste("ncdump -h ",tsdir,"/MOD35_",tile,"_ymonmean.nc ",sep=""))
|
| 112 |
b83b17ae
|
Adam M. Wilson
|
|
| 113 |
35d59dc1
|
Adam M. Wilson @ pfe
|
|
| 114 |
|
|
## Monthly standard deviation
|
| 115 |
be64daa8
|
Adam M. Wilson @ pfe
|
if(verbose) print(paste("Calculating the monthly SD:",tile))
|
| 116 |
|
|
system(paste("cdo -O -b I8 sorttimestamp -setyear,",myear," -setday,15 -ymonstd -mulc,-1 -subc,100 -monmean ",
|
| 117 |
aba23d60
|
Adam M. Wilson @ pfe
|
outdir2,"/MOD35_",tile,"_daily.nc ",
|
| 118 |
|
|
tsdir,"/MOD35_",tile,"_ymonstd.nc",sep=""))
|
| 119 |
1523a1e7
|
Adam M. Wilson
|
system(paste(ncopath,"ncrename -v PClear,PCloud_sd ",tsdir,"/MOD35_",tile,"_ymonstd.nc",sep=""))
|
| 120 |
35d59dc1
|
Adam M. Wilson @ pfe
|
system(paste(ncopath,"ncatted ",
|
| 121 |
1523a1e7
|
Adam M. Wilson
|
" -a long_name,PCloud_sd,o,c,\"Standard Deviation of p(cloud)\" ",
|
| 122 |
35d59dc1
|
Adam M. Wilson @ pfe
|
tsdir,"/MOD35_",tile,"_ymonstd.nc",sep=""))
|
| 123 |
|
|
|
| 124 |
1523a1e7
|
Adam M. Wilson
|
## frequency of cloud days p(clear<90%)
|
| 125 |
be64daa8
|
Adam M. Wilson @ pfe
|
if(verbose) print(paste("Calculating the proportion of cloudy and probably cloudy days:",tile))
|
| 126 |
c24e32a8
|
Adam M. Wilson @ pfe
|
system(paste("cdo -O -b I8 sorttimestamp -setyear,",myear," -setday,15 -ymonmean -mulc,100 -lec,90 -selvar,PClear ",outdir2,"/MOD35_",tile,"_daily.nc ",tsdir,"/MOD35_",tile,"_ymoncld01.nc",sep=""))
|
| 127 |
1523a1e7
|
Adam M. Wilson
|
system(paste(ncopath,"ncrename -v PClear,CF ",tsdir,"/MOD35_",tile,"_ymoncld01.nc",sep=""))
|
| 128 |
35d59dc1
|
Adam M. Wilson @ pfe
|
system(paste(ncopath,"ncatted ",
|
| 129 |
1523a1e7
|
Adam M. Wilson
|
" -a long_name,CF,o,c,\"Cloud Frequency: Proportion of Days with probability of clear < 90%\" ",
|
| 130 |
|
|
" -a units,CF,o,c,\"Proportion (%)\" ",
|
| 131 |
35d59dc1
|
Adam M. Wilson @ pfe
|
tsdir,"/MOD35_",tile,"_ymoncld01.nc",sep=""))
|
| 132 |
|
|
|
| 133 |
|
|
## number of observations
|
| 134 |
be64daa8
|
Adam M. Wilson @ pfe
|
if(verbose) print(paste("Calculating the number of missing variables:",tile))
|
| 135 |
c24e32a8
|
Adam M. Wilson @ pfe
|
system(paste("cdo -O -b I8 sorttimestamp -setyear,",myear," -setday,15 -ymonmean -mulc,100 -eqc,9999 -setmisstoc,9999 -selvar,PClear ",outdir2,"/MOD35_",tile,"_daily.nc ",tsdir,"/MOD35_",tile,"_ymonmiss.nc",sep=""))
|
| 136 |
78f94710
|
Adam M. Wilson
|
system(paste(ncopath,"ncrename -v PClear,Pmiss ",tsdir,"/MOD35_",tile,"_ymonmiss.nc",sep=""))
|
| 137 |
35d59dc1
|
Adam M. Wilson @ pfe
|
system(paste(ncopath,"ncatted ",
|
| 138 |
78f94710
|
Adam M. Wilson
|
" -a long_name,Pmiss,o,c,\"Proportion of Days with missing data\" ",
|
| 139 |
|
|
" -a units,Pmiss,o,c,\"Proportion (%)\" ",
|
| 140 |
|
|
tsdir,"/MOD35_",tile,"_ymonmiss.nc",sep=""))
|
| 141 |
35d59dc1
|
Adam M. Wilson @ pfe
|
|
| 142 |
|
|
## TODO: fix projection information so GDAL can read it correctly.
|
| 143 |
|
|
## clean up variables?
|
| 144 |
|
|
|
| 145 |
|
|
## append variables to a single file
|
| 146 |
be64daa8
|
Adam M. Wilson @ pfe
|
if(verbose) print(paste("Append all monthly climatologies into a single file:",tile))
|
| 147 |
35d59dc1
|
Adam M. Wilson @ pfe
|
system(paste(ncopath,"ncks -O ",tsdir,"/MOD35_",tile,"_ymonmean.nc ",tsdir,"/MOD35_",tile,"_ymon.nc",sep=""))
|
| 148 |
|
|
system(paste(ncopath,"ncks -A ",tsdir,"/MOD35_",tile,"_ymonstd.nc ",tsdir,"/MOD35_",tile,"_ymon.nc",sep=""))
|
| 149 |
|
|
system(paste(ncopath,"ncks -A ",tsdir,"/MOD35_",tile,"_ymoncld01.nc ",tsdir,"/MOD35_",tile,"_ymon.nc",sep=""))
|
| 150 |
|
|
system(paste(ncopath,"ncks -A ",tsdir,"/MOD35_",tile,"_ymonmiss.nc ",tsdir,"/MOD35_",tile,"_ymon.nc",sep=""))
|
| 151 |
|
|
|
| 152 |
|
|
## append sinusoidal grid from one of input files as CDO doesn't transfer all attributes
|
| 153 |
be64daa8
|
Adam M. Wilson @ pfe
|
if(verbose) print(paste("Clean up file (update attributes, flip latitudes, add grid description:",tile))
|
| 154 |
35d59dc1
|
Adam M. Wilson @ pfe
|
|
| 155 |
|
|
#system(paste(ncopath,"ncea -d time,0,1 -v sinusoidal ",list.files(outdir,full=T,pattern="[.]nc$")[1]," ",tsdir,"/sinusoidal.nc",sep=""))
|
| 156 |
|
|
#system(paste(ncopath,"ncks -A -d time,0,1 -v sinusoidal ",list.files(outdir,full=T,pattern="[.]nc$")[1]," ",tsdir,"/MOD35_",tile,"_ymon.nc",sep=""))
|
| 157 |
|
|
|
| 158 |
|
|
## invert latitude so it plays nicely with gdal
|
| 159 |
|
|
system(paste(ncopath,"ncpdq -O -a -y ",tsdir,"/MOD35_",tile,"_ymon.nc ",outdir2,"/MOD35_",tile,".nc",sep=""))
|
| 160 |
|
|
|
| 161 |
b716a3f0
|
Adam M. Wilson @ pfe
|
## proj string taken from GDAL-written MODIS tile
|
| 162 |
|
|
projstring="PROJCS[\"Sinusoidal (Sanson-Flamsteed)\",GEOGCS[\"wgs84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS_1984\",6378137,298.257223563],TOWGS84[0,0,0,0,0,0,0]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433]],PROJECTION[\"Sinusoidal\"],PARAMETER[\"longitude_of_center\",0],PARAMETER[\"false_easting\",0],PARAMETER[\"false_northing\",0],UNIT[\"Meter\",1]]"
|
| 163 |
|
|
|
| 164 |
35d59dc1
|
Adam M. Wilson @ pfe
|
## update attributes
|
| 165 |
|
|
system(paste(ncopath,"ncatted ",
|
| 166 |
b716a3f0
|
Adam M. Wilson @ pfe
|
" -a false_easting,sinusoidal,o,d,0. ",
|
| 167 |
|
|
" -a false_northing,sinusoidal,o,d,0. ",
|
| 168 |
|
|
" -a longitude_of_central_meridian,sinusoidal,o,d,0. ",
|
| 169 |
|
|
" -a longitude_of_prime_meridian,sinusoidal,o,d,0. ",
|
| 170 |
|
|
" -a semi_major_axis,sinusoidal,o,d,6378137. ",
|
| 171 |
|
|
" -a inverse_flattening,sinusoidal,o,d,298.257223563 ",
|
| 172 |
|
|
" -a spatial_ref,sinusoidal,o,c,",projstring,
|
| 173 |
|
|
" -a GeoTransform,sinusoidal,o,c,\"-7783653.638366 926.6254331391661 0 1111950.519767 0 -926.6254331391667\" ",
|
| 174 |
35d59dc1
|
Adam M. Wilson @ pfe
|
" -a units,time,o,c,\"days since 2000-1-1 0:0:0\" ",
|
| 175 |
|
|
" -a title,global,o,c,\"MODIS Cloud Product (MOD35) Climatology\" ",
|
| 176 |
b716a3f0
|
Adam M. Wilson @ pfe
|
" -a tile,global,o,c,\"",tile,"\" ",
|
| 177 |
35d59dc1
|
Adam M. Wilson @ pfe
|
" -a institution,global,o,c,\"Yale University\" ",
|
| 178 |
1523a1e7
|
Adam M. Wilson
|
" -a source,global,o,c,\"MODIS Cloud Product (MOD35) Collection 6\" ",
|
| 179 |
35d59dc1
|
Adam M. Wilson @ pfe
|
" -a comment,global,o,c,\"Compiled by Adam M. Wilson (adam.wilson@yale.edu)\" ",
|
| 180 |
|
|
outdir2,"/MOD35_",tile,".nc",sep=""))
|
| 181 |
|
|
|
| 182 |
|
|
|
| 183 |
|
|
print(paste("############################### Processed ",nrow(fdly),"days for tile:",tile," ###################################################"))
|
| 184 |
|
|
print("Years:")
|
| 185 |
|
|
print(table(fdly$fyear))
|
| 186 |
|
|
print("Months:")
|
| 187 |
|
|
print(table(fdly$fmonth))
|
| 188 |
1523a1e7
|
Adam M. Wilson
|
|
| 189 |
35d59dc1
|
Adam M. Wilson @ pfe
|
## quit R
|
| 190 |
|
|
q("no")
|
| 191 |
|
|
|
