Project

General

Profile

« Previous | Next » 

Revision d86b0a4a

Added by Adam Wilson almost 11 years ago

updates to MOD09 visualizations

View differences:

climate/procedures/MOD09_Visualize.R
14 14 

  		




  15
  15
  
    

  		




  16
  16
  
    
#### Evaluate MOD35 Cloud data


  		




  17
  
  
    
mod09=brick("~/acrobates/adamw/projects/cloud/data/mod09.nc",varname="CF")


  		




  
  17
  
    
mc=brick("~/acrobates/adamw/projects/cloud/data/mod09.nc",varname="CF")


  		




  18
  18
  
    
NAvalue(mod09)=-1


  		




  19
  19
  
    

  		




  20
  20
  
    
cols=colorRampPalette(c("#000000","#00FF00","#FF0000"))#"black","blue","red"))


  		




  21
  21
  
    
for(i in 1:156){


  		




  22
  22
  
    
png(paste("output/mod09_fullanimation_",i,".png",sep=""),width=2000,height=1000)


  		




  23
  23
  
    
  print(i)


  		




  24
  
  
    
  r=mod09[[i]]


  		




  
  24
  
    
  r=mm[[i]]


  		




  25
  25
  
    
  print(levelplot(r,col.regions=cols(100),at=seq(0,100,len=100),margin=F,maxpixels=1e6,ylim=c(-60,70),main=paste(names(mod09)[i])))+


  		




  26
  26
  
    
    layer(sp.lines(coast))


  		




  27
  27
  
    
dev.off()


  		




  28
  28
  
    
}


  		




  29
  29
  
    

  		




  30
  30
  
    
#### Evaluate MOD35 Cloud data


  		




  31
  
  
    
mod09=brick("~/acrobates/adamw/projects/cloud/data/mod09_clim.nc",varname="CF")


  		




  32
  
  
    
NAvalue(mod09)=-1


  		




  
  31
  
    
mmc=brick("~/acrobates/adamw/projects/cloud/data/mod09_clim_mean.nc",varname="CF")


  		




  
  32
  
    
names(mmc)=month.name


  		




  
  33
  
    
NAvalue(mmc)=-1


  		




  33
  34
  
    

  		




  34
  35
  
    
cols=colorRampPalette(c("#000000","#00FF00","#FF0000"))#"black","blue","red"))


  		




  35
  36
  
    
for(i in 1:12){


  		




  36
  37
  
    
png(paste("output/mod09_animation_",i,".png",sep=""),width=2000,height=1000)


  		




  37
  38
  
    
  print(i)


  		




  38
  
  
    
  r=mod09[[i]]


  		




  39
  
  
    
  print(levelplot(r,col.regions=cols(100),at=seq(0,100,len=100),margin=F,maxpixels=1e6,ylim=c(-60,70),main=paste(month.name[i])))+


  		




  
  39
  
    
  r=mmc[[i]]


  		




  
  40
  
    
  print(levelplot(r,col.regions=cols(100),at=seq(0,100,len=100),margin=F,maxpixels=1e7,ylim=c(-60,70),


  		




  
  41
  
    
                  main=paste(month.name[i]),cex.main=3,scales=list(draw=F),cuts=99))+


  		




  40
  42
  
    
    layer(sp.lines(coast))


  		




  41
  43
  
    
dev.off()


  		




  42
  44
  
    
}


  		




  43
  45
  
    

  		




  44
  46
  
    

  		




  45
  47
  
    
## climatologies


  		




  46
  
  
    
mod09a=brick("~/acrobates/adamw/projects/cloud/data/mod09_clim2.nc",varname="CF")


  		




  
  48
  
    
mac=brick("~/acrobates/adamw/projects/cloud/data/mod09_clim_mac.nc",varname="CF_annual")


  		




  47
  49
  
    

  		




  48
  50
  
    
pdf("output/mod09_climatology.pdf",width=11,height=8.5)


  		




  49
  
  
    
