Project

General

Profile

« Previous | Next » 

Revision 9a19743f

Added by Adam Wilson almost 11 years ago

addied initial MOD09 visualization script and updates to NDP script to use the new mod09 data

View differences:

climate/procedures/NDP-026D.R
7 7 
library(doMC)
8 8 
library(rasterVis)
9 9 
library(rgdal)
10 
library(reshape)
11 
library(hexbin)
10 12 
## register parallel processing
11 13 
registerDoMC(10)
12 14 

  		




  ......




  20
  22
  
    
st$lat=st$LAT/100


  		




  21
  23
  
    
st$lon=st$LON/100


  		




  22
  24
  
    
st$lon[st$lon>180]=st$lon[st$lon>180]-360


  		




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


  		




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


  		




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


  		




  23
  28
  
    

  		




  24
  29
  
    
## download data


  		




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


  		




  26
  
  
    
system("gunzip data/*.Z")


  		




  27
  31
  
    

  		




  28
  
  
    
## get monthly mean cloud amount MMCF


  		




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


  		




  30
  
  
    
#system("gunzip data/*.Z")


  		




  31
  
  
    
#f121=c(6,6,6,7,6,7,6,2) #format 121


  		




  32
  
  
    
#c121=c("StaID","NobD","AvgDy","NobN","AvgNt","NobDN","AvgDN","Acode")


  		




  33
  
  
    
#cld=do.call(rbind.data.frame,lapply(sprintf("%02d",1:12),function(m) {


  		




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


  		




  35
  
  
    
#  colnames(d)=c121


  		




  36
  
  
    
#  d$month=as.numeric(m)


  		




  37
  
  
    
#  return(d)}


  		




  38
  
  
    
#  ))


  		




  
  32
  
    
system("gunzip data/*.Z")


  		




  39
  33
  
    

  		




  40
  34
  
    
## define FWF widths


  		




  41
  35
  
    
f162=c(5,5,4,7,7,7,4) #format 162


  		




  ......




  78
  72
  
    
             cldsd=sd(x$Amt[x$Nobs>10]/100,na.rm=T))}))


  		




  79
  73
  
    
cldy[,c("lat","lon")]=st[match(cldy$StaID,st$StaID),c("lat","lon")]


  		




  80
  74
  
    

  		




  
  75
  
    
## add the MOD09 data to cld


  		




  
  76
  
    
#### Evaluate MOD35 Cloud data


  		




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


  		




  81
  78
  
    

  		




  82
  
  
    
#cldm=foreach(m=unique(cld$month),.combine='rbind')%:%


  		




  83
  
  
    
#  foreach(s=unique(cld$StaID),.combine="rbind") %dopar% {


  		




  84
  
  
    
#    x=cld[cld$month==m&cld$StaID==s,]


  		




  85
  
  
    
#    data.frame(


  		




  86
  
  
    
#               month=x$month[1],


  		




  87
  
  
    
#               StaID=x$StaID[1],


  		




  88
  
  
    
#               Amt=mean(x$Amt[x$Nobs>10],na.rm=T)/100)}


  		




  89
  
  
    
 


  		




  
  79
  
    
## overlay the data with 5km radius buffer


  		




  
  80
  
    
mod09st=extract(mod09,st,buffer=5000,fun=mean,na.rm=T,df=T)


  		




  
  81
  
    
mod09st$id=st$id


  		




  
  82
  
    
mod09stl=melt(mod09st[,-1],id.vars="id")


  		




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


  		




  
  84
  
    
## add it to cld


  		




  
  85
  
    
cld$mod09=mod09stl$value[match(paste(cld$StaID,cld$YR,cld$month),paste(mod09stl$id,mod09stl$year,as.numeric(mod09stl$month)))]


  		




  90
  86
  
    

  		




  91
  87
  
    
## write out the tables


  		




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


  		




  92
  89
  
    
write.csv(cldy,file="cldy.csv")


  		




  93
  90
  
    
write.csv(cldm,file="cldm.csv")


  		




  94
  91
  
    

  		




  95
  92
  
    
