Project

General

Profile

« Previous | Next » 

Revision c1352601

Added by Adam Wilson about 11 years ago

Updated MOD35_Climatology to use percentages instead of classes, but there are too many zeros in result

View differences:

climate/procedures/MOD35_Climatology.r
76 76 

  		




  77
  77
  
    
## merge all daily files to create a single file with all dates


  		




  78
  78
  
    
system(paste(ncopath,"ncrcat -O ",outdir,"/*nc ",outdir2,"/MOD35_",tile,"_daily.nc",sep=""))


  		




  79
  
  
    

  		




  
  79
  
    
 


  		




  80
  80
  
    
## Update attributes


  		




  81
  81
  
    
system(paste(ncopath,"ncatted ",


  		




  82
  82
  
    
" -a units,time,o,c,\"days since 2000-1-1 0:0:0\" ",


  		




  ......




  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 -mulc,1.0 ",outdir2,"/MOD35_",tile,"_daily.nc ",tsdir,"/MOD35_",tile,"_ymonmean.nc",sep=""),wait=T)


  		




  
  100
  
    
system(paste("cdo -O sorttimestamp -setyear,",myear," -setday,15 -ymonmean ",outdir2,"/MOD35_",tile,"_daily.nc ",tsdir,"/MOD35_",tile,"_ymonmean.nc",sep=""),wait=T)


  		




  101
  101
  
    

  		




  102
  102
  
    
## Monthly standard deviation


  		




  103
  103
  
    
if(verbose) print("Calculating the monthly SD")


  		




  104
  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 CLD,CLD_sd ",tsdir,"/MOD35_",tile,"_ymonstd.nc",sep=""))


  		




  
  105
  
    
system(paste(ncopath,"ncrename -v PClear,PClear_sd ",tsdir,"/MOD35_",tile,"_ymonstd.nc",sep=""))


  		




  106
  106
  
    
system(paste(ncopath,"ncatted ",


  		




  107
  
  
    
#" -a long_name,CER_sd,o,c,\"Cloud Particle Effective Radius (standard deviation of daily observations)\" ",


  		




  108
  
  
    
" -a long_name,CLD_sd,o,c,\"Cloud Mask (standard deviation of daily observations)\" ",


  		




  109
  
  
    
#" -a long_name,COT_sd,o,c,\"Cloud Optical Thickness (standard deviation of daily observations)\" ",


  		




  
  107
  
    
" -a long_name,PClear_sd,o,c,\"Standard Deviation of p(clear)\" ",


  		




  110
  108
  
    
tsdir,"/MOD35_",tile,"_ymonstd.nc",sep=""))


  		




  111
  109
  
    

  		




  112
  
  
    
## cld == 0


  		




  113
  
  
    
if(verbose) print("Calculating the proportion of cloudy days")


  		




  114
  
  
    
system(paste("cdo -O  sorttimestamp -setyear,",myear," -setday,15 -nint -mulc,100 -ymonmean -mulc,1.0 -eqc,0 -setctomiss,1 -selvar,CLD2 ",outdir2,"/MOD35_",tile,"_daily.nc ",tsdir,"/MOD35_",tile,"_ymoncld0.nc",sep=""))


  		




  115
  
  
    
system(paste(ncopath,"ncrename -v CLD2,CLD0 ",tsdir,"/MOD35_",tile,"_ymoncld0.nc",sep=""))


  		




  116
  
  
    
system(paste(ncopath,"ncatted ",


  		




  117
  
  
    
" -a long_name,CLD0,o,c,\"Proportion of Days with Cloud Mask == 0\" ",


  		




  118
  
  
    
" -a units,CLD0,o,c,\"Proportion\" ",


  		




  119
  
  
    
tsdir,"/MOD35_",tile,"_ymoncld0.nc",sep=""))


  		




  120
  
  
    

  		




  121
  
  
    
## cld == 0|1


  		




  122
  
  
    
if(verbose) print("Calculating the proportion of cloudy days")


  		




  123
  
  
    
system(paste("cdo -O  sorttimestamp -setyear,",myear," -setday,15 -nint -mulc,100 -ymonmean  -mulc,1.0 -lec,1 -selvar,CLD2 ",outdir2,"/MOD35_",tile,"_daily.nc ",tsdir,"/MOD35_",tile,"_ymoncld01.nc",sep=""))


  		




  124
  
  
    
system(paste(ncopath,"ncrename -v CLD2,CLD01 ",tsdir,"/MOD35_",tile,"_ymoncld01.nc",sep=""))


  		




  
  110
  
    
## frequency of cloud days p(clear<90%)  


  		




  
  111
  
    
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=""))


  		




  
  113
  
    
