Project

General

Profile

« Previous | Next » 

Revision 9a19743f

Added by Adam Wilson about 11 years ago

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

View differences:

climate/procedures/MOD09_Visualize.R
1 
## script to visualize cloud frequency data
2 

  		




  
  3
  
    
setwd("~/acrobates/adamw/projects/cloud/")


  		




  
  4
  
    

  		




  
  5
  
    
library(rasterVis)


  		




  
  6
  
    

  		




  
  7
  
    

  		




  
  8
  
    
#### Evaluate MOD35 Cloud data


  		




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


  		




  
  10
  
    
NAvalue(mod09)=-1


  		




  
  11
  
    

  		




  
  12
  
    
cols=colorRampPalette(c("black","blue","red"))


  		




  
  13
  
    

  		




  
  14
  
    
r=mod09[[10]]


  		




  
  15
  
    
levelplot(r,col.regions=cols(100),at=seq(0,100,len=100),margin=F,maxpixels=1e6)


  		




  
  16
  
    

  		




  
  17
  
    

  		




  
  18
  
    

  		




  
  19
  
    
## climatologies


  		




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


  		




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


  		




  
  22
  
    

  		




  
  23
  
    

  		




  
  24
  
    

  		




  
  25
  
    

  		












  

    
      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
  
    

  		












  

    
      climate/procedures/ee_compile.R
    
  





  
  1
  
    
###  Script to compile the monthly cloud data from earth engine into a netcdf file for further processing


  		




  
  2
  
    

  		




  
  3
  
    
setwd("~/acrobates/adamw/projects/cloud")


  		




  
  4
  
    
library(raster)


  		




  
  5
  
    
library(doMC)


  		




  
  6
  
    

  		




  
  7
  
    
registerDoMC(10)


  		




  
  8
  
    

  		




  
  9
  
    
tempdir="tmp"


  		




  
  10
  
    
if(!file.exists(tempdir)) dir.create(tempdir)


  		




  
  11
  
    

  		




  
  12
  
    
##  Get list of available files


  		




  
  13
  
    
df=data.frame(path=list.files("data/mod09",pattern="*.tif$",full=T),stringsAsFactors=F)


  		




  
  14
  
    
df[,c("region","year","month")]=do.call(rbind,strsplit(basename(df$path),"_|[.]"))[,c(2,3,4)]


  		




  
  15
  
    
df$date=as.Date(paste(df$year,"_",df$month,"_15",sep=""),"%Y_%m_%d")


  		




  
  16
  
    

  		




  
  17
  
    
table(df$year,df$month)


  		




  
  18
  
    
## drop some if not complete


  		




  
  19
  
    
df=df[df$year<=2006,]


  		




  
  20
  
    

  		




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


  		




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


  		




  
  23
  
    
## get date


  		




  
  24
  
    
  print(date)


  		




  
  25
  
    
  ## Define output and check if it already exists


  		




  
  26
  
    
  ncfile=paste(tempdir,"/mod09_",date,".nc",sep="")


  		




  
  27
  
    
  if(file.exists(ncfile)) next


  		




  
  28
  
    
  ## merge regions to a new netcdf file


  		




  
  29
  
    
  system(paste("gdal_merge.py -o ",ncfile," -of netCDF -ot Byte -n 0 ",paste(df$path[df$date==date],collapse=" ")))


  		




  
  30
  
    
  system(paste("ncecat -O -u time ",ncfile," ",ncfile,sep=""))


  		




  
  31
  
    
## create temporary nc file with time information to append to MOD06 data


  		




  
  32
  
    
  cat(paste("


  		




  
  33
  
    
    netcdf time {


  		




  
  34
  
    
      dimensions:


  		




  
  35
  
    
        time = 1 ;


  		




  
  36
  
    
      variables:


  		




  
  37
  
    
        int time(time) ;


  		




  
  38
  
    
      time:units = \"days since 2000-01-01 00:00:00\" ;


  		




  
  39
  
    
      time:calendar = \"gregorian\";


  		




  
  40
  
    
      time:long_name = \"time of observation\"; 


  		




  
  41
  
    
    data:


  		




  
  42
  
    
      time=",as.integer(date-as.Date("2000-01-01")),";


  		




  
  43
  
    
    }"),file=paste(tempdir,"/",date,"_time.cdl",sep=""))


  		




  
  44
  
    
system(paste("ncgen -o ",tempdir,"/",date,"_time.nc ",tempdir,"/",date,"_time.cdl",sep=""))


  		




  
  45
  
    
system(paste("ncks -A ",tempdir,"/",date,"_time.nc ",ncfile,sep=""))


  		




  
  46
  
    
## add other attributes


  		




  
  47
  
    
  system(paste("ncrename -v Band1,CF ",ncfile,sep=""))


  		




  
  48
  
    
  system(paste("ncatted ",


  		




  
  49
  
    
" -a units,CF,o,c,\"Proportion Days Cloudy\" ",


  		




  
  50
  
    
" -a missing_value,CF,o,b,255 ",


  		




  
  51
  
    
" -a _FillValue,CF,d,, ", 


  		




  
  52
  
    
" -a valid_range,CF,o,b,\"0,100\" ",


  		




  
  53
  
    
" -a long_name,CF,o,c,\"Proportion cloudy days (%)\" ",


  		




  
  54
  
    
ncfile,sep=""))


  		




  
  55
  
    
## add the fillvalue attribute back (without changing the actual values)


  		




  
  56
  
    
system(paste("ncatted -a _FillValue,CF,o,b,255 ",ncfile,sep=""))


  		




  
  57
  
    

  		




  
  58
  
    
if(as.numeric(system(paste("cdo -s ntime ",ncfile),intern=T))<1) {


  		




  
  59
  
    
  print(paste(ncfile," has no time, deleting"))


  		




  
  60
  
    
  file.remove(ntime)


  		




  
  61
  
    
}


  		




  
  62
  
    
  print(paste(basename(ncfile)," Finished"))


  		




  
  63
  
    

  		




  
  64
  
    
}


  		




  
  65
  
    

  		




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


  		




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


  		




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


  		




  
  69
  
    

  		




  
  70
  
    

  		




  
  71
  
    
### generate the monthly mean and sd


  		




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


  		




  
  73
  
    

  		




  
  74
  
    

  		












Also available in:  Unified diff