#########################################################################


  		




  96
  93
  
    
##################


  		




  97
  94
  
    
###


  		




  
  95
  
    
cld=read.csv("cld.csv")


  		




  98
  96
  
    
cldm=read.csv("cldm.csv")


  		




  99
  97
  
    
cldy=read.csv("cldy.csv")


  		




  
  98
  
    
st=read.csv("stations.csv")


  		




  100
  99
  
    

  		




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


  		




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


  		




  101
  102
  
    

  		




  102
  103
  
    
##make spatial object


  		




  103
  104
  
    
cldms=cldm


  		




  ......




  110
  111
  
    
projection(cldys)=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")


  		




  111
  112
  
    

  		




  112
  113
  
    
#### Evaluate MOD35 Cloud data


  		




  113
  
  
    
mod35c6=brick("~/acrobates/adamw/projects/MOD35C5/data/MOD35C6_2009_new.tif")


  		




  114
  
  
    
#projection(mod35c6)="+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs"


  		




  115
  
  
    

  		




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


  		




  116
  115
  
    

  		




  117
  
  
    
### use data from google earth engine


  		




  118
  
  
    
mod35c5=raster("../modis/mod09/global_2009/MOD35_2009.tif")


  		




  119
  
  
    
mod09=raster("../modis/mod09/global_2009/MOD09_2009.tif")


  		




  120
  116
  
    

  		




  121
  117
  
    
## LULC


  		




  122
  118
  
    
#system(paste("gdalwarp -r near -co \"COMPRESS=LZW\" -tr ",paste(res(mod09),collapse=" ",sep=""),


  		




  123
  119
  
    
#             "-tap -multi -t_srs \"",   projection(mod09),"\" /mnt/data/jetzlab/Data/environ/global/landcover/MODIS/MCD12Q1_IGBP_2005_v51.tif ../modis/mod12/MCD12Q1_IGBP_2005_v51.tif"))


  		




  124
  120
  
    
lulc=raster("../modis/mod12/MCD12Q1_IGBP_2005_v51.tif")


  		




  125
  
  
    
lulc=ratify(lulc)


  		




  
  121
  
    
#lulc=ratify(lulc)


  		




  126
  122
  
    
require(plotKML); data(worldgrids_pal)  #load IGBP palette


  		




  127
  123
  
    
IGBP=data.frame(ID=0:16,col=worldgrids_pal$IGBP[-c(18,19)],stringsAsFactors=F)


  		




  128
  124
  
    
IGBP$class=rownames(IGBP);rownames(IGBP)=1:nrow(IGBP)


  		




  129
  125
  
    
levels(lulc)=list(IGBP)


  		




  130
  
  
    
lulc=crop(lulc,mod09)


  		




  
  126
  
    
#lulc=crop(lulc,mod09)


  		




  131
  127
  
    

  		




  132
  128
  
    
n=100


  		




  133
  129
  
    
at=seq(0,100,length=n)


  		




  134
  130
  
    
colr=colorRampPalette(c("black","green","red"))


  		




  135
  131
  
    
cols=colr(n)


  		




  136
  132
  
    

  		




  137
  
  
    
#dif=mod35-mod09


  		




  138
  
  
    
#bwplot(dif~as.factor(lulc))


  		




  139
  133
  
    

  		




  140
  
  
    
#levelplot(mod35,col.regions=cols,at=at,margins=F,maxpixels=1e6)#,xlim=c(-100,-50),ylim=c(0,10))


  		




  141
  
  
    
#levelplot(lulc,att="class",col.regions=levels(lulc)[[1]]$col,margin=F,maxpixels=1e6)


  		




  
  134
  
    
hexbinplot(Amt/100~mod09,data=cld[cld$Nobs>100,])+


  		




  
  135
  
    
  layer(panel.abline(lm(y~x),col="blue"))+


  		




  
  136
  
    
  layer(panel.abline(0,1,col="red"))


  		




  142
  137
  
    

  		




  143
  
  
    