system(paste(ncopath,"ncrename -v PClear,CF ",tsdir,"/MOD35_",tile,"_ymoncld01.nc",sep=""))


  		




  125
  114
  
    
system(paste(ncopath,"ncatted ",


  		




  126
  
  
    
" -a long_name,CLD01,o,c,\"Proportion of Days with Cloud Mask == 0|1\" ",


  		




  127
  
  
    
" -a units,CLD01,o,c,\"Proportion\" ",


  		




  
  115
  
    
" -a long_name,CF,o,c,\"Proportion of Days with probability of clear < 90%\" ",


  		




  
  116
  
    
" -a units,CF,o,c,\"Proportion\" ",


  		




  128
  117
  
    
tsdir,"/MOD35_",tile,"_ymoncld01.nc",sep=""))


  		




  129
  118
  
    

  		




  130
  
  
    
## cld == 1


  		




  131
  
  
    
if(verbose) print("Calculating the proportion of uncertain days")


  		




  132
  
  
    
system(paste("cdo -O  sorttimestamp -setyear,",myear," -setday,15 -nint -mulc,100 -ymonmean  -mulc,1.0 -eqc,1 -selvar,CLD2 ",outdir2,"/MOD35_",tile,"_daily.nc ",tsdir,"/MOD35_",tile,"_ymoncld1.nc",sep=""))


  		




  133
  
  
    
system(paste(ncopath,"ncrename -v CLD2,CLD1 ",tsdir,"/MOD35_",tile,"_ymoncld1.nc",sep=""))


  		




  134
  
  
    
system(paste(ncopath,"ncatted ",


  		




  135
  
  
    
" -a long_name,CLD1,o,c,\"Proportion of Days with Cloud Mask == 1 (uncertain)\" ",


  		




  136
  
  
    
" -a units,CLD1,o,c,\"Proportion\" ",


  		




  137
  
  
    
tsdir,"/MOD35_",tile,"_ymoncld1.nc",sep=""))


  		




  138
  
  
    

  		




  139
  
  
    

  		




  140
  
  
    
## cld >= 2 (setting cld==01 to missing because 'uncertain')


  		




  141
  
  
    
if(verbose) print("Calculating the proportion of clear days")


  		




  142
  
  
    
system(paste("cdo -O  sorttimestamp -setyear,",myear," -setday,15 -nint -mulc,100 -ymonmean  -mulc,1.0 -gtc,1 -setctomiss,1 -selvar,CLD2 ",outdir2,"/MOD35_",tile,"_daily.nc ",tsdir,"/MOD35_",tile,"_ymoncld2.nc",sep=""))


  		




  143
  
  
    
system(paste(ncopath,"ncrename -v CLD2,CLD23 ",tsdir,"/MOD35_",tile,"_ymoncld2.nc",sep=""))


  		




  144
  
  
    
system(paste(ncopath,"ncatted ",


  		




  145
  
  
    
" -a long_name,CLD23,o,c,\"Proportion of Days with Cloud Mask >= 2 (Probably Clear or Certainly Clear)\" ",


  		




  146
  
  
    
" -a units,CLD23,o,c,\"Proportion\" ",


  		




  147
  
  
    
tsdir,"/MOD35_",tile,"_ymoncld2.nc",sep=""))


  		




  148
  
  
    

  		




  149
  
  
    
## cld >= 1


  		




  150
  
  
    
if(verbose) print("Calculating the proportion of clear days")


  		




  151
  
  
    
system(paste("cdo -O  sorttimestamp -setyear,",myear," -setday,15 -nint -mulc,100 -ymonmean  -mulc,1.0 -gec,1 -selvar,CLD2 ",outdir2,"/MOD35_",tile,"_daily.nc ",tsdir,"/MOD35_",tile,"_ymoncld13.nc",sep=""))


  		




  152
  
  
    
system(paste(ncopath,"ncrename -v CLD2,CLD13 ",tsdir,"/MOD35_",tile,"_ymoncld13.nc",sep=""))


  		




  153
  
  
    
