| 97 |
97 |
|
| 98 |
98 |
## Monthly means
|
| 99 |
99 |
if(verbose) print("Calculating the monthly means")
|
| 100 |
|
system(paste("cdo -O sorttimestamp -setyear,",myear," -setday,15 -ymonmean ",outdir2,"/MOD35_",tile,"_daily.nc ",tsdir,"/MOD35_",tile,"_ymonmean.nc",sep=""),wait=T)
|
|
100 |
system(paste("cdo -O sorttimestamp -setyear,",myear," -setday,15 -mulc,-1 -subc,100 -ymonmean ",outdir2,"/MOD35_",tile,"_daily.nc ",tsdir,"/MOD35_",tile,"_ymonmean.nc",sep=""),wait=T)
|
|
101 |
system(paste(ncopath,"ncrename -v PClear,PCloud ",tsdir,"/MOD35_",tile,"_ymonmean.nc",sep=""))
|
|
102 |
system(paste(ncopath,"ncatted ",
|
|
103 |
" -a long_name,PCloud,o,c,\"Mean Probability of Cloud\" ",
|
|
104 |
tsdir,"/MOD35_",tile,"_ymonmean.nc",sep=""))
|
| 101 |
105 |
|
| 102 |
106 |
## Monthly standard deviation
|
| 103 |
107 |
if(verbose) print("Calculating the monthly SD")
|
| 104 |
|
system(paste("cdo -O sorttimestamp -setyear,",myear," -setday,15 -ymonstd ",outdir2,"/MOD35_",tile,"_daily.nc ",tsdir,"/MOD35_",tile,"_ymonstd.nc",sep=""))
|
| 105 |
|
system(paste(ncopath,"ncrename -v PClear,PClear_sd ",tsdir,"/MOD35_",tile,"_ymonstd.nc",sep=""))
|
|
108 |
system(paste("cdo -O sorttimestamp -setyear,",myear," -setday,15 -ymonstd ",outdir2,"/MOD35_",tile,"_daily.nc ",tsdir,"/MOD35_",tile,"_ymonstd.nc",sep=""))
|
|
109 |
system(paste(ncopath,"ncrename -v PClear,PCloud_sd ",tsdir,"/MOD35_",tile,"_ymonstd.nc",sep=""))
|
| 106 |
110 |
system(paste(ncopath,"ncatted ",
|
| 107 |
|
" -a long_name,PClear_sd,o,c,\"Standard Deviation of p(clear)\" ",
|
|
111 |
" -a long_name,PCloud_sd,o,c,\"Standard Deviation of p(cloud)\" ",
|
| 108 |
112 |
tsdir,"/MOD35_",tile,"_ymonstd.nc",sep=""))
|
| 109 |
113 |
|
| 110 |
114 |
## frequency of cloud days p(clear<90%)
|
| 111 |
115 |
if(verbose) print("Calculating the proportion of cloudy and probably cloudy days")
|
| 112 |
|
system(paste("cdo -O sorttimestamp -setyear,",myear," -setday,15 -nint -mulc,100 -ymonmean -lec,90 -setrtomiss,100,255 -selvar,PClear ",outdir2,"/MOD35_",tile,"_daily.nc ",tsdir,"/MOD35_",tile,"_ymoncld01.nc",sep=""))
|
|
116 |
system(paste("cdo -O sorttimestamp -setyear,",myear," -setday,15 -ymonmean -mulc,100 -lec,90 -selvar,PClear ",outdir2,"/MOD35_",tile,"_daily.nc ",tsdir,"/MOD35_",tile,"_ymoncld01.nc",sep=""))
|
| 113 |
117 |
system(paste(ncopath,"ncrename -v PClear,CF ",tsdir,"/MOD35_",tile,"_ymoncld01.nc",sep=""))
|
| 114 |
118 |
system(paste(ncopath,"ncatted ",
|
| 115 |
|
" -a long_name,CF,o,c,\"Proportion of Days with probability of clear < 90%\" ",
|
| 116 |
|
" -a units,CF,o,c,\"Proportion\" ",
|
|
119 |
" -a long_name,CF,o,c,\"Cloud Frequency: Proportion of Days with probability of clear < 90%\" ",
|
|
120 |
" -a units,CF,o,c,\"Proportion (%)\" ",
|
| 117 |
121 |
tsdir,"/MOD35_",tile,"_ymoncld01.nc",sep=""))
|
| 118 |
122 |
|
| 119 |
123 |
## number of observations
|
| 120 |
|
#if(verbose) print("Calculating the number of missing variables")
|
| 121 |
|
#system(paste("cdo -O sorttimestamp -setyear,",myear," -setday,15 -nint -mulc,100 -ymonmean -mulc,1.0 -eqc,9999 -setmisstoc,9999 -selvar,CLD ",outdir2,"/MOD35_",tile,"_daily.nc ",tsdir,"/MOD35_",tile,"_ymonmiss.nc",sep=""))
|
| 122 |
|
#system(paste(ncopath,"ncrename -v CLD,CLD_pmiss ",tsdir,"/MOD35_",tile,"_ymonmiss.nc",sep=""))
|
| 123 |
|
#system(paste(ncopath,"ncatted ",
|
| 124 |
|
#" -a long_name,CLD_pmiss,o,c,\"Proportion of Days with missing data for CLD\" ",
|
| 125 |
|
#" -a scale_factor,CLD_pmiss,o,d,0.01 ",
|
| 126 |
|
#" -a units,CLD_pmiss,o,c,\"Proportion\" ",
|
| 127 |
|
#tsdir,"/MOD35_",tile,"_ymonmiss.nc",sep=""))
|
|
124 |
if(verbose) print("Calculating the number of missing variables")
|
|
125 |
system(paste("cdo -O sorttimestamp -setyear,",myear," -setday,15 -nint -ymonmean -mulc,100 -eqc,9999 -setmisstoc,9999 -selvar,CLD ",outdir2,"/MOD35_",tile,"_daily.nc ",tsdir,"/MOD35_",tile,"_ymonmiss.nc",sep=""))
|
|
126 |
system(paste(ncopath,"ncrename -v PClear,CF_pmiss ",tsdir,"/MOD35_",tile,"_ymonmiss.nc",sep=""))
|
|
127 |
system(paste(ncopath,"ncatted ",
|
|
128 |
" -a long_name,CF_pmiss,o,c,\"Proportion of Days with missing data for CF\" ",
|
|
129 |
" -a units,CF_pmiss,o,c,\"Proportion (%)\" ",
|
|
130 |
tsdir,"/MOD35_",tile,"_ymonmiss.nc",sep=""))
|
| 128 |
131 |
|
| 129 |
132 |
## TODO: fix projection information so GDAL can read it correctly.
|
| 130 |
133 |
## clean up variables?
|
| ... | ... | |
| 134 |
137 |
system(paste(ncopath,"ncks -O ",tsdir,"/MOD35_",tile,"_ymonmean.nc ",tsdir,"/MOD35_",tile,"_ymon.nc",sep=""))
|
| 135 |
138 |
system(paste(ncopath,"ncks -A ",tsdir,"/MOD35_",tile,"_ymonstd.nc ",tsdir,"/MOD35_",tile,"_ymon.nc",sep=""))
|
| 136 |
139 |
system(paste(ncopath,"ncks -A ",tsdir,"/MOD35_",tile,"_ymoncld01.nc ",tsdir,"/MOD35_",tile,"_ymon.nc",sep=""))
|
|
140 |
system(paste(ncopath,"ncks -A ",tsdir,"/MOD35_",tile,"_ymonmiss.nc ",tsdir,"/MOD35_",tile,"_ymon.nc",sep=""))
|
| 137 |
141 |
|
| 138 |
142 |
## append sinusoidal grid from one of input files as CDO doesn't transfer all attributes
|
| 139 |
143 |
if(verbose) print("Clean up file (update attributes, flip latitudes, add grid description")
|
| 140 |
144 |
|
| 141 |
|
#system(paste(ncopath,"ncea -d time,0,1 -v sinusoidal ",list.files(outdir,full=T,pattern="[.]nc$")[1]," ",tsdir,"/sinusoidal.nc",sep=""))
|
| 142 |
|
#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=""))
|
| 143 |
|
|
| 144 |
145 |
## invert latitude so it plays nicely with gdal
|
| 145 |
146 |
system(paste(ncopath,"ncpdq -O -a -y ",tsdir,"/MOD35_",tile,"_ymon.nc ",outdir2,"/MOD35_",tile,".nc",sep=""))
|
| 146 |
147 |
|
|
148 |
## proj string taken from GDAL-written MODIS tile
|
|
149 |
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]]"
|
|
150 |
|
| 147 |
151 |
## update attributes
|
| 148 |
152 |
system(paste(ncopath,"ncatted ",
|
| 149 |
|
#" -a standard_parallel,sinusoidal,o,c,\"0\" ",
|
| 150 |
|
#" -a longitude_of_central_meridian,sinusoidal,o,c,\"0\" ",
|
| 151 |
|
#" -a latitude_of_central_meridian,sinusoidal,o,c,\"0\" ",
|
|
153 |
" -a false_easting,sinusoidal,o,d,0. ",
|
|
154 |
" -a false_northing,sinusoidal,o,d,0. ",
|
|
155 |
" -a longitude_of_central_meridian,sinusoidal,o,d,0. ",
|
|
156 |
" -a longitude_of_prime_meridian,sinusoidal,o,d,0. ",
|
|
157 |
" -a semi_major_axis,sinusoidal,o,d,6378137. ",
|
|
158 |
" -a inverse_flattening,sinusoidal,o,d,298.257223563 ",
|
|
159 |
" -a spatial_ref,sinusoidal,o,c,",projstring,
|
|
160 |
" -a GeoTransform,sinusoidal,o,c,\"-7783653.638366 926.6254331391661 0 1111950.519767 0 -926.6254331391667\" ",
|
| 152 |
161 |
" -a units,time,o,c,\"days since 2000-1-1 0:0:0\" ",
|
| 153 |
162 |
" -a title,global,o,c,\"MODIS Cloud Product (MOD35) Climatology\" ",
|
|
163 |
" -a tile,global,o,c,\"",tile,"\" ",
|
| 154 |
164 |
" -a institution,global,o,c,\"Yale University\" ",
|
| 155 |
|
" -a source,global,o,c,\"MODIS Cloud Product (MOD35)\" ",
|
|
165 |
" -a source,global,o,c,\"MODIS Cloud Product (MOD35) Collection 6\" ",
|
| 156 |
166 |
" -a comment,global,o,c,\"Compiled by Adam M. Wilson (adam.wilson@yale.edu)\" ",
|
| 157 |
167 |
outdir2,"/MOD35_",tile,".nc",sep=""))
|
| 158 |
168 |
|
| ... | ... | |
| 162 |
172 |
print(table(fdly$fyear))
|
| 163 |
173 |
print("Months:")
|
| 164 |
174 |
print(table(fdly$fmonth))
|
| 165 |
|
|
|
175 |
|
| 166 |
176 |
## quit R
|
| 167 |
177 |
q("no")
|
| 168 |
178 |
|
Removed some unsed code from MOD35_climatology. Still having problem of missing data stripes in February means.