levelplot(mod09c,col.regions=cols(100),at=seq(0,100,len=100),margin=F,max.pixels=1e7)


  		




  
  51
  
    
levelplot(mac,col.regions=grey(seq(0,1,len=100)),cuts=99,margin=F,max.pixels=1e6)+


  		




  
  52
  
    
    layer(sp.lines(coast,lwd=.5,col="black"))


  		




  50
  53
  
    
dev.off()


  		




  51
  54
  
    

  		




  52
  55
  
    

  		




  53
  
  
    
## Compare with worldclim


  		




  
  56
  
    
## Compare with worldclim and NPP


  		




  54
  57
  
    
wc=stack(as.list(paste("/mnt/data/jetzlab/Data/environ/global/worldclim/prec_",1:12,".bil",sep="")))


  		




  55
  
  
    
reg=extent(c(-83,-45,-5,13))


  		




  56
  
  
    
reg2=extent(c(-81,-70,-4,10))


  		




  57
  
  
    

  		




  
  58
  
    
wc_map=stack(as.list(paste("/mnt/data/jetzlab/Data/environ/global/worldclim/bio_12.bil",sep="")))


  		




  
  59
  
    
npp=raster("/mnt/data/jetzlab/Data/environ/global/MODIS/MOD17A3/MOD17A3_Science_NPP_mean_00_12.tif",sep="")


  		




  58
  60
  
    

  		




  59
  
  
    
wc_map=mean(crop(wc,reg))


  		




  60
  61
  
    

  		




  61
  
  
    
twc_12=crop(wc[[12]],reg)


  		




  62
  
  
    
mod09_12=crop(mod09[[12]],reg)


  		




  
  62
  
    
pdf("output/mod09_worldclim.pdf",width=11,height=8.5)


  		




  
  63
  
    
p1=levelplot(wc_map,col.regions=grey(seq(0,1,len=100)),cuts=99,margin=F,maxpixels=1e5,colorkey=list("top"))


  		




  
  64
  
    
p2=levelplot(mac,col.regions=grey(seq(0,1,len=100)),cuts=99,margin=F,maxpixels=1e5,colorkey=list("bottom"))


  		




  
  65
  
    
p3=c(p1,p2,x.same=T,y.same=T,merge.legends=T)


  		




  
  66
  
    
print(p3)


  		




  63
  67
  
    

  		




  64
  68
  
    

  		




  65
  
  
    
p1=levelplot(wc_map,col.regions=grey(seq(0,1,len=100)),cuts=99,margin=F,colorkey=list("top"))


  		




  66
  
  
    
p2=levelplot(mod09_12,col.regions=grey(seq(0,1,len=100)),cuts=99,margin=F,colorkey=list("bottom"))


  		




  
  69
  
    
regs=list(


  		




  
  70
  
    
  Hawaii=extent(c(-156.5,-154,18.75,20.5)),


  		




  
  71
  
    
  Boliva=extent(c(-71,-63,-20,-15)),


  		




  
  72
  
    
  Venezuela=extent(c(-69,-59,0,7)),


  		




  
  73
  
    
  reg=extent(c(-83,-45,-5,13)),


  		




  
  74
  
    
  reg2=extent(c(-81,-70,-4,10))


  		




  
  75
  
    
  )


  		




  67
  76
  
    

  		




  68
  
  
    
p3=c(p1,p2,x.same=T,y.same=T,merge.legends=F)


  		




  
  77
  
    
r=3  #set region


  		




  
  78
  
    
p_map=levelplot(crop(wc_map,regs[[r]]),col.regions=grey(seq(0,1,len=100)),cuts=99,margin=F,maxpixels=1e5,colorkey=list(space="bottom",height=.75),xlab="",ylab="",main=names(regs)[r],useRaster=T)


  		




  
  79
  
    
