Project

General

Profile

« Previous | Next » 

Revision 8009a075

Added by Adam Wilson almost 11 years ago

Updating validation to use new cloud data

View differences:

climate/procedures/MOD09_CloudFigures.R
41 41 
projection(cldys)=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
42 42 

  		




  43
  43
  
    
#### Evaluate MOD35 Cloud data


  		




  44
  
  
    
mod09=brick("data/mod09.nc")


  		




  45
  
  
    
mod09c=brick("data/mod09_clim_mean.nc",varname="CF");names(mod09c)=month.name


  		




  46
  
  
    
mod09a=brick("data/mod09_clim_mac.nc",varname="CF_annual")#;names(mod09c)=month.name


  		




  47
  
  
    

  		




  48
  
  
    
mod09min=raster("data/mod09_metrics.nc",varname="CFmin")


  		




  49
  
  
    
mod09max=raster("data/mod09_metrics.nc",varname="CFmax")


  		




  50
  
  
    
mod09sd=raster("data/mod09_metrics.nc",varname="CFsd")


  		




  51
  
  
    
mod09mean=raster("data/mod09_clim_mac.nc")


  		




  52
  
  
    

  		




  53
  
  
    
names(mod09d)=c("Mean","Minimum","Maximum","Standard Deviation")


  		




  
  44
  
    
mod09=brick("data/cloud_daily.nc")


  		




  
  45
  
    
mod09c=brick("data/cloud_ymonmean.nc",varname="CF");names(mod09c)=month.name


  		




  
  46
  
    
mod09a=brick("data/cloud_mean.nc",varname="CF_annual")#;names(mod09c)=month.name


  		




  
  47
  
    

  		




  
  48
  
    
mod09min=raster("data/cloud_min.nc",varname="CFmin")


  		




  
  49
  
    
mod09max=raster("data/cloud_max.nc",varname="CFmax")


  		




  
  50
  
    
mod09sd=raster("data/cloud_std.nc",varname="CFsd")


  		




  
  51
  
    
#mod09mean=raster("data/mod09_clim_mac.nc")


  		




  
  52
  
    
#names(mod09d)=c("Mean","Minimum","Maximum","Standard Deviation")


  		




  54
  53
  
    

  		




  55
  54
  
    
#plot(mod09a,layers=1,margin=F,maxpixels=100)


  		




  56
  55
  
    

  		












  

    
      climate/procedures/NDP-026D.R
    
  





  6
  6
  
    
library(doMC)


  		




  7
  7
  
    
library(rasterVis)


  		




  8
  8
  
    
library(rgdal)


  		




  9
  
  
    

  		




  
  9
  
    
library(reshape)


  		




  10
  10
  
    

  		




  11
  11
  
    

  		




  12
  12
  
    
## Data available here http://cdiac.ornl.gov/epubs/ndp/ndp026d/ndp026d.html


  		




  ......




  22
  22
  
    
colnames(st)=c("id","elev","lat","lon")


  		




  23
  23
  
    
write.csv(st,"stations.csv",row.names=F)


  		




  24
  24
  
    
coordinates(st)=c("lon","lat")


  		




  
  25
  
    
projection(st)="+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"


  		




  
  26
  
    

  		




  25
  27
  
    
## download data


  		




  26
  28
  
    
system("wget -N -nd ftp://cdiac.ornl.gov/pub/ndp026d/cat67_78/* -A '.tc.Z' -P data/")


  		




  27
  29
  
    

  		




  ......




  33
  35
  
    

  		




  34
  36
  
    
## use monthly timeseries


  		




  35
  37
  
    
