Project

General

Profile

« Previous | Next » 

Revision ddd9a810

Added by Adam Wilson over 11 years ago

Updated script to process MOD35 Processing Path using HEG to interpolate sensor zenith observations. Also working on figures for MOD35 C5 vs C6 comparison in MOD35C5_Evaluation script

View differences:

climate/procedures/NDP-026D.R
5 5 
library(multicore)
6 6 
library(latticeExtra)
7 7 
library(doMC)
8 
library(raster)
8 
library(rasterVis)
9 9 
library(rgdal)
10 10 
## register parallel processing
11 11 
registerDoMC(20)
......
105 105 
coordinates(cldms)=c("lon","lat")
106 106 
projection(cldms)=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
107 107 

  		




  
  108
  
    
##make spatial object


  		




  
  109
  
    
cldys=cldy


  		




  
  110
  
    
coordinates(cldys)=c("lon","lat")


  		




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


  		




  
  112
  
    

  		




  108
  113
  
    
#### Evaluate MOD35 Cloud data


  		




  109
  114
  
    
mod35=brick("../modis/mod35/MOD35_h11v08.nc",varname="CLD01")


  		




  110
  115
  
    
mod35sd=brick("../modis/mod35/MOD35_h11v08.nc",varname="CLD_sd")


  		




  ......




  112
  117
  
    

  		




  113
  118
  
    

  		




  114
  119
  
    
### use data from google earth engine


  		




  115
  
  
    
mod35=brick("../modis/mod09/global_2009/")


  		




  
  120
  
    
mod35=raster("../modis/mod09/global_2009/MOD35_2009.tif")


  		




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


  		




  117
  122
  
    

  		




  
  123
  
    
## LULC


  		




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


  		




  
  125
  
    
             "-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"))


  		




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


  		




  
  127
  
    
lulc=ratify(lulc)


  		




  
  128
  
    
require(plotKML); data(worldgrids_pal)  #load IGBP palette


  		




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


  		




  
  130
  
    
IGBP$class=rownames(IGBP);rownames(IGBP)=1:nrow(IGBP)


  		




  
  131
  
    
levels(lulc)=list(IGBP)


  		




  
  132
  
    
lulc=crop(lulc,mod09)


  		




  
  133
  
    

  		




  118
  134
  
    
n=100


  		




  119
  135
  
    
at=seq(0,100,length=n)


  		




  120
  136
  
    
colr=colorRampPalette(c("black","green","red"))


  		




  121
  137
  
    
cols=colr(n)


  		




  122
  138
  
    

  		




  123
  
  
    
levelplot(mod09,col.regions=cols,at=at)


  		




  
  139
  
    
dif=mod35-mod09


  		




  
  140
  
    
bwplot(dif~as.factor(lulc))


  		




  124
  141
  
    

  		




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


  		




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


  		




  125
  144
  
    

  		




  126
  
  
    
cldms=spTransform(cldms,CRS(projection(mod35)))


  		




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


  		




  127
  146
  
    

  		




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


  		




  129
  148
  
    
  dr=subset(mod35,subset=m);projection(dr)=projection(mod35)


  		




  ......




  134
  153
  
    
  print(m)


  		




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


  		




  136
  155
  
    

  		




  
  156
  
    
y=2009


  		




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


  		




  
  158
  
    

  		




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


  		




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


  		




  
  161
  
    
d$dif=d$mod35_10-d$mod09_10


  		




  
  162
  
    
d$dif2=d$mod35_10-d$cld


  		




  
  163
  
    

  		




  
  164
  
    
d$lulc=unlist(extract(lulc,d))


  		




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


  		




  
  166
  
    
d$lulc=factor(d$lulc,labels=IGBP$class)


  		




  
  167
  
    

  		




  
  168
  
    
save(d,file="annualsummary.Rdata")


  		




  
  169
  
    

  		




  
  170
  
    
## quick model to explore fit


  		




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


  		




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


  		




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


  		




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


  		




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


  		




  
  176
  
    

  		




  
  177
  
    
### exploratory plots


  		




  
  178
  
    
xyplot(cld~mod09_10,groups=lulc,data=d@data,pch=16,cex=.5)+layer(panel.abline(0,1,col="red"))


  		




  
  179
  
    
xyplot(cld~mod09_10+mod35_10|as.factor(lulc),data=d@data,type=c("p","r"),pch=16,cex=.25,auto.key=T)+layer(panel.abline(0,1,col="green"))


  		




  
  180
  
    
xyplot(cld~mod35_10|as.factor(lulc),data=d@data,pch=16,cex=.5)+layer(panel.abline(0,1,col="red"))


  		




  
  181
  
    
xyplot(mod35_10~mod09_10|as.factor(lulc),data=d@data,pch=16,cex=.5)+layer(panel.abline(0,1,col="red"))


  		




  
  182
  
    

  		




  
  183
  
    
densityplot(stack(mod35,mod09))


  		




  
  184
  
    
boxplot(mod35,lulc)


  		




  
  185
  
    

  		




  
  186
  
    
bwplot(mod09~mod35|cut(y,5),data=stack(mod09,mod35))


  		




  
  187
  
    

  		




  137
  188
  
    
## month factors


  		




  138
  189
  
    
cldm$month2=factor(cldm$month,labels=month.name)


  		




  139
  190
  
    
## add a color key


  		












Also available in:  Unified diff