1
|
###################################################################################
|
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(rgdal)
|
11
|
library(reshape)
|
12
|
library(ncdf4)
|
13
|
library(geosphere)
|
14
|
library(rgeos)
|
15
|
library(multicore)
|
16
|
library(raster)
|
17
|
library(lattice)
|
18
|
library(rgl)
|
19
|
library(hdf5)
|
20
|
|
21
|
X11.options(type="Xlib")
|
22
|
ncores=20 #number of threads to use
|
23
|
|
24
|
setwd("/home/adamw/personal/projects/interp")
|
25
|
setwd("/home/adamw/acrobates/projects/interp")
|
26
|
|
27
|
roi=readOGR("data/regions/Test_sites/Oregon.shp","Oregon")
|
28
|
roi=spTransform(roi,CRS(" +proj=sinu +lon_0=0 +x_0=0 +y_0=0"))
|
29
|
|
30
|
|
31
|
### Get processed MOD06 nc files
|
32
|
datadir="data/modis/MOD06_nc"
|
33
|
|
34
|
fs=data.frame(
|
35
|
path=list.files(datadir,full=T,recursive=T,pattern="nc$"),
|
36
|
file=basename(list.files(datadir,full=F,recursive=T,pattern="nc$")))
|
37
|
fs$date=as.Date(substr(fs$file,11,17),"%Y%j")
|
38
|
fs$time=substr(fs$file,19,22)
|
39
|
fs$datetime=as.POSIXct(strptime(paste(substr(fs$file,11,17),substr(fs$file,19,22)), '%Y%j %H%M'))
|
40
|
fs$path=as.character(fs$path)
|
41
|
fs$file=as.character(fs$file)
|
42
|
|
43
|
|
44
|
### get station data
|
45
|
load("data/ghcn/roi_ghcn.Rdata")
|
46
|
load("data/allstations.Rdata")
|
47
|
|
48
|
|
49
|
d2=d[d$variable=="ppt"&d$date%in%fs$date,]
|
50
|
d2=d2[,-grep("variable",colnames(d2)),]
|
51
|
st2=st[st$id%in%d$id,]
|
52
|
d2[,c("lon","lat")]=coordinates(st2)[match(d2$id,st2$id),]
|
53
|
|
54
|
## generate list split by day to merge with MOD06 data
|
55
|
d2l=split(d2,d2$date)
|
56
|
|
57
|
ds=do.call(rbind.data.frame,mclapply(d2l,function(td){
|
58
|
print(td$date[1])
|
59
|
tfs=fs[fs$date==td$date[1],]
|
60
|
## make spatial object
|
61
|
tdsp=td
|
62
|
coordinates(tdsp)=c("lon","lat")
|
63
|
projection(tdsp)=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
|
64
|
tdsp=spTransform(tdsp, CRS(" +proj=sinu +lon_0=0 +x_0=0 +y_0=0"))
|
65
|
|
66
|
td2=do.call(rbind.data.frame,lapply(1:nrow(tfs),function(i){
|
67
|
tnc1=brick(
|
68
|
raster(tfs$path[i],varname="Cloud_Water_Path"),
|
69
|
raster(tfs$path[i],varname="Cloud_Effective_Radius"),
|
70
|
raster(tfs$path[i],varname="Cloud_Optical_Thickness"))
|
71
|
projection(tnc1)=CRS(" +proj=sinu +lon_0=0 +x_0=0 +y_0=0")
|
72
|
td3=as.data.frame(extract(tnc1,tdsp))
|
73
|
if(ncol(td3)==1) return(NULL)
|
74
|
colnames(td3)=
|
75
|
c("Cloud_Water_Path","Cloud_Effective_Radius","Cloud_Optical_Thickness")
|
76
|
## drop negative values (need to check why these exist)
|
77
|
td3[td3<0]=NA
|
78
|
td3[,c("date","id")]=td[,c("date","id")]
|
79
|
return(td3)
|
80
|
}))
|
81
|
td2=melt(td2,id.vars=c("date","id"))
|
82
|
td2=cast(td2,date+id~variable,fun=mean,na.rm=T)
|
83
|
td2=merge(td,td2)
|
84
|
td2$id=as.character(td2$id)
|
85
|
td2$date=as.character(td2$date)
|
86
|
return(td2)
|
87
|
}))
|
88
|
|
89
|
colnames(ds)[grep("value",colnames(ds))]="ppt"
|
90
|
ds$ppt=as.numeric(ds$ppt)
|
91
|
ds$Cloud_Water_Path=as.numeric(ds$Cloud_Water_Path)
|
92
|
ds$Cloud_Effective_Radius=as.numeric(ds$Cloud_Effective_Radius)
|
93
|
ds$Cloud_Optical_Thickness=as.numeric(ds$Cloud_Optical_Thickness)
|
94
|
ds$date=as.Date(ds$date)
|
95
|
ds$year=as.numeric(ds$year)
|
96
|
ds$lon=as.numeric(ds$lon)
|
97
|
ds$lat=as.numeric(ds$lat)
|
98
|
ds=ds[!is.na(ds$ppt),]
|
99
|
|
100
|
save(ds,file="data/modis/pointsummary.Rdata")
|
101
|
|
102
|
load("data/modis/pointsummary.Rdata")
|
103
|
|
104
|
|
105
|
dsl=melt(ds,id.vars=c("id","date","ppt","lon","lat"),measure.vars= c("Cloud_Water_Path","Cloud_Effective_Radius","Cloud_Optical_Thickness"))
|
106
|
|
107
|
dsl=dsl[!is.nan(dsl$value),]
|
108
|
|
109
|
|
110
|
summary(lm(ppt~Cloud_Effective_Radius,data=ds))
|
111
|
summary(lm(ppt~Cloud_Water_Path,data=ds))
|
112
|
summary(lm(ppt~Cloud_Optical_Thickness,data=ds))
|
113
|
summary(lm(ppt~Cloud_Effective_Radius+Cloud_Water_Path+Cloud_Optical_Thickness,data=ds))
|
114
|
|
115
|
|
116
|
####
|
117
|
## mean annual precip
|
118
|
dp=d[d$variable=="ppt",]
|
119
|
dp$year=format(dp$date,"%Y")
|
120
|
dm=tapply(dp$value,list(id=dp$id,year=dp$year),sum,na.rm=T)
|
121
|
dms=apply(dm,1,mean,na.rm=T)
|
122
|
dms=data.frame(id=names(dms),ppt=dms/10)
|
123
|
|
124
|
dslm=tapply(dsl$value,list(id=dsl$id,variable=dsl$variable),mean,na.rm=T)
|
125
|
dslm=data.frame(id=rownames(dslm),dslm)
|
126
|
|
127
|
dms=merge(dms,dslm)
|
128
|
dmsl=melt(dms,id.vars=c("id","ppt"))
|
129
|
|
130
|
summary(lm(ppt~Cloud_Effective_Radius,data=dms))
|
131
|
summary(lm(ppt~Cloud_Water_Path,data=dms))
|
132
|
summary(lm(ppt~Cloud_Optical_Thickness,data=dms))
|
133
|
summary(lm(ppt~Cloud_Effective_Radius+Cloud_Water_Path+Cloud_Optical_Thickness,data=dms))
|
134
|
|
135
|
|
136
|
#### draw some plots
|
137
|
#pdf("output/MOD06_summary.pdf",width=11,height=8.5)
|
138
|
png("output/MOD06_summary_%d.png",width=1024,height=780)
|
139
|
|
140
|
## daily data
|
141
|
xyplot(value~ppt/10|variable,data=dsl,
|
142
|
scales=list(relation="free"),type=c("p","r"),
|
143
|
pch=16,cex=.5,layout=c(3,1))
|
144
|
|
145
|
|
146
|
densityplot(~value|variable,groups=cut(dsl$ppt,c(0,50,100,500)),data=dsl,auto.key=T,
|
147
|
scales=list(relation="free"),plot.points=F)
|
148
|
|
149
|
## annual means
|
150
|
|
151
|
xyplot(value~ppt|variable,data=dmsl,
|
152
|
scales=list(relation="free"),type=c("p","r"),pch=16,cex=0.5,layout=c(3,1),
|
153
|
xlab="Mean Annual Precipitation (mm)",ylab="Mean value")
|
154
|
|
155
|
densityplot(~value|variable,groups=cut(dsl$ppt,c(0,50,100,500)),data=dmsl,auto.key=T,
|
156
|
scales=list(relation="free"),plot.points=F)
|
157
|
|
158
|
|
159
|
## plot some swaths
|
160
|
|
161
|
nc1=raster(fs$path[3],varname="Cloud_Effective_Radius")
|
162
|
nc2=raster(fs$path[4],varname="Cloud_Effective_Radius")
|
163
|
nc3=raster(fs$path[5],varname="Cloud_Effective_Radius")
|
164
|
|
165
|
nc1[nc1<=0]=NA
|
166
|
nc2[nc2<=0]=NA
|
167
|
nc3[nc3<=0]=NA
|
168
|
|
169
|
plot(roi)
|
170
|
plot(nc3)
|
171
|
|
172
|
plot(nc1,add=T)
|
173
|
plot(nc2,add=T)
|
174
|
|
175
|
|
176
|
dev.off()
|
177
|
|
178
|
|