Project

General

Profile

« Previous | Next » 

Revision d6ddb3e2

Added by Adam Wilson almost 11 years ago

udpating validation to full globe

View differences:

climate/procedures/NDP-026D.R
62 62 
## overlay the data with 32km diameter (16km radius) buffer
63 63 
## buffer size from Dybbroe, et al. (2005) doi:10.1175/JAM-2189.1.
64 64 
buf=16000
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 
#mod09sta=lapply(1:nlayers(mod09),function(l) {print(l); extract(mod09[[l]],st,buffer=buf,fun=mean,na.rm=T,df=T)[,2]})
66 
bins=cut(1:nrow(st),100)
67 
mod09sta=lapply(levels(bins),function(lb) {
68 
  l=which(bins==lb)
69 
  td=extract(mod09,st[l,],buffer=buf,fun=mean,na.rm=T,df=T)
70 
  td$id=st$id[l]
71 
  print(lb)#as.vector(c(l,td[,1:4])))
72 
  write.table(td,"valid.csv",append=T,col.names=F,quote=F,sep=",",row.names=F)
73 
  td
74 
})#,mc.cores=3)
75 

  		




  
  76
  
    
#mod09sta=extract(mod09,st,buffer=buf,fun=mean,na.rm=T,df=T)


  		




  
  77
  
    
mod09st=read.csv("valid.csv",header=F)[,-c(1,2)]


  		




  
  78
  
    

  		




  
  79
  
    
#mod09st=do.call(rbind.data.frame,mod09sta)


  		




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


  		




  68
  
  
    
colnames(mod09st)=names(mod09)


  		




  69
  
  
    
mod09st$id=st$id


  		




  
  81
  
    
colnames(mod09st)=c(names(mod09),"id")


  		




  
  82
  
    
#mod09st$id=st$id


  		




  70
  83
  
    
mod09stl=melt(mod09st,id.vars="id")


  		




  71
  84
  
    
mod09stl[,c("year","month")]=do.call(rbind,strsplit(sub("X","",mod09stl$variable),"[.]"))[,1:2]


  		




  72
  85
  
    
## add it to cld


  		




  ......




  85
  98
  
    
#lulc=crop(lulc,mod09)


  		




  86
  99
  
    
  Mode <- function(x) {


  		




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


  		




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


  		




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


  		




  89
  102
  
    
      }


  		




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


  		




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


  		




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


  		




  92
  105
  
    
## add it to cld


  		




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


  		




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


  		




  
  107
  
    
#cld$lulc=factor(as.integer(cld$lulc),labels=IGBP$class[sort(unique(cld$lulc))])


  		




  95
  108
  
    

  		




  96
  109
  
    
## update cld column names


  		




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


  		












Also available in:  Unified diff