cld=do.call(rbind.data.frame,mclapply(sprintf("%02d",1:12),function(m) {


  		




  36
  
  
    
  d=read.fwf(list.files("data",pattern=paste("MNYDC.",m,".tc",sep=""),full=T),skip=1,widths=f162)


  		




  
  38
  
    
  d=read.fwf(list.files("data",pattern=paste("MNYDC.",m,".tc$",sep=""),full=T),skip=1,widths=f162)


  		




  37
  39
  
    
  colnames(d)=c162


  		




  38
  40
  
    
  d$month=as.numeric(m)


  		




  39
  41
  
    
  print(m)


  		




  ......




  53
  55
  
    

  		




  54
  56
  
    
## add the MOD09 data to cld


  		




  55
  57
  
    
#### Evaluate MOD35 Cloud data


  		




  56
  
  
    
mod09=brick("~/acrobates/adamw/projects/cloud/data/mod09.nc")


  		




  
  58
  
    
mod09=brick("~/acrobates/adamw/projects/cloud/data/cloud_ymonmean.nc")


  		




  57
  59
  
    

  		




  58
  60
  
    
## overlay the data with 32km diameter (16km radius) buffer


  		




  59
  61
  
    
## buffer size from Dybbroe, et al. (2005) doi:10.1175/JAM-2189.1.


  		




  60
  62
  
    
buf=16000


  		




  61
  63
  
    
bins=cut(1:nrow(st),100)


  		




  
  64
  
    
if(file.exists("valid.csv")) file.remove("valid.csv")


  		




  62
  65
  
    
mod09sta=lapply(levels(bins),function(lb) {


  		




  63
  66
  
    
  l=which(bins==lb)


  		




  64
  67
  
    
  td=extract(mod09,st[l,],buffer=buf,fun=mean,na.rm=T,df=T)


  		












  

    
      climate/procedures/ee_compile.R
    
  





  36
  36
  
    

  		




  37
  37
  
    
## drop some if not complete


  		




  38
  38
  
    
#df=df[df$month%in%1:9&df$year%in%c(2001:2012),]


  		




  39
  
  
    
rerun=T  # set to true to recalculate all dates even if file already exists


  		




  
  39
  
    
rerun=F  # set to true to recalculate all dates even if file already exists


  		




  40
  40
  
    

  		




  41
  41
  
    
## Loop over existing months to build composite netcdf files


  		




  42
  42
  
    
foreach(date=unique(df$date)) %dopar% {


  		




  ......




  106
  106
  
    

  		




  107
  107
  
    
### merge all the tiles to a single global composite


  		




  108
  108
  
    
#system(paste("ncdump -h ",list.files(tempdir,pattern="mod09.*.nc$",full=T)[10]))


  		




  109
  
  
    
system(paste("cdo -O  mergetime ",paste(list.files(tempdir,pattern="mod09.*.nc$",full=T),collapse=" ")," data/cloud_daily.nc"))


  		




  
  109
  
    
system(paste("cdo -O  mergetime -setctomiss,-32768 ",paste(list.files(tempdir,pattern="mod09.*.nc$",full=T),collapse=" ")," data/cloud_monthly.nc"))


  		




  110
  110
  
    

  		




  111
  111
  
    
#  Overall mean


  		




  112
  
  
    
system(paste("cdo -O  -chname,CF,CF_annual -timmean data/cloud_daily.nc  data/cloud_mean.nc"))


  		




  
  112
  
    
system(paste("cdo -O  -chname,CF,CF_annual -timmean data/cloud_monthly.nc  data/cloud_mean.nc"))


  		




  113
  113
  
    

  		




  114
  114
  
    
### generate the monthly mean and sd


  		




  115
  115
  
    
#system(paste("cdo -P 10 -O merge -ymonmean data/mod09.nc -chname,CF,CF_sd -ymonstd data/mod09.nc data/mod09_clim.nc"))


  		




  116
  
  
    
system(paste("cdo  -O -ymonmean data/cloud_daily.nc data/cloud_ymonmean.nc"))


  		




  117
  
  
    
system(paste("cdo  -O -chname,CF,CF_sd -ymonstd data/cloud_daily.nc data/cloud_ymonsd.nc"))


  		




  
  116
  
    
system(paste("cdo  -f nc4c -O -ymonmean data/cloud_monthly.nc data/cloud_ymonmean.nc"))


  		




  
  117
  
    
system(paste("cdo  -f nc4c -O -chname,CF,CF_sd -ymonstd data/cloud_monthly.nc data/cloud_ymonsd.nc"))


  		




  118
  118
  
    

  		




  119
  119
  
    
#if(!file.exists("data/mod09_metrics.nc")) {


  		




  120
  
  
    
#    system("cdo -chname,CF,CFmin -timmin data/mod09_clim_mean.nc data/mod09_min.nc")


  		




  121
  
  
    
#    system("cdo -chname,CF,CFmax -timmax data/mod09_clim_mean.nc data/mod09_max.nc")


  		




  122
  
  
    
#    system("cdo -chname,CF,CFsd -timstd data/mod09_clim_mean.nc data/mod09_std.nc")


  		




  123
  
  
    
#    system("cdo -f nc2 merge data/mod09_std.nc data/mod09_min.nc data/mod09_max.nc data/mod09_metrics.nc")


  		




  124
  
  
    
    system("cdo merge -chname,CF,CFmin -timmin data/cloud_clim_mean.nc -chname,CF,CFmax -timmax data/cloud_clim_mean.nc  -chname,CF,CFsd -timstd data/cloud_clim_mean.nc  data/cloud_metrics.nc")


  		




  
  120
  
    
    system("cdo -f nc4c -chname,CF,CFmin -timmin data/cloud_ymonmean.nc data/cloud_min.nc")


  		




  
  121
  
    
    system("cdo -f nc4c -chname,CF,CFmax -timmax data/cloud_ymonmean.nc data/cloud_max.nc")


  		




  
  122
  
    
    system("cdo -f nc4c -chname,CF,CFsd -timstd data/cloud_ymonmean.nc data/cloud_std.nc")


  		




  
  123
  
    
#    system("cdo -f nc2 merge data/mod09_std.nc data/mod09_min.nc data/cloud_max.nc data/cloud_metrics.nc")


  		




  
  124
  
    
#    system("cdo merge -chname,CF,CFmin -timmin data/cloud_ymonmean.nc -chname,CF,CFmax -timmax data/cloud_ymonmean.nc  -chname,CF,CFsd -timstd data/cloud_ymonmean.nc  data/cloud_metrics.nc")


  		




  125
  125
  
    
#}


  		




  126
  126
  
    

  		




  127
  127
  
    

  		




  ......




  163
  163
  
    

  		




  164
  164
  
    

  		




  165
  165
  
    
## read in monthly dataset


  		




  166
  
  
    
mod09=brick("data/mod09_clim_mean.nc",varname="CF")


  		




  
  166
  
    
mod09=brick("data/cloud_ymonmean.nc",varname="CF")


  		




  167
  167
  
    
plot(mod09[1])


  		




  168
  168
  
    

  		




  169
  169
  
    
mod09_seas=calc(mod09,seasconc,return.Pc=T,return.thetat=F,overwrite=T,filename="data/mod09_seas.nc",NAflag=255,datatype="INT1U")


  		












Also available in:  Unified diff