Revision f2b2e7d5
Added by Adam M. Wilson about 12 years ago
| climate/procedures/MOD06_L2_data_summary.r | ||
|---|---|---|
| 39 | 39 |
########################## |
| 40 | 40 |
#### explore the data |
| 41 | 41 |
|
| 42 |
## load data |
|
| 43 | 42 |
months=seq(as.Date("2000-01-15"),as.Date("2000-12-15"),by="month")
|
| 43 |
|
|
| 44 |
|
|
| 45 |
## load data |
|
| 44 | 46 |
cerfiles=list.files(summarydatadir,pattern="CER_mean_.*tif$",full=T); cerfiles |
| 45 | 47 |
cer=brick(stack(cerfiles)) |
| 46 | 48 |
setZ(cer,months,name="time") |
| ... | ... | |
| 71 | 73 |
cldm=mean(cld,na.rm=T) |
| 72 | 74 |
### TODO: change to bilinear if reprojecting! |
| 73 | 75 |
|
| 76 |
cer20files=list.files(summarydatadir,pattern="CER_P20um_.*tif$",full=T); cer20files |
|
| 77 |
cer20=brick(stack(cer20files)) |
|
| 78 |
setZ(cer20,months,name="time") |
|
| 79 |
cer20@z=list(months) |
|
| 80 |
cer20@zname="time" |
|
| 81 |
layerNames(cer20) <- as.character(format(months,"%b")) |
|
| 82 |
#cot=projectRaster(from=cot,crs="+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0",method="ngb") |
|
| 83 |
cotm=mean(cot,na.rm=T) |
|
| 84 |
### TODO: change to bilinear! |
|
| 85 |
|
|
| 86 |
|
|
| 74 | 87 |
### load PRISM data for comparison |
| 75 | 88 |
prism=brick("data/prism/prism_climate.nc",varname="ppt")
|
| 76 | 89 |
## project to sinusoidal |
| ... | ... | |
| 83 | 96 |
layerNames(prism) <- as.character(format(months,"%b")) |
| 84 | 97 |
|
| 85 | 98 |
#### build a pixel by variable matrix |
| 86 |
vars=c("cer","cld","cot","prism")
|
|
| 99 |
vars=c("cer","cer20","cld","cot","prism")
|
|
| 87 | 100 |
bd=melt(as.matrix(vars[1])) |
| 88 | 101 |
colnames(bd)=c("cell","month",vars[1])
|
| 89 | 102 |
for(v in vars[-1]) {print(v); bd[,v]=melt(as.matrix(get(v)))$value}
|
| ... | ... | |
| 116 | 129 |
|
| 117 | 130 |
### extract MOD06 data for each station |
| 118 | 131 |
stcer=extract(cer,st2)#;colnames(stcer)=paste("cer_mean_",1:12,sep="")
|
| 132 |
stcer20=extract(cer20,st2)#;colnames(stcer)=paste("cer_mean_",1:12,sep="")
|
|
| 119 | 133 |
stcot=extract(cot,st2)#;colnames(stcot)=paste("cot_mean_",1:12,sep="")
|
| 120 | 134 |
stcld=extract(cld,st2)#;colnames(stcld)=paste("cld_mean_",1:12,sep="")
|
| 121 |
mod06=cbind.data.frame(id=st2$id,lat=st2$lat,lon=st2$lon,stcer,stcot,stcld) |
|
| 135 |
mod06=cbind.data.frame(id=st2$id,lat=st2$lat,lon=st2$lon,stcer,stcer20,stcot,stcld)
|
|
| 122 | 136 |
mod06l=melt(mod06,id.vars=c("id","lon","lat"))
|
| 123 | 137 |
mod06l[,c("variable","moment","month")]=do.call(rbind,strsplit(as.character(mod06l$variable),"_"))
|
| 124 | 138 |
mod06l=as.data.frame(cast(mod06l,id+lon+lat+month~variable+moment,value="value")) |
| ... | ... | |
| 158 | 172 |
|
| 159 | 173 |
## melt it |
| 160 | 174 |
mod06sl=melt(mod06s[,!grepl("lppt",colnames(mod06s))],id.vars=c("id","lon","lat","elev","month","ppt","glon","glon2","gelev","gbin"))
|
| 161 |
levels(mod06sl$variable)=c("Effective Radius (um)","Cloudy Days (%)","Optical Thickness (%)","Liquid Water Path")
|
|
| 175 |
levels(mod06sl$variable)=c("Effective Radius (um)","Very Cloudy Days (%)","Cloudy Days (%)","Optical Thickness (%)","Liquid Water Path")
|
|
| 162 | 176 |
|
| 163 | 177 |
################################################################### |
| 164 | 178 |
################################################################### |
| ... | ... | |
| 234 | 248 |
print(p) |
| 235 | 249 |
|
| 236 | 250 |
### monthly comparisons of variables |
| 237 |
mod06sl=melt(mod06s,measure.vars=c("value","COT_mean","CER_mean"))
|
|
| 238 |
bwplot(value~month|variable,data=mod06sl,cex=.5,pch=16,col="black",scales=list(y=list(relation="free")),layout=c(1,3)) |
|
| 239 |
splom(mod06s[grep("CER|COT|CLD",colnames(mod06s))],cex=.2,pch=16,main="Scatterplot matrix of MOD06 products")
|
|
| 251 |
#mod06sl=melt(mod06s,measure.vars=c("ppt","COT_mean","CER_mean","CER_P20um"))
|
|
| 252 |
#bwplot(value~month|variable,data=mod06sl,cex=.5,pch=16,col="black",scales=list(y=list(relation="free")),layout=c(1,3))
|
|
| 253 |
#splom(mod06s[grep("CER|COT|CLD|ppt",colnames(mod06s))],cex=.2,pch=16,main="Scatterplot matrix of MOD06 products")
|
|
| 240 | 254 |
|
| 241 | 255 |
### run some regressions |
| 242 | 256 |
#plot(log(ppt)~COT_mean,data=mod06s) |
| ... | ... | |
| 246 | 260 |
xyplot(ppt~value|variable,groups=glon,data=mod06sl, |
| 247 | 261 |
scales=list(y=list(log=T),x=list(relation="free",log=F)), |
| 248 | 262 |
par.settings = list(superpose.symbol = list(col=bgyr(5),pch=16,cex=.5)),auto.key=list(space="top",title="Station Longitude"), |
| 249 |
main="Comparison of MOD06_L2 and Precipitation Monthly Climatologies",ylab="Mean Monthly Station Precipitation (mm)",xlab="MOD06_L2 Product",layout=c(4,1))+
|
|
| 263 |
main="Comparison of MOD06_L2 and Precipitation Monthly Climatologies",ylab="Mean Monthly Station Precipitation (mm)",xlab="MOD06_L2 Product",layout=c(5,1))+
|
|
| 250 | 264 |
layer(panel.text(9,2.5,label="Coastal stations",srt=30,cex=1.3,col="blue"),columns=1)+ |
| 251 |
layer(panel.text(13,.9,label="Inland stations",srt=10,cex=1.3,col="red"),columns=1) |
|
| 265 |
layer(panel.text(13,.9,label="Inland stations",srt=10,cex=1.3,col="red"),columns=1)+ |
|
| 266 |
layer(panel.abline(lm(y~x),col="red"))+ |
|
| 267 |
layer(panel.text(0,0,paste("R2=",round(summary(lm(y~x))$r.squared,2)),pos=4,cex=.5,col="grey"))
|
|
| 268 |
|
|
| 252 | 269 |
|
| 253 | 270 |
## ppt~metric with longitude bins |
| 254 | 271 |
#CairoPNG("output/COT.png",width=10,height=5,units="in",dpi=300,pointsize=20)
|
| ... | ... | |
| 257 | 274 |
at=quantile(as.matrix(subset(m01,subset=3)),seq(0,1,len=100),na.rm=T) |
| 258 | 275 |
p1=levelplot(subset(m01,subset=3),xlab.top="Optical Thickness (%)",at=at,col.regions=bgyr(length(at)),margin=F, |
| 259 | 276 |
)+layer(sp.lines(roi_geo, lwd=1.2, col='black'))+layer(sp.points(st2, cex=.5,col='black')) |
| 277 |
at=quantile(as.matrix(subset(m01,subset=3)),seq(0,1,len=100),na.rm=T) |
|
| 278 |
|
|
| 260 | 279 |
at=quantile(as.matrix(subset(m01,subset=4)),seq(0,1,len=100),na.rm=T) |
| 261 | 280 |
p2=levelplot(subset(m01,subset=4),xlab.top="PRISM MAP",at=at,col.regions=bgyr(length(at)),margin=F, |
| 262 | 281 |
)+layer(sp.lines(roi_geo, lwd=1.2, col='black'))+layer(sp.points(st2, cex=.5, col='black')) |
| ... | ... | |
| 268 | 287 |
layer(panel.text(9,2.6,label="Coastal stations",srt=10,cex=1.3,col="blue"),columns=1)+ |
| 269 | 288 |
layer(panel.text(13,.95,label="Inland stations",srt=10,cex=1.3,col="red"),columns=1) |
| 270 | 289 |
|
| 271 |
save(p1,p2,p3,file="plotdata.Rdata") |
|
| 272 |
load("plotdata.Rdata")
|
|
| 290 |
p4=xyplot(ppt~value,groups=glon,data=mod06sl[mod06sl$variable=="Very Cloudy Days (%)",], |
|
| 291 |
scales=list(y=list(log=T),x=list(relation="free",log=F)), |
|
| 292 |
par.settings = list(superpose.symbol = list(col=bgyr(5),pch=16,cex=.3)),auto.key=list(space="right",title="Station \n Longitude"), |
|
| 293 |
ylab="Mean Monthly Station Precipitation (mm)",xlab="Proportion days with Cloud Effective Radius >20um from MOD06_L2 (%)",layout=c(1,1))+ |
|
| 294 |
layer(panel.text(9,2.6,label="Coastal stations",srt=10,cex=1.3,col="blue"),columns=1)+ |
|
| 295 |
layer(panel.text(13,.95,label="Inland stations",srt=10,cex=1.3,col="red"),columns=1) |
|
| 296 |
|
|
| 297 |
#save(p1,p2,p3,file="plotdata.Rdata") |
|
| 298 |
#load("plotdata.Rdata")
|
|
| 273 | 299 |
|
| 274 |
CairoPDF("output/MOD06_Summaryfig.pdf",width=11,height=8.5)
|
|
| 300 |
#CairoPDF("output/MOD06_Summaryfig.pdf",width=11,height=8.5)
|
|
| 275 | 301 |
print(p3,position=c(0,0,1,.5),save.object=F) |
| 302 |
#print(p4,position=c(0,0,1,.5),save.object=F) |
|
| 276 | 303 |
print(p1,split=c(1,1,2,2),new=F) |
| 277 | 304 |
print(p2,split=c(2,1,2,2),new=F) |
| 278 |
dev.off() |
|
| 279 |
system("convert output/MOD06_Summaryfig.pdf output/MOD06_Summaryfig.png")
|
|
| 305 |
#dev.off()
|
|
| 306 |
#system("convert output/MOD06_Summaryfig.pdf output/MOD06_Summaryfig.png")
|
|
| 280 | 307 |
#dev.off() |
| 281 | 308 |
|
| 282 | 309 |
## with elevation |
Also available in: Unified diff
Added initial version of LST_Landcover exploration (not finished). Also moved prior version of interpolation procedure (bayesian krig) into repository.