p_mac=levelplot(crop(mac,regs[[r]]),col.regions=grey(seq(0,1,len=100)),cuts=99,margin=F,maxpixels=1e5,colorkey=list(space="bottom",height=.75),useRaster=T)


  		




  
  80
  
    
p_npp=levelplot(crop(npp,regs[[r]]),col.regions=grey(seq(0,1,len=100)),cuts=99,margin=F,maxpixels=1e5,colorkey=list(space="bottom",height=.5),zscaleLog=T,useRaster=T)


  		




  
  81
  
    
p3=c("WorldClim Mean Annual Precip (mm)"=p_map,"MOD09 Cloud Frequency (%)"=p_mac,"NPP"=p_npp,x.same=T,y.same=T,merge.legends=T,layout=c(3,1))


  		




  69
  82
  
    
print(p3)


  		




  70
  83
  
    

  		




  71
  
  
    
plot(mod09_12)


  		




  
  84
  
    
dev.off()


  		




  72
  85
  
    

  		




  73
  86
  
    

  		




  74
  87
  
    
## reduced resolution


  		




  75
  
  
    
mod09_8km=aggregate(mod09_12,8)


  		




  76
  
  
    
mod09_1deg=aggregate(mod09_12,110)


  		




  
  88
  
    
mod09_8km=aggregate(mod09_mac,8)


  		




  
  89
  
    
mod09_1deg=aggregate(mod09_mac,110)


  		




  77
  90
  
    

  		




  78
  91
  
    
pdf("output/mod09_resolution.pdf",width=11,height=8.5)


  		




  79
  
  
    
p1=levelplot(mod09_12,col.regions=grey(seq(0,1,len=100)),cuts=99,margin=F,max.pixels=1e6)


  		




  80
  
  
    
p2=levelplot(mod09_8km,col.regions=grey(seq(0,1,len=100)),cuts=99,margin=F,max.pixels=1e6)


  		




  81
  
  
    
p3=levelplot(mod09_1deg,col.regions=grey(seq(0,1,len=100)),cuts=99,margin=F,max.pixels=1e6)


  		




  82
  
  
    
print(c(p1,p2,p3))


  		




  
  92
  
    
p1=levelplot(mod09_mac,col.regions=grey(seq(0,1,len=100)),cuts=99,margin=F,max.pixels=1e5)


  		




  
  93
  
    
p2=levelplot(mod09_8km,col.regions=grey(seq(0,1,len=100)),cuts=99,margin=F,max.pixels=1e5)


  		




  
  94
  
    
p3=levelplot(mod09_1deg,col.regions=grey(seq(0,1,len=100)),cuts=99,margin=F,max.pixels=1e5)


  		




  
  95
  
    
print(c(p1,p2,p3,x.same=T,y.same=T,merge.legends=F))


  		




  83
  96
  
    

  		




  84
  97
  
    

  		




  85
  
  
    
p1=levelplot(crop(mod09_12,reg2),col.regions=grey(seq(0,1,len=100)),cuts=99,margin=F,max.pixels=1e6)


  		




  86
  
  
    
p2=levelplot(crop(mod09_8km,reg2),col.regions=grey(seq(0,1,len=100)),cuts=99,margin=F,max.pixels=1e6)


  		




  87
  
  
    
p3=levelplot(crop(mod09_1deg,reg2),col.regions=grey(seq(0,1,len=100)),cuts=99,margin=F,max.pixels=1e6)


  		




  88
  
  
    
print(c(p1,p2,p3),x.same=T,y.same=T)


  		




  
  98
  
    
p1=levelplot(crop(mod09_mac,reg2),col.regions=grey(seq(0,1,len=100)),cuts=99,margin=F,max.pixels=1e5)


  		




  
  99
  
    
p2=levelplot(crop(mod09_8km,reg2),col.regions=grey(seq(0,1,len=100)),cuts=99,margin=F,max.pixels=1e5)


  		




  
  100
  
    
p3=levelplot(crop(mod09_1deg,reg2),col.regions=grey(seq(0,1,len=100)),cuts=99,margin=F,max.pixels=1e5)


  		




  
  101
  
    
