Project

General

Profile

« Previous | Next » 

Revision a6d44f2d

Added by Adam Wilson over 11 years ago

Added script to process NDP-026D station cloud climatologies

View differences:

climate/procedures/MOD35_Explore.r
1 
## explore the MOD35 data downloaded and gridded by the DAAC
2 
setwd("~/acrobates/projects/interp/data/modis/mod35")
3 

  		




  
  4
  
    
library(raster)


  		




  
  5
  
    
library(rgdal)


  		




  
  6
  
    

  		




  
  7
  
    
f=list.files(pattern="Cloud_Mask_1_1")


  		




  
  8
  
    

  		




  
  9
  
    
GDALinfo(f[1])


  		




  
  10
  
    

  		




  
  11
  
    
## get tile


  		




  
  12
  
    
tile=raster("~/acrobates/projects/interp/data/modis/mod06/summary/MOD06_h09v04.nc",varname="CER")


  		




  
  13
  
    
h11v08=extent(tile)


  		




  
  14
  
    

  		




  
  15
  
    
r=raster(f[1])


  		




  
  16
  
    
extent(r)


  		




  
  17
  
    

  		




  
  18
  
    

  		




  
  19
  
    
st=lapply(f[1:10],raster)


  		




  
  20
  
    
str=lapply(2:length(st),function(i) union(extent(st[[i-1]]),extent(st[[i]])))[[length(st)-1]]


  		




  
  21
  
    
str=union(extent(h11v08),str)


  		




  
  22
  
    

  		




  
  23
  
    
b1=brick(lapply(st,function(stt) {


  		




  
  24
  
    
  x=crop(alignExtent(stt,str),h11v08)


  		




  
  25
  
    
  return(x)


  		




  
  26
  
    
}))


  		




  
  27
  
    

  		




  
  28
  
    

  		




  
  29
  
    

  		




  
  30
  
    
c=brick(f[1:10])


  		












  

    
      climate/procedures/NDP-026D.R
    
  





  6
  6
  
    

  		




  7
  7
  
    

  		




  8
  8
  
    
## Get station locations


  		




  9
  
  
    
system("wget -nd http://cdiac.ornl.gov/ftp/ndp026d/cat01/01_STID -P data/")


  		




  
  9
  
    
system("wget -N -nd http://cdiac.ornl.gov/ftp/ndp026d/cat01/01_STID -P data/")


  		




  10
  10
  
    
st=read.table("data/01_STID",skip=1)


  		




  11
  11
  
    
colnames(st)=c("StaID","LAT","LON","ELEV","ny1","fy1","ly1","ny7","fy7","ly7","SDC","b5c")


  		




  12
  12
  
    
st$lat=st$LAT/100


  		




  ......




  35
  35
  
    

  		




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


  		




  37
  37
  
    

  		




  38
  
  
    
## calculate means


  		




  39
  
  
    
cldm=by(cld,list(as.factor(cld$month),as.factor(cld$StaID)),function(x){


  		




  40
  
  
    
  data.frame(Amt=mean(x$Amt[x$Nobs>20],na.rm=T))})


  		




  41
  
  
    

  		




  42
  
  
    

  		




  43
  
  
    
## add color val


  		




  
  38
  
    
## drop missing values


  		




  44
  39
  
    
cld$Amt[cld$Amt<0]=NA


  		




  45
  40
  
    
cld$Fq[cld$Fq<0]=NA


  		




  46
  41
  
    
cld$AWP[cld$AWP<0]=NA


  		




  47
  42
  
    
cld$NC[cld$NC<0]=NA


  		




  48
  43
  
    

  		




  49
  
  
    
cld$col=cut(cld$Amt/100,quantile(cld$Amt/100,seq(0,1,len=5),na.rm=T))


  		




  
  44
  
    
## calculate means


  		




  
  45
  
    
cldm=do.call(rbind.data.frame,by(cld,list(month=as.factor(cld$month),StaID=as.factor(cld$StaID)),function(x){


  		




  
  46
  
    
  data.frame(month=x$month[1],StaID=x$StaID[1],Amt=mean(x$Amt[x$Nobs>20],na.rm=T))}))


  		




  
  47
  
    
cldm[,c("lat","lon")]=st[match(st$StaID,cldm$StaID),c("lat","lon")]


  		




  
  48
  
    

  		




  
  49
  
    

  		




  
  50
  
    

  		




  
  51
  
    
## write out the table


  		




  
  52
  
    
write.csv(cldm,file="cldm.csv")


  		




  
  53
  
    

  		




  
  54
  
    

  		




  
  55
  
    
##################


  		




  
  56
  
    
###


  		




  
  57
  
    
cldm=read.csv("cldm.csv")


  		




  
  58
  
    

  		




  
  59
  
    
## add a color key


  		




  
  60
  
    
cldm$col=cut(cldm$Amt/100,quantile(cldm$Amt/100,seq(0,1,len=5),na.rm=T))


  		




  50
  61
  
    

  		




  51
  62
  
    
library(lattice)


  		




  52
  
  
    
xyplot(lat~lon|as.factor(YR)+as.factor(month),groups=col,data=cld,pch=16,cex=.2,auto.key=T)


  		




  
  63
  
    
xyplot(lat~lon|+month,groups=col,data=cldm,pch=16,cex=.2,auto.key=T)


  		




  
  64
  
    

  		












Also available in:  Unified diff