system(paste(ncopath,"ncatted ",


  		




  154
  
  
    
" -a long_name,CLD13,o,c,\"Proportion of Days with Cloud Mask >= 1\" ",


  		




  155
  
  
    
" -a units,CLD13,o,c,\"Proportion\" ",


  		




  156
  
  
    
tsdir,"/MOD35_",tile,"_ymoncld13.nc",sep=""))


  		




  157
  
  
    

  		




  158
  119
  
    
## number of observations


  		




  159
  
  
    
if(verbose) print("Calculating the number of missing variables")


  		




  160
  
  
    
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=""))


  		




  161
  
  
    
system(paste(ncopath,"ncrename -v CLD,CLD_pmiss ",tsdir,"/MOD35_",tile,"_ymonmiss.nc",sep=""))


  		




  162
  
  
    
system(paste(ncopath,"ncatted ",


  		




  163
  
  
    
" -a long_name,CLD_pmiss,o,c,\"Proportion of Days with missing data for CLD\" ",


  		




  164
  
  
    
" -a scale_factor,CLD_pmiss,o,d,0.01 ",


  		




  165
  
  
    
" -a units,CLD_pmiss,o,c,\"Proportion\" ",


  		




  166
  
  
    
tsdir,"/MOD35_",tile,"_ymonmiss.nc",sep=""))


  		




  
  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=""))


  		




  167
  128
  
    

  		




  168
  129
  
    
## TODO: fix projection information so GDAL can read it correctly.


  		




  169
  130
  
    
## clean up variables?


  		




  ......




  172
  133
  
    
if(verbose) print("Append all monthly climatologies into a single file")


  		




  173
  134
  
    
system(paste(ncopath,"ncks -O ",tsdir,"/MOD35_",tile,"_ymonmean.nc  ",tsdir,"/MOD35_",tile,"_ymon.nc",sep=""))


  		




  174
  135
  
    
system(paste(ncopath,"ncks -A ",tsdir,"/MOD35_",tile,"_ymonstd.nc  ",tsdir,"/MOD35_",tile,"_ymon.nc",sep=""))


  		




  175
  
  
    
system(paste(ncopath,"ncks -A ",tsdir,"/MOD35_",tile,"_ymoncld0.nc  ",tsdir,"/MOD35_",tile,"_ymon.nc",sep=""))


  		




  176
  136
  
    
system(paste(ncopath,"ncks -A ",tsdir,"/MOD35_",tile,"_ymoncld01.nc  ",tsdir,"/MOD35_",tile,"_ymon.nc",sep=""))


  		




  177
  
  
    
system(paste(ncopath,"ncks -A ",tsdir,"/MOD35_",tile,"_ymoncld1.nc  ",tsdir,"/MOD35_",tile,"_ymon.nc",sep=""))


  		




  178
  
  
    
system(paste(ncopath,"ncks -A ",tsdir,"/MOD35_",tile,"_ymoncld2.nc  ",tsdir,"/MOD35_",tile,"_ymon.nc",sep=""))


  		




  179
  
  
    
system(paste(ncopath,"ncks -A ",tsdir,"/MOD35_",tile,"_ymonmiss.nc  ",tsdir,"/MOD35_",tile,"_ymon.nc",sep=""))


  		




  180
  137
  
    

  		




  181
  138
  
    
## append sinusoidal grid from one of input files as CDO doesn't transfer all attributes


  		




  182
  139
  
    
if(verbose) print("Clean up file (update attributes, flip latitudes, add grid description")


  		












  

    
      climate/procedures/MOD35_L2_process.r
    
  





  182
  182
  
    
## Process the gridded files to align exactly with MODLAND tile and produce a daily summary of multiple swaths


  		




  183
  183
  
    
  


  		




  184
  184
  
    
## Identify output file


  		




  185
  
  
    
  ncfile=paste(outdir,"/MOD35_",tile,"_",date,".nc",sep="")  #this is the 'final' daily output file


  		




  
  185
  
    
ncfile=paste(outdir,"/MOD35_",tile,"_",date,".nc",sep="")  #this is the 'final' daily output file


  		




  186
  186
  
    

  		




  187
  187
  
    
## function to convert binary to decimal to assist in identifying correct values


  		




  188
  188
  
    