print(c(p1,p2,p3,x.same=T,y.same=T,merge.legends=F))


  		




  89
  102
  
    

  		




  90
  
  
    
dev.off()


  		




  
  103
  
    
p1=levelplot(crop(mod09_mac,reg3),col.regions=grey(seq(0,1,len=100)),cuts=99,margin=F,max.pixels=1e5)


  		




  
  104
  
    
p2=levelplot(crop(mod09_8km,reg3),col.regions=grey(seq(0,1,len=100)),cuts=99,margin=F,max.pixels=1e5)


  		




  
  105
  
    
p3=levelplot(crop(mod09_1deg,reg3),col.regions=grey(seq(0,1,len=100)),cuts=99,margin=F,max.pixels=1e5)


  		




  
  106
  
    
print(c(p1,p2,p3,x.same=T,y.same=T,merge.legends=F))


  		




  91
  107
  
    

  		




  
  108
  
    
dev.off()


  		




  92
  109
  
    

  		












  

    
      climate/procedures/NDP-026D.R
    
  





  25
  25
  
    
st=st[,c("StaID","ELEV","lat","lon")]


  		




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


  		




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


  		




  28
  
  
    

  		




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


  		




  29
  29
  
    
## download data


  		




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


  		




  31
  31
  
    

  		




  ......




  45
  45
  
    
  ))


  		




  46
  46
  
    

  		




  47
  47
  
    
## add lat/lon


  		




  48
  
  
    
cld[,c("lat","lon")]=st[match(cld$StaID,st$StaID),c("lat","lon")]


  		




  
  48
  
    
cld[,c("lat","lon")]=st[match(cld$StaID,st$id),c("lat","lon")]


  		




  49
  49
  
    

  		




  50
  50
  
    
## drop missing values


  		




  
  51
  
    
cld=cld[,!grepl("Fq|AWP|NC",colnames(cld))]


  		




  51
  52
  
    
cld$Amt[cld$Amt<0]=NA


  		




  52
  
  
    
cld$Fq[cld$Fq<0]=NA


  		




  53
  
  
    
cld$AWP[cld$AWP<0]=NA


  		




  54
  
  
    
cld$NC[cld$NC<0]=NA


  		




  55
  
  
    
cld=cld[cld$Nobs>0,]


  		




  
  53
  
    
#cld$Fq[cld$Fq<0]=NA


  		




  
  54
  
    
#cld$AWP[cld$AWP<0]=NA


  		




  
  55
  
    
#cld$NC[cld$NC<0]=NA


  		




  
  56
  
    
#cld=cld[cld$Nobs>0,]


  		




  56
  57
  
    

  		




  57
  58
  
    
## add the MOD09 data to cld


  		




  58
  59
  
    
#### Evaluate MOD35 Cloud data


  		




  ......




  61
  62
  
    
## overlay the data with 32km diameter (16km radius) buffer


  		




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


  		




  63
  64
  
    
buf=16000


  		




  64
  
  
    
mod09st=do.call(cbind.data.frame,mclapply(1:nlayers(mod09),function(l) extract(mod09[[l]],st,buffer=buf,fun=mean,na.rm=T,df=T)[,2]))


  		




  
  65
  
    
mod09sta=lapply(1:nlayers(mod09),function(l) {print(l); extract(mod09[[l]],st,buffer=buf,fun=mean,na.rm=T,df=T)[,2]})


  		




  
  66
  
    
mod09st=do.call(cbind.data.frame,mod09sta)


  		




  65
  67
  
    
#mod09st=mod09st[,!is.na(colnames(mod09st))]


  		




  66
  68
  
    
colnames(mod09st)=names(mod09)


  		




  67
  69
  
    
mod09st$id=st$id


  		




  ......




  82
  84
  
    
levels(lulc)=list(IGBP)


  		




  83
  85
  
    