#cldys=spTransform(cldys,CRS(projection(mod35)))


  		




  
  138
  
    
xyplot(Amt/100~mod09,grpups="month",data=cld[cld$Nobs>75,],cex=.2,pch=16)+


  		




  
  139
  
    
  layer(panel.abline(lm(y~x),col="blue"))+


  		




  
  140
  
    
#  layer(panel.lines(x,predict(lm(y~x),type="prediction")))+


  		




  
  141
  
    
  layer(panel.abline(0,1,col="red"))


  		




  144
  142
  
    

  		




  145
  
  
    
#mod35v=foreach(m=unique(cldm$month),.combine="rbind") %do% {


  		




  146
  
  
    
#  dr=subset(mod35,subset=m);projection(dr)=projection(mod35)


  		




  147
  
  
    
#  dr2=subset(mod35sd,subset=m);projection(dr2)=projection(mod35)


  		




  148
  
  
    
#  ds=cldms[cldms$month==m,]


  		




  149
  
  
    
#  ds$mod35=unlist(extract(dr,ds,buffer=10,fun=mean,na.rm=T))


  		




  150
  
  
    
#  ds$mod35sd=extract(dr2,ds,buffer=10)


  		




  151
  
  
    
#  print(m)


  		




  152
  
  
    
#  return(ds@data[!is.na(ds$mod35),])}


  		




  
  143
  
    
xyplot(Amt/100~mod09|month,data=cld[cld$Nobs>75,],cex=.2,pch=16)+


  		




  
  144
  
    
  layer(panel.abline(lm(y~x),col="blue"))+


  		




  
  145
  
    
#  layer(panel.lines(x,predict(lm(y~x),type="prediction")))+


  		




  
  146
  
    
  layer(panel.abline(0,1,col="red"))


  		




  153
  147
  
    

  		




  154
  
  
    
y=2009


  		




  155
  
  
    
d=cldys[cldys$year==y,]


  		




  156
  
  
    

  		




  157
  
  
    
d$mod35c6_10=unlist(extract(mod35c6,d,buffer=10000,fun=mean,na.rm=T))


  		




  158
  
  
    
d$mod35c5_10=unlist(extract(mod35c5,d,buffer=10000,fun=mean,na.rm=T))


  		




  159
  
  
    
d$mod09_10=unlist(extract(mod09,d,buffer=10000,fun=mean,na.rm=T))


  		




  160
  
  
    
#d$dif=d$mod35_10-d$mod09_10


  		




  161
  
  
    
#d$dif2=d$mod35_10-d$cld


  		




  162
  148
  
    

  		




  163
  149
  
    
d$lulc=unlist(extract(lulc,d))


  		




  164
  150
  
    
d$lulc_10=unlist(extract(lulc,d,buffer=10000,fun=mode,na.rm=T))


  		




  ......




  166
  152
  
    

  		




  167
  153
  
    
save(d,file="annualsummary.Rdata")


  		




  168
  154
  
    

  		




  
  155
  
    

  		




  
  156
  
    

  		




  
  157
  
    
load("annualsummary.Rdata")


  		




  
  158
  
    

  		




  169
  159
  
    
## quick model to explore fit


  		




  170
  
  
    
plot(cld~mod35,groups=lulc,data=d)


  		




  171
  
  
    
summary(lm(cld~mod35+as.factor(lulc),data=d))


  		




  172
  
  
    
summary(lm(cld~mod09_10,data=d))


  		




  
  160
  
    
xyplot(cld~mod35c5_10,groups=lulc,data=d@data)


  		




  
  161
  
    
summary(lm(cld~mod35c5_10+as.factor(lulc),data=d@data))


  		




  
  162
  
    
summary(lm(Amt~mod09,data=cld))


  		




  173
  163
  
    
summary(lm(cld~mod09_10+as.factor(lulc),data=d))


  		




  174
  164
  
    
summary(lm(cld~mod09_10+as.factor(lulc),data=d))


  		




  175
  165
  
    

  		












Also available in:  Unified diff