## this is helpful when defining QA handling below, but isn't used in processing


  		




  ......




  243
  243
  
    
                CM_cloud_",i," =  if((CM1_",i," / 2^0) % 2==1,(CM1_",i," / 2^1) % 2^2,-9999)


  		




  244
  244
  
    
                Pclear1_",i," = if(CM_cloud_",i,"==0,0,if(CM_cloud_",i,"==1,66,if(CM_cloud_",i,"==2,95,if(CM_cloud_",i,"==3,99,-9999))))


  		




  245
  245
  
    
                Pclear_",i," = if(QA_useful_",i,"==1,Pclear1_",i,",-9999)


  		




  246
  
  
    
                CM_path_",i," =   ((CM1_",i," / 2^6) % 2^2) 


  		




  247
  246
  
    
EOF",sep=""))


  		




  248
  247
  
    

  		




  
  248
  
    
     #                CM_path_",i," =   ((CM1_",i," / 2^6) % 2^2) 


  		




  
  249
  
    

  		




  249
  250
  
    
     execGRASS("r.null",map=paste("Pclear_",i,sep=""),setnull="-9999")


  		




  250
  251
  
    
     execGRASS("r.null",map=paste("CM_cloud_",i,sep=""),setnull="-9999")


  		




  251
  252
  
    
   


  		




  252
  253
  
    
    


  		




  253
  254
  
    
 } #end loop through sub daily files


  		




  254
  255
  
    

  		




  255
  
  
    
#### Now generate daily averages (or maximum in case of cloud flag)


  		




  256
  
  
    
  


  		




  
  256
  
    
#### Now generate daily minimum p(clear)


  		




  257
  257
  
    
  system(paste("r.mapcalc <<EOF


  		




  258
  258
  
    
         Pclear_daily=int((min(",paste("if(isnull(Pclear_",1:nfs,"),9999,Pclear_",1:nfs,")",sep="",collapse=","),"))) 


  		




  259
  259
  
    
EOF",sep=""))


  		




  260
  260
  
    

  		




  261
  
  
    
#         CLD_daily=int((min(",paste("if(isnull(CM_cloud_",1:nfs,"),9999,CM_cloud_",1:nfs,")",sep="",collapse=","),"))) 


  		




  262
  261
  
    

  		




  263
  262
  
    
## reset null values


  		




  264
  
  
    
#execGRASS("r.null",map="CLD_daily",setnull="9999")


  		




  265
  263
  
    
execGRASS("r.null",map="Pclear_daily",setnull="9999")


  		




  266
  264
  
    

  		




  267
  265
  
    

  		




  ......




  292
  290
  
    
system(paste(ncopath,"ncks -A ",tempdir(),"/time.nc ",ncfile,sep=""))


  		




  293
  291
  
    
## add other attributes


  		




  294
  292
  
    
  system(paste(ncopath,"ncrename -v Band1,PClear ",ncfile,sep=""))


  		




  295
  
  
    
  system(paste(ncopath,"ncatted -a scale_factor,PClear,o,d,1 -a units,PClear,o,c,\"Probability (%)\" -a missing_value,PClear,o,d,255 -a long_name,PClear,o,c,\"Probability of Clear Sky\" ",ncfile,sep=""))


  		




  
  293
  
    
  system(paste(ncopath,"ncatted -a scale_factor,PClear,o,d,1 -a units,PClear,o,c,\"Probability (%)\" -a missing_value,PClear,o,d,255 -a _FillValue,PClear,o,d,255 -a long_name,PClear,o,c,\"Probability of Clear Sky\" ",ncfile,sep=""))


  		




  296
  294
  
    

  		




  297
  295
  
    
                                        #  system(paste(ncopath,"ncatted -a sourcecode,global,o,c,",script," ",ncfile,sep=""))


  		




  298
  296
  
    
   


  		












  

    
      climate/procedures/Pleiades_MOD35.R
    
  





  1
  1
  
    
#### Script to facilitate processing of MOD06 data


  		




  2
  2
  
    
  


  		




  3
  
  
    
  setwd("/nobackupp1/awilso10/mod35")


  		




  
  3
  
    
setwd("/nobackupp1/awilso10/mod35")


  		




  4
  4
  
    

  		




  5
  5
  
    
library(rgdal)


  		




  6
  6
  
    
library(raster)


  		












Also available in:  Unified diff