#lulc=crop(lulc,mod09)


  		




  84
  86
  
    
  Mode <- function(x) {


  		




  85
  
  
    
      ux <- unique(x)


  		




  86
  
  
    
        ux[which.max(tabulate(match(x, ux)))]


  		




  
  87
  
    
      ux <- na.omit(unique(x))


  		




  
  88
  
    
        ux[which.max(tabulate(match(x, ux)),na.rm=T)]


  		




  87
  89
  
    
      }


  		




  88
  
  
    
lulcst=extract(lulc,st,fun=Mode,buffer=buf,df=T)


  		




  
  90
  
    
lulcst=extract(lulc,st,fun=Mode,na.rm=T,buffer=buf,df=T)


  		




  89
  91
  
    
colnames(lulcst)=c("id","lulc")


  		




  90
  92
  
    
## add it to cld


  		




  91
  93
  
    
cld$lulc=lulcst$lulc[match(cld$StaID,lulcst$id)]


  		




  92
  
  
    
cld$lulc=factor(cld$lulc,labels=IGBP$class)


  		




  
  94
  
    
cld$lulc=factor(as.integer(cld$lulc),labels=IGBP$class)


  		




  93
  95
  
    

  		




  94
  96
  
    
## update cld column names


  		




  95
  97
  
    
colnames(cld)[grep("Amt",colnames(cld))]="cld"


  		




  ......




  135
  137
  
    
write.csv(cld,file="cld.csv",row.names=F)


  		




  136
  138
  
    
write.csv(cldy,file="cldy.csv",row.names=F)


  		




  137
  139
  
    
write.csv(cldm,file="cldm.csv",row.names=F)


  		




  138
  
  
    
write.csv(clda,file="clda.csv",row.names=F)


  		




  139
  
  
    

  		




  
  140
  
    
write.csv(clda,file="clda.csv",row.names=F


  		




  
  141
  
    
)


  		




  140
  142
  
    
#########################################################################


  		




  141
  143
  
    
##################


  		




  142
  144
  
    
###


  		












  

    
      climate/procedures/ee_compile.R
    
  





  76
  76
  
    

  		




  77
  77
  
    

  		




  78
  78
  
    
### generate the monthly mean and sd


  		




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


  		




  80
  
  
    
system(paste("cdo -O -ymonmean data/mod09.nc data/mod09_clim2.nc"))


  		




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


  		




  
  80
  
    
system(paste("cdo  -O -ymonmean data/mod09.nc data/mod09_clim_mean.nc"))


  		




  
  81
  
    
system(paste("cdo  -O -chname,CF,CF_sd -ymonstd data/mod09.nc data/mod09_clim_sd.nc"))


  		




  81
  82
  
    

  		




  82
  83
  
    
#  Overall mean


  		




  83
  
  
    
system(paste("cdo -O  -chname,CF,CF_annual -timmean data/mod09.nc  data/mod09_clim2.nc"))


  		




  
  84
  
    
system(paste("cdo -O  -chname,CF,CF_annual -timmean data/mod09.nc  data/mod09_clim_mac.nc"))


  		




  84
  85
  
    

  		




  85
  86
  
    
### Long term summaries


  		




  86
  87
  
    
seasconc <- function(x,return.Pc=T,return.thetat=F) {


  		




  ......




  114
  115
  
    

  		




  115
  116
  
    

  		




  116
  117
  
    
## read in monthly dataset


  		




  117
  
  
    
mod09=brick("data/mod09_clim.nc",varname="CF")


  		




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


  		




  118
  119
  
    
plot(mod09[1])


  		




  119
  120
  
    

  		




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


  		




  
  122
  
    
mod09_seas2=calc(mod09,seasconc,return.Pc=F,return.thetat=T,overwrite=T,filename="data/mod09_seas_theta.nc",datatype="INT1U")


  		




  121
  123
  
    

  		




  122
  124
  
    
plot(mod09_seas)


  		












Also available in:  Unified diff