| 1 |
e1df6e22
|
Adam M. Wilson
|
###################################################################################
|
| 2 |
|
|
### R code to aquire and process MOD06_L2 cloud data from the MODIS platform
|
| 3 |
|
|
|
| 4 |
|
|
|
| 5 |
|
|
## connect to server of choice
|
| 6 |
|
|
#system("ssh litoria")
|
| 7 |
|
|
#R
|
| 8 |
|
|
|
| 9 |
|
|
library(sp)
|
| 10 |
|
|
library(spgrass6)
|
| 11 |
|
|
library(rgdal)
|
| 12 |
|
|
library(reshape)
|
| 13 |
|
|
library(ncdf4)
|
| 14 |
|
|
library(geosphere)
|
| 15 |
|
|
library(rgeos)
|
| 16 |
|
|
library(multicore)
|
| 17 |
|
|
library(raster)
|
| 18 |
|
|
library(lattice)
|
| 19 |
|
|
library(rgl)
|
| 20 |
|
|
library(hdf5)
|
| 21 |
|
|
library(rasterVis)
|
| 22 |
|
|
library(heR.Misc)
|
| 23 |
|
|
library(car)
|
| 24 |
|
|
|
| 25 |
|
|
X11.options(type="Xlib")
|
| 26 |
|
|
ncores=20 #number of threads to use
|
| 27 |
|
|
|
| 28 |
|
|
setwd("/home/adamw/acrobates/projects/interp")
|
| 29 |
|
|
|
| 30 |
|
|
psin=CRS("+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs")
|
| 31 |
|
|
|
| 32 |
|
|
## get MODLAND tile information
|
| 33 |
|
|
tb=read.table("http://landweb.nascom.nasa.gov/developers/sn_tiles/sn_bound_10deg.txt",skip=6,nrows=648,header=T)
|
| 34 |
|
|
tb$tile=paste("h",sprintf("%02d",tb$ih),"v",sprintf("%02d",tb$iv),sep="")
|
| 35 |
|
|
save(tb,file="modlandTiles.Rdata")
|
| 36 |
|
|
|
| 37 |
|
|
tile="h11v08" #can move this to submit script if needed
|
| 38 |
|
|
#tile="h09v04" #oregon
|
| 39 |
|
|
|
| 40 |
|
|
tile_bb=tb[tb$tile==tile,] ## identify tile of interest
|
| 41 |
|
|
roi_ll=extent(tile_bb$lon_min,tile_bb$lon_max,tile_bb$lat_min,tile_bb$lat_max)
|
| 42 |
|
|
#roi=spTransform(roi,psin)
|
| 43 |
|
|
#roil=as(roi,"SpatialLines")
|
| 44 |
|
|
|
| 45 |
|
|
dmod06="data/modis/mod06/summary"
|
| 46 |
|
|
|
| 47 |
|
|
|
| 48 |
|
|
##########################
|
| 49 |
|
|
#### explore the data
|
| 50 |
|
|
|
| 51 |
|
|
months=seq(as.Date("2000-01-15"),as.Date("2000-12-15"),by="month")
|
| 52 |
|
|
|
| 53 |
|
|
getmod06<-function(variable){
|
| 54 |
|
|
d=brick(list.files(dmod06,pattern=paste("MOD06_",tile,".nc",sep=""),full=T),varname=toupper(variable))
|
| 55 |
|
|
projection(d)=psin
|
| 56 |
|
|
setZ(d,format(as.Date(d@z$Date),"%m"),name="time")
|
| 57 |
|
|
# d@z=list(months)
|
| 58 |
|
|
layerNames(d) <- as.character(format(as.Date(d@z$Date),"%b"))
|
| 59 |
|
|
return(d)
|
| 60 |
|
|
}
|
| 61 |
|
|
|
| 62 |
|
|
cer=getmod06("cer")
|
| 63 |
|
|
cld=getmod06("cld")
|
| 64 |
|
|
cot=getmod06("cot")
|
| 65 |
|
|
|
| 66 |
|
|
pcol=colorRampPalette(c("brown","red","yellow","darkgreen"))
|
| 67 |
|
|
#levelplot(cer,col.regions=pcol(20))
|
| 68 |
|
|
|
| 69 |
|
|
## load WorldClim data for comparison (download then uncompress)
|
| 70 |
|
|
#system("wget -P data/worldclim/ http://biogeo.ucdavis.edu/data/climate/worldclim/1_4/grid/cur/prec_30s_bil.zip",wait=F)
|
| 71 |
|
|
#system("wget -P data/worldclim/ http://biogeo.ucdavis.edu/data/climate/worldclim/1_4/grid/cur/alt_30s_bil.zip",wait=F)
|
| 72 |
|
|
|
| 73 |
|
|
### load WORLDCLIM data for comparison
|
| 74 |
|
|
wc=stack(list.files("data/worldclim/prec_30s_bil/",pattern="bil$",full=T)[c(4:12,1:3)])
|
| 75 |
|
|
projection(wc)=CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")
|
| 76 |
|
|
wc=crop(wc,roi_ll)
|
| 77 |
|
|
wc[wc==55537]=NA
|
| 78 |
|
|
wc=projectRaster(wc,cer)#crs=projection(psin))
|
| 79 |
|
|
setZ(wc,months,name="time")
|
| 80 |
|
|
wc@z=list(months)
|
| 81 |
|
|
layerNames(wc) <- as.character(format(months,"%b"))
|
| 82 |
|
|
writeRaster(wc,file=paste("data/tiles/",tile,"/worldclim_",tile,".tif",sep=""),format="GTiff")
|
| 83 |
|
|
|
| 84 |
|
|
### load WORLDCLIM elevation
|
| 85 |
|
|
dem=raster(list.files("data/worldclim/alt_30s_bil/",pattern="bil$",full=T))
|
| 86 |
|
|
projection(dem)=CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")
|
| 87 |
|
|
dem=crop(dem,roi_ll)
|
| 88 |
|
|
dem[dem>60000]=NA
|
| 89 |
|
|
dem=projectRaster(dem,cer)
|
| 90 |
|
|
writeRaster(dem,file=paste("data/tiles/",tile,"/dem_",tile,".tif",sep=""),format="GTiff")
|
| 91 |
|
|
|
| 92 |
|
|
|
| 93 |
|
|
### get station data, subset to stations in region, and transform to sinusoidal
|
| 94 |
|
|
dm=readOGR(paste("data/tiles/",tile,sep=""),paste("station_monthly_",tile,"_PRCP",sep=""))
|
| 95 |
|
|
xyplot(latitude~longitude|month,data=dm@data)
|
| 96 |
|
|
dm2=spTransform(dm,CRS(projection(cer)))
|
| 97 |
|
|
dm2@data[,c("x","y")]=coordinates(dm2)
|
| 98 |
|
|
|
| 99 |
|
|
### extract MOD06 data for each station
|
| 100 |
|
|
stcer=extract(cer,dm2,fun=mean);colnames(stcer)=paste("cer_mean_",as.numeric(format(as.Date(cer@z$Date),"%m")),sep="")
|
| 101 |
|
|
#stcer20=extract(cer20,st2)#;colnames(stcer)=paste("cer_mean_",1:12,sep="")
|
| 102 |
|
|
stcot=extract(cot,dm2);colnames(stcot)=paste("cot_mean_",as.numeric(format(as.Date(cot@z$Date),"%m")),sep="")
|
| 103 |
|
|
stcld=extract(cld,dm2);colnames(stcld)=paste("cld_mean_",as.numeric(format(as.Date(cld@z$Date),"%m")),sep="")
|
| 104 |
|
|
stdem=extract(dem,dm2)
|
| 105 |
|
|
mod06=cbind.data.frame(station=dm$station,stcer,stcot,stcld)
|
| 106 |
|
|
mod06l=melt(mod06,id.vars=c("station"));colnames(mod06l)[grep("value",colnames(mod06l))]="mod06"
|
| 107 |
|
|
mod06l[,c("variable","moment","month")]=do.call(rbind,strsplit(as.character(mod06l$variable),"_"))
|
| 108 |
|
|
mod06l=unique(mod06l)
|
| 109 |
|
|
mod06l=cast(mod06l,station+moment+month~variable,value="mod06")
|
| 110 |
|
|
mod06l=merge(dm2@data,mod06l,by=c("station","month"))
|
| 111 |
|
|
mod06l=mod06l[!is.na(mod06l$cer),]
|
| 112 |
|
|
|
| 113 |
|
|
#mod06l=melt(mod06,id.vars=c("station","longitude","latitude","elevation","month","count","value"))
|
| 114 |
|
|
#mod06l[,c("variable","moment","month2")]=do.call(rbind,strsplit(as.character(mod06l$variable),"_"))
|
| 115 |
|
|
#mod06l=as.data.frame(cast(mod06l,station+longitude+latitude+month~variable,value="value.1"))
|
| 116 |
|
|
#mod06l=mod06l[mod06l$month==mod06l$month2&!is.na(mod06l$value.1),]
|
| 117 |
|
|
mod06l=mod06l[order(mod06l$month),]
|
| 118 |
|
|
|
| 119 |
|
|
xyplot(value~cer|month,data=mod06l,scales=list(relation="free"),pch=16,cex=.5)
|
| 120 |
|
|
xyplot(value~cer|station,data=mod06l[mod06l$count>400,],pch=16,cex=.5)
|
| 121 |
|
|
|
| 122 |
|
|
xyplot(cot~month,groups=station,data=mod06l,type="l")
|
| 123 |
|
|
|
| 124 |
|
|
|
| 125 |
|
|
## explore fit of simple model
|
| 126 |
|
|
m=11
|
| 127 |
|
|
cor(mod06l[mod06l$month==m,c("value","cer","cld","cot")])
|
| 128 |
|
|
lm1=lm(value~latitude+longitude+elevation+cer+cld+cot,data=mod06l[mod06l$month==m,])
|
| 129 |
|
|
summary(lm1)
|
| 130 |
|
|
crPlots(lm1)
|
| 131 |
|
|
|
| 132 |
|
|
plot(mod06l$value[mod06l$month==m],as.vector(predict(lm1,data=mod06l[mod06l$month==m,])))
|
| 133 |
|
|
|
| 134 |
|
|
### draw some plots
|
| 135 |
|
|
gq=function(x,n=10,cut=F) {
|
| 136 |
|
|
if(!cut) return(unique(quantile(x,seq(0,1,len=n+1),na.rm=T)))
|
| 137 |
|
|
if(cut) return(cut(x,unique(quantile(x,seq(0,1,len=n+1),na.rm=T))))
|
| 138 |
|
|
}
|
| 139 |
|
|
|
| 140 |
|
|
### add some additional variables
|
| 141 |
|
|
mod06s$month=factor(mod06s$month,labels=format(as.Date(paste("2000",1:12,"15",sep="-")),"%b"))
|
| 142 |
|
|
mod06s$lppt=log(mod06s$ppt)
|
| 143 |
|
|
mod06s$glon=cut(mod06s$lon,gq(mod06s$lon,n=5),include.lowest=T,ordered=T)#gq(mod06s$lon,n=3))
|
| 144 |
|
|
mod06s$glon2=cut(mod06s$lon,breaks=c(-125,-122,-115),labels=c("Coastal","Inland"),include.lowest=T,ordered=T)#gq(mod06s$lon,n=3))
|
| 145 |
|
|
mod06s$gelev=cut(mod06s$elev,breaks=gq(mod06s$elev,n=3),labels=c("Low","Mid","High"),include.lowest=T,ordered=T)
|
| 146 |
|
|
mod06s$gbin=factor(paste(mod06s$gelev,mod06s$glon2,sep="_"),levels=c("Low_Coastal","Mid_Coastal","High_Coastal","Low_Inland","Mid_Inland","High_Inland"),ordered=T)
|
| 147 |
|
|
mod06s$LWP_mean=(2/3)*mod06s$CER_mean*mod06s$COT_mean
|
| 148 |
|
|
|
| 149 |
|
|
## melt it
|
| 150 |
|
|
mod06sl=melt(mod06s[,!grepl("lppt",colnames(mod06s))],id.vars=c("id","lon","lat","elev","month","ppt","glon","glon2","gelev","gbin"))
|
| 151 |
|
|
levels(mod06sl$variable)=c("Effective Radius (um)","Very Cloudy Days (%)","Cloudy Days (%)","Optical Thickness (%)","Liquid Water Path")
|
| 152 |
|
|
|
| 153 |
|
|
###################################################################
|
| 154 |
|
|
###################################################################
|
| 155 |
|
|
|
| 156 |
|
|
bgyr=colorRampPalette(c("blue","green","yellow","red"))
|
| 157 |
|
|
|
| 158 |
|
|
X11.options(type="cairo")
|
| 159 |
|
|
pdf("output/MOD06_summary.pdf",width=11,height=8.5)
|
| 160 |
|
|
|
| 161 |
|
|
# % cloudy maps
|
| 162 |
|
|
title="Cloudiness (% cloudy days) "
|
| 163 |
|
|
at=unique(quantile(as.matrix(cld),seq(0,1,len=100),na.rm=T))
|
| 164 |
|
|
p=levelplot(cld,xlab.top=title,at=at,col.regions=bgyr(length(at)))+layer(sp.lines(roil, lwd=1.2, col='black'))
|
| 165 |
|
|
print(p)
|
| 166 |
|
|
#bwplot(cer,main=title,ylab="Cloud Effective Radius (microns)")
|
| 167 |
|
|
|
| 168 |
|
|
# CER maps
|
| 169 |
|
|
title="Cloud Effective Radius (microns)"
|
| 170 |
|
|
at=quantile(as.matrix(cer),seq(0,1,len=100),na.rm=T)
|
| 171 |
|
|
p=levelplot(cer,xlab.top=title,at=at,col.regions=bgyr(length(at)))+layer(sp.lines(roil, lwd=1.2, col='black'))
|
| 172 |
|
|
print(p)
|
| 173 |
|
|
#bwplot(cer,main=title,ylab="Cloud Effective Radius (microns)")
|
| 174 |
|
|
|
| 175 |
|
|
# COT maps
|
| 176 |
|
|
title="Cloud Optical Thickness (%)"
|
| 177 |
|
|
at=quantile(as.matrix(cot),seq(0,1,len=100),na.rm=T)
|
| 178 |
|
|
p=levelplot(cot,xlab.top=title,at=at,col.regions=bgyr(length(at)))+layer(sp.lines(roil, lwd=0.8, col='black'))
|
| 179 |
|
|
print(p)
|
| 180 |
|
|
#bwplot(cot,xlab.top=title,ylab="Cloud Optical Thickness (%)")
|
| 181 |
|
|
|
| 182 |
|
|
###########################################
|
| 183 |
|
|
#### compare with PRISM data
|
| 184 |
|
|
|
| 185 |
|
|
at=quantile(as.matrix(subset(m01,subset=1)),seq(0,1,len=100),na.rm=T)
|
| 186 |
|
|
p1=levelplot(subset(m01,subset=1),xlab.top="Effective Radius (um)",at=at,col.regions=bgyr(length(at)),margin=F,
|
| 187 |
|
|
)+layer(sp.lines(roi_geo, lwd=1.2, col='black'))
|
| 188 |
|
|
at=quantile(as.matrix(subset(m01,subset=2)),seq(0,1,len=100),na.rm=T)
|
| 189 |
|
|
p2=levelplot(subset(m01,subset=2),xlab.top="Cloudy days (%)",at=at,col.regions=bgyr(length(at)),margin=F,
|
| 190 |
|
|
)+layer(sp.lines(roi_geo, lwd=1.2, col='black'))
|
| 191 |
|
|
at=quantile(as.matrix(subset(m01,subset=3)),seq(0,1,len=100),na.rm=T)
|
| 192 |
|
|
p3=levelplot(subset(m01,subset=3),xlab.top="Optical Thickness (%)",at=at,col.regions=bgyr(length(at)),margin=F,
|
| 193 |
|
|
)+layer(sp.lines(roi_geo, lwd=1.2, col='black'))
|
| 194 |
|
|
at=quantile(as.matrix(subset(m01,subset=4)),seq(0,1,len=100),na.rm=T)
|
| 195 |
|
|
p4=levelplot(subset(m01,subset=4),xlab.top="PRISM MAP",at=at,col.regions=bgyr(length(at)),margin=F,
|
| 196 |
|
|
)+layer(sp.lines(roi_geo, lwd=1.2, col='black'))
|
| 197 |
|
|
|
| 198 |
|
|
print(p1,split=c(1,1,2,2))
|
| 199 |
|
|
print(p2,split=c(1,2,2,2),new=F)
|
| 200 |
|
|
print(p3,split=c(2,1,2,2),new=F)
|
| 201 |
|
|
print(p4,split=c(2,2,2,2),new=F)
|
| 202 |
|
|
|
| 203 |
|
|
### compare COT and PRISM
|
| 204 |
|
|
print(p3,split=c(1,1,2,1))
|
| 205 |
|
|
print(p4,split=c(2,1,2,1),new=F)
|
| 206 |
|
|
|
| 207 |
|
|
|
| 208 |
|
|
### sample to speed up processing
|
| 209 |
|
|
s=sample(1:nrow(bd),10000)
|
| 210 |
|
|
|
| 211 |
|
|
## melt it to ease comparisons
|
| 212 |
|
|
bdl=melt(bd[s,],measure.vars=c("cer","cld","cot"))
|
| 213 |
|
|
|
| 214 |
|
|
combineLimits(useOuterStrips(xyplot(prism~value|variable+month,data=bdl,pch=16,cex=.2,scales=list(y=list(log=T),x=list(relation="free")),
|
| 215 |
|
|
ylab="PRISM Monthly mean precipitation (mm)",xlab="MOD06 metric",main="PRISM vs. MOD06 (mean monthly ppt)")))+
|
| 216 |
|
|
layer(panel.abline(lm(y~x),col="red"))+
|
| 217 |
|
|
layer(panel.text(0,2.5,paste("R2=",round(summary(lm(y~x))$r.squared,2))))
|
| 218 |
|
|
|
| 219 |
|
|
|
| 220 |
|
|
### Comparison at station values
|
| 221 |
|
|
at=quantile(as.matrix(cotm),seq(0,1,len=100),na.rm=T)
|
| 222 |
|
|
p=levelplot(cotm, layers=1,at=at,col.regions=bgyr(length(at)),main="Mean Annual Cloud Optical Thickness",FUN.margin=function(x) 0)+
|
| 223 |
|
|
layer(sp.lines(roil, lwd=1.2, col='black'))+layer(sp.points(st2_sin, pch=16, col='black'))
|
| 224 |
|
|
print(p)
|
| 225 |
|
|
|
| 226 |
|
|
### monthly comparisons of variables
|
| 227 |
|
|
#mod06sl=melt(mod06s,measure.vars=c("ppt","COT_mean","CER_mean","CER_P20um"))
|
| 228 |
|
|
#bwplot(value~month|variable,data=mod06sl,cex=.5,pch=16,col="black",scales=list(y=list(relation="free")),layout=c(1,3))
|
| 229 |
|
|
#splom(mod06s[grep("CER|COT|CLD|ppt",colnames(mod06s))],cex=.2,pch=16,main="Scatterplot matrix of MOD06 products")
|
| 230 |
|
|
|
| 231 |
|
|
### run some regressions
|
| 232 |
|
|
#plot(log(ppt)~COT_mean,data=mod06s)
|
| 233 |
|
|
#summary(lm(log(ppt)~COT_mean*month,data=mod06s))
|
| 234 |
|
|
|
| 235 |
|
|
## ppt~metric with longitude bins
|
| 236 |
|
|
xyplot(ppt~value|variable,groups=glon,data=mod06sl,
|
| 237 |
|
|
scales=list(y=list(log=T),x=list(relation="free",log=F)),
|
| 238 |
|
|
par.settings = list(superpose.symbol = list(col=bgyr(5),pch=16,cex=.5)),auto.key=list(space="top",title="Station Longitude"),
|
| 239 |
|
|
main="Comparison of MOD06_L2 and Precipitation Monthly Climatologies",ylab="Mean Monthly Station Precipitation (mm)",xlab="MOD06_L2 Product",layout=c(5,1))+
|
| 240 |
|
|
layer(panel.text(9,2.5,label="Coastal stations",srt=30,cex=1.3,col="blue"),columns=1)+
|
| 241 |
|
|
layer(panel.text(13,.9,label="Inland stations",srt=10,cex=1.3,col="red"),columns=1)+
|
| 242 |
|
|
layer(panel.abline(lm(y~x),col="red"))+
|
| 243 |
|
|
layer(panel.text(0,0,paste("R2=",round(summary(lm(y~x))$r.squared,2)),pos=4,cex=.5,col="grey"))
|
| 244 |
|
|
|
| 245 |
|
|
|
| 246 |
|
|
## ppt~metric with longitude bins
|
| 247 |
|
|
#CairoPNG("output/COT.png",width=10,height=5,units="in",dpi=300,pointsize=20)
|
| 248 |
|
|
#png("output/COT.png",width=10,height=5,units="in",res=150)
|
| 249 |
|
|
#trellis.par.set("fontsize",12)
|
| 250 |
|
|
at=quantile(as.matrix(subset(m01,subset=3)),seq(0,1,len=100),na.rm=T)
|
| 251 |
|
|
p1=levelplot(subset(m01,subset=3),xlab.top="Optical Thickness (%)",at=at,col.regions=bgyr(length(at)),margin=F,
|
| 252 |
|
|
)+layer(sp.lines(roi_geo, lwd=1.2, col='black'))+layer(sp.points(st2, cex=.5,col='black'))
|
| 253 |
|
|
at=quantile(as.matrix(subset(m01,subset=3)),seq(0,1,len=100),na.rm=T)
|
| 254 |
|
|
|
| 255 |
|
|
at=quantile(as.matrix(subset(m01,subset=4)),seq(0,1,len=100),na.rm=T)
|
| 256 |
|
|
p2=levelplot(subset(m01,subset=4),xlab.top="PRISM MAP",at=at,col.regions=bgyr(length(at)),margin=F,
|
| 257 |
|
|
)+layer(sp.lines(roi_geo, lwd=1.2, col='black'))+layer(sp.points(st2, cex=.5, col='black'))
|
| 258 |
|
|
|
| 259 |
|
|
p3=xyplot(ppt~value,groups=glon,data=mod06sl[mod06sl$variable=="Optical Thickness (%)",],
|
| 260 |
|
|
scales=list(y=list(log=T),x=list(relation="free",log=F)),
|
| 261 |
|
|
par.settings = list(superpose.symbol = list(col=bgyr(5),pch=16,cex=.3)),auto.key=list(space="right",title="Station \n Longitude"),
|
| 262 |
|
|
ylab="Mean Monthly Station Precipitation (mm)",xlab="Cloud Optical Thickness from MOD06_L2 (%)",layout=c(1,1))+
|
| 263 |
|
|
layer(panel.text(9,2.6,label="Coastal stations",srt=10,cex=1.3,col="blue"),columns=1)+
|
| 264 |
|
|
layer(panel.text(13,.95,label="Inland stations",srt=10,cex=1.3,col="red"),columns=1)
|
| 265 |
|
|
|
| 266 |
|
|
p4=xyplot(ppt~value,groups=glon,data=mod06sl[mod06sl$variable=="Very Cloudy Days (%)",],
|
| 267 |
|
|
scales=list(y=list(log=T),x=list(relation="free",log=F)),
|
| 268 |
|
|
par.settings = list(superpose.symbol = list(col=bgyr(5),pch=16,cex=.3)),auto.key=list(space="right",title="Station \n Longitude"),
|
| 269 |
|
|
ylab="Mean Monthly Station Precipitation (mm)",xlab="Proportion days with Cloud Effective Radius >20um from MOD06_L2 (%)",layout=c(1,1))+
|
| 270 |
|
|
layer(panel.text(9,2.6,label="Coastal stations",srt=10,cex=1.3,col="blue"),columns=1)+
|
| 271 |
|
|
layer(panel.text(13,.95,label="Inland stations",srt=10,cex=1.3,col="red"),columns=1)
|
| 272 |
|
|
|
| 273 |
|
|
#save(p1,p2,p3,file="plotdata.Rdata")
|
| 274 |
|
|
#load("plotdata.Rdata")
|
| 275 |
|
|
|
| 276 |
|
|
#CairoPDF("output/MOD06_Summaryfig.pdf",width=11,height=8.5)
|
| 277 |
|
|
print(p3,position=c(0,0,1,.5),save.object=F)
|
| 278 |
|
|
#print(p4,position=c(0,0,1,.5),save.object=F)
|
| 279 |
|
|
print(p1,split=c(1,1,2,2),new=F)
|
| 280 |
|
|
print(p2,split=c(2,1,2,2),new=F)
|
| 281 |
|
|
#dev.off()
|
| 282 |
|
|
#system("convert output/MOD06_Summaryfig.pdf output/MOD06_Summaryfig.png")
|
| 283 |
|
|
#dev.off()
|
| 284 |
|
|
|
| 285 |
|
|
## with elevation
|
| 286 |
|
|
# xyplot(ppt~value|variable,groups=gbin,data=mod06sl,
|
| 287 |
|
|
# scales=list(y=list(log=T),x=list(relation="free")),
|
| 288 |
|
|
# par.settings = list(superpose.symbol = list(col=c(rep("blue",3),rep("red",3)),pch=rep(c(3,4,8),2),cex=.5)),auto.key=list(space="right",title="Station Longitude"),
|
| 289 |
|
|
# main="Comparison of MOD06_L2 and Precipitation Monthly Climatologies",ylab="Precipitation",xlab="MOD06_L2 Product",layout=c(3,1))+
|
| 290 |
|
|
# layer(panel.text(9,2.5,label="Coastal stations",srt=30,cex=1.3,col="blue"),columns=1)+
|
| 291 |
|
|
# layer(panel.text(13,.9,label="Inland stations",srt=10,cex=1.3,col="red"),columns=1)
|
| 292 |
|
|
|
| 293 |
|
|
## with elevation and longitude bins
|
| 294 |
|
|
combineLimits(useOuterStrips(xyplot(ppt~value|variable+gbin,data=mod06sl,
|
| 295 |
|
|
scales=list(y=list(log=T),x=list(relation="free")),col="black",pch=16,cex=.5,type=c("p","r"),
|
| 296 |
|
|
main="Comparison of MOD06_L2 and Precipitation Monthly Climatologies",ylab="Precipitation",xlab="MOD06_L2 Product")))+
|
| 297 |
|
|
layer(panel.xyplot(x,y,type="r",col="red"))
|
| 298 |
|
|
|
| 299 |
|
|
## *** MOD06 vars vs precipitation by month, colored by longitude
|
| 300 |
|
|
combineLimits(useOuterStrips(xyplot(ppt~value|month+variable,groups=glon,data=mod06sl,cex=.5,pch=16,
|
| 301 |
|
|
scales=list(y=list(log=T),x=list(relation="free")),
|
| 302 |
|
|
par.settings = list(superpose.symbol = list(pch =16, col=bgyr(5),cex=1)),auto.key=list(space="top",title="Station Longitude"),
|
| 303 |
|
|
main="Comparison of MOD06_L2 and Precipitation Monthly Climatologies",ylab="Precipitation",xlab="MOD06_L2 Product")),
|
| 304 |
|
|
margin.x=1,adjust.labels=F)+
|
| 305 |
|
|
layer(panel.abline(lm(y~x),col="red"))+
|
| 306 |
|
|
layer(panel.text(0,0,paste("R2=",round(summary(lm(y~x))$r.squared,2)),pos=4,cex=.5,col="grey"))
|
| 307 |
|
|
|
| 308 |
|
|
|
| 309 |
|
|
xyplot(ppt~CLD_mean|id,data=mod06s,panel=function(x,y,group){
|
| 310 |
|
|
panel.xyplot(x,y,type=c("r"),cex=.5,pch=16,col="red")
|
| 311 |
|
|
panel.xyplot(x,y,type=c("p"),cex=.5,pch=16,col="black")
|
| 312 |
|
|
} ,scales=list(y=list(log=T)),strip=F,main="Monthly Mean Precipitation and % Cloudy by station",
|
| 313 |
|
|
sub="Each panel is a station, each point is a monthly mean",
|
| 314 |
|
|
ylab="Precipitation (mm, log axis)",xlab="% of Cloudy Days")+
|
| 315 |
|
|
layer(panel.text(.5,.5,round(summary(lm(y~x))$r.squared,2),pos=4,cex=.75,col="grey"))
|
| 316 |
|
|
|
| 317 |
|
|
xyplot(ppt~CER_mean|id,data=mod06s,panel=function(x,y,group){
|
| 318 |
|
|
panel.xyplot(x,y,type=c("r"),cex=.5,pch=16,col="red")
|
| 319 |
|
|
panel.xyplot(x,y,type=c("p"),cex=.5,pch=16,col="black")
|
| 320 |
|
|
} ,scales=list(y=list(log=T)),strip=F,main="Monthly Mean Precipitation and Cloud Effective Radius by station",sub="Each panel is a station, each point is a monthly mean",ylab="Precipitation (mm, log axis)",xlab="Mean Monthly Cloud Effective Radius (mm)")+
|
| 321 |
|
|
layer(panel.text(10,.5,round(summary(lm(y~x))$r.squared,2),pos=4,cex=.75,col="grey"))
|
| 322 |
|
|
|
| 323 |
|
|
xyplot(ppt~COT_mean|id,data=mod06s,panel=function(x,y,group){
|
| 324 |
|
|
panel.xyplot(x,y,type=c("r"),cex=.5,pch=16,col="red")
|
| 325 |
|
|
panel.xyplot(x,y,type=c("p"),cex=.5,pch=16,col="black")
|
| 326 |
|
|
} ,scales=list(y=list(log=T)),strip=F,main="Monthly Mean Precipitation and Cloud Optical Thickness by station",
|
| 327 |
|
|
sub="Each panel is a station, each point is a monthly mean \n Number in lower right of each panel is R^2",
|
| 328 |
|
|
ylab="Precipitation (mm, log axis)",xlab="Mean Monthly Cloud Optical Thickness (%)")+
|
| 329 |
|
|
layer(panel.text(10,.5,round(summary(lm(y~x))$r.squared,2),pos=4,cex=.75,col="grey"))
|
| 330 |
|
|
|
| 331 |
|
|
xyplot(ppt~LWP_mean|id,data=mod06s,panel=function(x,y,group){
|
| 332 |
|
|
panel.xyplot(x,y,type=c("r"),cex=.5,pch=16,col="red")
|
| 333 |
|
|
panel.xyplot(x,y,type=c("p"),cex=.5,pch=16,col="black")
|
| 334 |
|
|
} ,scales=list(y=list(log=T)),strip=F,main="Monthly Mean Precipitation and Liquid Water Path by station",
|
| 335 |
|
|
sub="Each panel is a station, each point is a monthly mean \n Number in lower right of each panel is R^2",
|
| 336 |
|
|
ylab="Precipitation (mm, log axis)",xlab="Mean Monthly Liquid Water Path")+
|
| 337 |
|
|
layer(panel.text(10,.5,round(summary(lm(y~x))$r.squared,2),pos=4,cex=.75,col="grey"))
|
| 338 |
|
|
|
| 339 |
|
|
### Calculate the slope of each line
|
| 340 |
|
|
mod06s.sl=dapply(mod06s,list(id=mod06s$id),function(x){
|
| 341 |
|
|
lm1=lm(log(x$ppt)~x$CER_mean,)
|
| 342 |
|
|
data.frame(lat=x$lat[1],lon=x$lon[1],elev=x$elev[1],intcpt=coefficients(lm1)[1],cer=coefficients(lm1)[2],r2=summary(lm1)$r.squared)
|
| 343 |
|
|
})
|
| 344 |
|
|
mod06s.sl$cex=gq(mod06s.sl$r2,n=5,cut=T)
|
| 345 |
|
|
mod06s.sl$cer.s=gq(mod06s.sl$cer,n=5,cut=T)
|
| 346 |
|
|
|
| 347 |
|
|
### and plot it on a map
|
| 348 |
|
|
xyplot(lat~lon,group=cer.s,data=mod06s.sl,par.settings = list(superpose.symbol = list(pch =16, col=bgyr(5),cex=1)),auto.key=list(space="right",title="Slope Coefficient"),asp=1,
|
| 349 |
|
|
main="Slopes of linear regressions {log(ppt)~CloudEffectiveRadius}")+
|
| 350 |
|
|
layer(sp.lines(roi_geo, lwd=1.2, col='black'))
|
| 351 |
|
|
|
| 352 |
|
|
### look for relationships with longitude
|
| 353 |
|
|
xyplot(cer~lon,group=cut(mod06s.sl$elev,gq(mod06s.sl$elev,n=5)),data=mod06s.sl,
|
| 354 |
|
|
par.settings = list(superpose.symbol = list(col=bgyr(5),pch=16,cex=1)),auto.key=list(space="right",title="Station Elevation"),
|
| 355 |
|
|
ylab="Slope of lm(ppt~EffectiveRadius)",xlab="Longitude",main="Precipitation~Effective Radius relationship by latitude")
|
| 356 |
|
|
|
| 357 |
|
|
|
| 358 |
|
|
############################################################
|
| 359 |
|
|
### simple regression to get spatial residuals
|
| 360 |
|
|
m="01"
|
| 361 |
|
|
mod06s2=mod06s#[mod06s$month==m,]
|
| 362 |
|
|
|
| 363 |
|
|
lm1=lm(log(ppt)~CER_mean*month*lon,data=mod06s2); summary(lm1)
|
| 364 |
|
|
mod06s2$pred=exp(predict(lm1,mod06s2))
|
| 365 |
|
|
mod06s2$resid=mod06s2$pred-mod06s2$ppt
|
| 366 |
|
|
mod06s2$residg=gq(mod06s2$resid,n=5,cut=T)
|
| 367 |
|
|
mod06s2$presid=mod06s2$resid/mod06s2$ppt
|
| 368 |
|
|
|
| 369 |
|
|
for(l in c(F,T)){
|
| 370 |
|
|
## all months
|
| 371 |
|
|
xyplot(pred~ppt,groups=gelev,data=mod06s2,
|
| 372 |
|
|
par.settings = list(superpose.symbol = list(col=bgyr(3),pch=16,cex=.75)),auto.key=list(space="right",title="Station Elevation"),
|
| 373 |
|
|
scales=list(log=l),
|
| 374 |
|
|
ylab="Predicted Mean Monthly Precipitation (mm)",xlab="Observed Mean Monthly Precipitation (mm)",main="Predicted vs. Observed for Simple Model",
|
| 375 |
|
|
sub="Red line is y=x")+
|
| 376 |
|
|
layer(panel.abline(0,1,col="red"))
|
| 377 |
|
|
|
| 378 |
|
|
## month by month
|
| 379 |
|
|
print(xyplot(pred~ppt|month,groups=gelev,data=mod06s2,
|
| 380 |
|
|
par.settings = list(superpose.symbol = list(col=bgyr(3),pch=16,cex=.75)),auto.key=list(space="right",title="Station Elevation"),
|
| 381 |
|
|
scales=list(log=l),
|
| 382 |
|
|
ylab="Predicted Mean Monthly Precipitation (mm)",xlab="Observed Mean Monthly Precipitation (mm)",main="Predicted vs. Observed for Simple Model",
|
| 383 |
|
|
sub="Red line is y=x")+
|
| 384 |
|
|
layer(panel.abline(0,1,col="red"))
|
| 385 |
|
|
)}
|
| 386 |
|
|
|
| 387 |
|
|
## residuals by month
|
| 388 |
|
|
xyplot(lat~lon|month,group=residg,data=mod06s2,
|
| 389 |
|
|
par.settings = list(superpose.symbol = list(pch =16, col=bgyr(5),cex=.5)),
|
| 390 |
|
|
auto.key=list(space="right",title="Residuals"),
|
| 391 |
|
|
main="Spatial plot of monthly residuals")+
|
| 392 |
|
|
layer(sp.lines(roi_geo, lwd=1.2, col='black'))
|
| 393 |
|
|
|
| 394 |
|
|
|
| 395 |
|
|
dev.off()
|
| 396 |
|
|
|
| 397 |
|
|
|
| 398 |
|
|
####################################
|
| 399 |
|
|
#### build table comparing various metrics
|
| 400 |
|
|
mods=data.frame(
|
| 401 |
|
|
models=c(
|
| 402 |
|
|
"log(ppt)~CER_mean",
|
| 403 |
|
|
"log(ppt)~CLD_mean",
|
| 404 |
|
|
"log(ppt)~COT_mean",
|
| 405 |
|
|
"log(ppt)~CER_mean*month",
|
| 406 |
|
|
"log(ppt)~CLD_mean*month",
|
| 407 |
|
|
"log(ppt)~COT_mean*month",
|
| 408 |
|
|
"log(ppt)~CER_mean*month*lon",
|
| 409 |
|
|
"log(ppt)~CLD_mean*month*lon",
|
| 410 |
|
|
"log(ppt)~COT_mean*month*lon",
|
| 411 |
|
|
"ppt~CER_mean*month*lon",
|
| 412 |
|
|
"ppt~CLD_mean*month*lon",
|
| 413 |
|
|
"ppt~COT_mean*month*lon"),stringsAsFactors=F)
|
| 414 |
|
|
|
| 415 |
|
|
mods$r2=
|
| 416 |
|
|
do.call(rbind,lapply(1:nrow(mods),function(i){
|
| 417 |
|
|
lm1=lm(as.formula(mods$models[i]),data=mod06s2)
|
| 418 |
|
|
summary(lm1)$r.squared}))
|
| 419 |
|
|
|
| 420 |
|
|
mods
|
| 421 |
|
|
|
| 422 |
|
|
|
| 423 |
|
|
|
| 424 |
|
|
|
| 425 |
|
|
|
| 426 |
|
|
|
| 427 |
|
|
|
| 428 |
|
|
|
| 429 |
|
|
load("data/modis/pointsummary.Rdata")
|
| 430 |
|
|
|
| 431 |
|
|
|
| 432 |
|
|
dsl=melt(ds,id.vars=c("id","date","ppt","lon","lat"),measure.vars= c("Cloud_Water_Path","Cloud_Effective_Radius","Cloud_Optical_Thickness"))
|
| 433 |
|
|
|
| 434 |
|
|
dsl=dsl[!is.nan(dsl$value),]
|
| 435 |
|
|
|
| 436 |
|
|
|
| 437 |
|
|
|
| 438 |
|
|
|
| 439 |
|
|
####
|
| 440 |
|
|
## mean annual precip
|
| 441 |
|
|
dp=d[d$variable=="ppt",]
|
| 442 |
|
|
dp$year=format(dp$date,"%Y")
|
| 443 |
|
|
dm=tapply(dp$value,list(id=dp$id,year=dp$year),sum,na.rm=T)
|
| 444 |
|
|
dms=apply(dm,1,mean,na.rm=T)
|
| 445 |
|
|
dms=data.frame(id=names(dms),ppt=dms/10)
|
| 446 |
|
|
|
| 447 |
|
|
dslm=tapply(dsl$value,list(id=dsl$id,variable=dsl$variable),mean,na.rm=T)
|
| 448 |
|
|
dslm=data.frame(id=rownames(dslm),dslm)
|
| 449 |
|
|
|
| 450 |
|
|
dms=merge(dms,dslm)
|
| 451 |
|
|
dmsl=melt(dms,id.vars=c("id","ppt"))
|
| 452 |
|
|
|
| 453 |
|
|
summary(lm(ppt~Cloud_Effective_Radius,data=dms))
|
| 454 |
|
|
summary(lm(ppt~Cloud_Water_Path,data=dms))
|
| 455 |
|
|
summary(lm(ppt~Cloud_Optical_Thickness,data=dms))
|
| 456 |
|
|
summary(lm(ppt~Cloud_Effective_Radius+Cloud_Water_Path+Cloud_Optical_Thickness,data=dms))
|
| 457 |
|
|
|
| 458 |
|
|
|
| 459 |
|
|
#### draw some plots
|
| 460 |
|
|
#pdf("output/MOD06_summary.pdf",width=11,height=8.5)
|
| 461 |
|
|
png("output/MOD06_summary_%d.png",width=1024,height=780)
|
| 462 |
|
|
|
| 463 |
|
|
## daily data
|
| 464 |
|
|
xyplot(value~ppt/10|variable,data=dsl,
|
| 465 |
|
|
scales=list(relation="free"),type=c("p","r"),
|
| 466 |
|
|
pch=16,cex=.5,layout=c(3,1))
|
| 467 |
|
|
|
| 468 |
|
|
|
| 469 |
|
|
densityplot(~value|variable,groups=cut(dsl$ppt,c(0,50,100,500)),data=dsl,auto.key=T,
|
| 470 |
|
|
scales=list(relation="free"),plot.points=F)
|
| 471 |
|
|
|
| 472 |
|
|
## annual means
|
| 473 |
|
|
|
| 474 |
|
|
xyplot(value~ppt|variable,data=dmsl,
|
| 475 |
|
|
scales=list(relation="free"),type=c("p","r"),pch=16,cex=0.5,layout=c(3,1),
|
| 476 |
|
|
xlab="Mean Annual Precipitation (mm)",ylab="Mean value")
|
| 477 |
|
|
|
| 478 |
|
|
densityplot(~value|variable,groups=cut(dsl$ppt,c(0,50,100,500)),data=dmsl,auto.key=T,
|
| 479 |
|
|
scales=list(relation="free"),plot.points=F)
|
| 480 |
|
|
|
| 481 |
|
|
|
| 482 |
|
|
## plot some swaths
|
| 483 |
|
|
|
| 484 |
|
|
nc1=raster(fs$path[3],varname="Cloud_Effective_Radius")
|
| 485 |
|
|
nc2=raster(fs$path[4],varname="Cloud_Effective_Radius")
|
| 486 |
|
|
nc3=raster(fs$path[5],varname="Cloud_Effective_Radius")
|
| 487 |
|
|
|
| 488 |
|
|
nc1[nc1<=0]=NA
|
| 489 |
|
|
nc2[nc2<=0]=NA
|
| 490 |
|
|
nc3[nc3<=0]=NA
|
| 491 |
|
|
|
| 492 |
|
|
plot(roi)
|
| 493 |
|
|
plot(nc3)
|
| 494 |
|
|
|
| 495 |
|
|
plot(nc1,add=T)
|
| 496 |
|
|
plot(nc2,add=T)
|
| 497 |
|
|
|
| 498 |
|
|
|
| 499 |
|
|
dev.off()
|
| 500 |
|
|
|
