|
1
|
## explore the MOD35 data downloaded and gridded by the DAAC
|
|
2
|
setwd("~/acrobates/adamw/projects/interp/data/modis/mod35")
|
|
3
|
setwd("/Users/adamw/GoogleDrive/Work/presentations/20130426_NASAMeeting/data")
|
|
4
|
library(raster)
|
|
5
|
library(rasterVis)
|
|
6
|
library(rgdal)
|
|
7
|
|
|
8
|
#f=list.files(pattern="*.hdf")
|
|
9
|
|
|
10
|
#Sys.setenv(GEOL_AS_GCPS = "PARTIAL")
|
|
11
|
## try swath-grid with gdal
|
|
12
|
#GDALinfo(f[1])
|
|
13
|
#system(paste("gdalinfo",f[2]))
|
|
14
|
#GDALinfo("HDF4_EOS:EOS_SWATH:\"MOD35_L2.A2000100.1445.006.2012252024758.hdf\":mod35:Cloud_Mask")
|
|
15
|
#system("gdalinfo HDF4_EOS:EOS_SWATH:\"data/modis/mod35/MOD35_L2.A2000100.1445.006.2012252024758.hdf\":mod35:Cloud_Mask")
|
|
16
|
#system("gdalwarp -overwrite -geoloc -order 2 -r near -s_srs \"EPSG:4326\" HDF4_EOS:EOS_SWATH:\"MOD35_L2.A2000100.1445.006.2012252024758.hdf\":mod35:Cloud_Mask cloudmask.tif")
|
|
17
|
#system("gdalwarp -overwrite -r near -s_srs \"EPSG:4326\" HDF4_EOS:EOS_SWATH:\"MOD35_L2.A2000100.1445.006.2012252024758.hdf\":mod35:Cloud_Mask:1 cloudmask2.tif")
|
|
18
|
#r=raster(f[1])
|
|
19
|
#extent(r)
|
|
20
|
#st=lapply(f[1:10],raster)
|
|
21
|
#str=lapply(2:length(st),function(i) union(extent(st[[i-1]]),extent(st[[i]])))[[length(st)-1]]
|
|
22
|
#str=union(extent(h11v08),str)
|
|
23
|
#b1=brick(lapply(st,function(stt) {
|
|
24
|
# x=crop(alignExtent(stt,str),h11v08)
|
|
25
|
# return(x)
|
|
26
|
#}))
|
|
27
|
#c=brick(f[1:10])
|
|
28
|
|
|
29
|
## get % cloudy
|
|
30
|
v5=stack(brick("../mod06/summary/MOD06_h11v08_ymoncld01.nc",varname="CLD01"))
|
|
31
|
v5=stack(brick("summary/MOD06_h11v08_ymoncld01.nc",varname="CLD01"))
|
|
32
|
projection(v5)="+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs"
|
|
33
|
v6=stack(brick("summary/MOD35_h11v08.nc",varname="PCloud"))
|
|
34
|
projection(v6)="+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs"
|
|
35
|
v6=setZ(v6,as.Date(paste("2011-",1:12,"-15",sep="")))
|
|
36
|
names(v6)=month.name
|
|
37
|
|
|
38
|
## generate means
|
|
39
|
v6m=mean(v6)
|
|
40
|
v5m=mean(v5)
|
|
41
|
|
|
42
|
|
|
43
|
## worldclim
|
|
44
|
wc=brick("summary/worldclim_h11v08.tif")
|
|
45
|
names(wc)=paste("WC",1:12,sep="_")
|
|
46
|
projection(wc)="+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs"
|
|
47
|
wcsp=projectRaster(wc,crs="+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
|
|
48
|
wcsp=readAll(wcsp)
|
|
49
|
wcsp=setZ(wcsp,paste("2011-",1:12,"-15",sep=""))
|
|
50
|
|
|
51
|
## landcover
|
|
52
|
lulc=raster("~/acrobatesroot/jetzlab/Data/environ/global/landcover/MODIS/MCD12Q1_IGBP_2005_v51.tif")
|
|
53
|
lulc=raster("~/acrobatesroot/jetzlab/Data/environ/global/landcover/MODIS/MCD12Q1_IGBP_2005_v51.tif")
|
|
54
|
|
|
55
|
projection(lulc)="+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs"
|
|
56
|
lulc=crop(lulc,v6)
|
|
57
|
|
|
58
|
Mode <- function(x,na.rm=T) { #get MODE
|
|
59
|
x=na.omit(x)
|
|
60
|
ux <- unique(x)
|
|
61
|
ux[which.max(tabulate(match(x, ux)))]
|
|
62
|
}
|
|
63
|
## aggregate to 1km resolution
|
|
64
|
lulc2=aggregate(lulc,2,fun=function(x,na.rm=T) Mode(x))
|
|
65
|
## convert to factor table
|
|
66
|
lulcf=lulc2
|
|
67
|
ratify(lulcf)
|
|
68
|
levels(lulcf)[[1]]
|
|
69
|
lulc_levels=c("Water","Evergreen Needleleaf forest","Evergreen Broadleaf forest","Deciduous Needleleaf forest","Deciduous Broadleaf forest","Mixed forest","Closed shrublands","Open shrublands","Woody savannas","Savannas","Grasslands","Permanent wetlands","Croplands","Urban and built-up","Cropland/Natural vegetation mosaic","Snow and ice","Barren or sparsely vegetated")
|
|
70
|
lulc_levels2=c("Water","Forest","Forest","Forest","Forest","Forest","Shrublands","Shrublands","Savannas","Savannas","Grasslands","Permanent wetlands","Croplands","Urban and built-up","Cropland/Natural vegetation mosaic","Snow and ice","Barren or sparsely vegetated")
|
|
71
|
|
|
72
|
levels(lulcf)=list(data.frame(ID=0:16,LULC=lulc_levels,LULC2=lulc_levels2))
|
|
73
|
|
|
74
|
|
|
75
|
### load WORLDCLIM elevation
|
|
76
|
dem=raster("../../tiles/h11v08/dem_h11v08.tif",format="GTiff")
|
|
77
|
|
|
78
|
projection(dem)="+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs"
|
|
79
|
|
|
80
|
dem=raster("summary/dem_h11v08.tif",format="GTiff")
|
|
81
|
dif=v6-v5
|
|
82
|
names(dif)=month.name
|
|
83
|
|
|
84
|
difm=v6m-v5m
|
|
85
|
v5v6compare=stack(v5m,v6m,difm)
|
|
86
|
names(v5v6compare)=c("Collection 5","Collection 6","Difference (C6-C5)")
|
|
87
|
|
|
88
|
tile=extent(v6)
|
|
89
|
|
|
90
|
### compare differences between v5 and v6 by landcover type
|
|
91
|
lulcm=as.matrix(lulc)
|
|
92
|
forest=lulcm>=1&lulcm<=5
|
|
93
|
|
|
94
|
|
|
95
|
#####################################
|
|
96
|
### compare MOD43 and MOD17 products
|
|
97
|
|
|
98
|
## MOD17
|
|
99
|
#mod17=raster("../MOD17/Npp_1km_C5.1_mean_00_to_06.tif",format="GTiff")
|
|
100
|
#mod17=crop(projectRaster(mod17,v6,method="bilinear"),v6)
|
|
101
|
#NAvalue(mod17)=32767
|
|
102
|
|
|
103
|
mod17qc=raster("../MOD17/Npp_QC_1km_C5.1_mean_00_to_06.tif",format="GTiff")
|
|
104
|
mod17qc=raster("summary/Npp_QC_1km_C5.1_mean_00_to_06.tif",format="GTiff")
|
|
105
|
|
|
106
|
mod17qc=crop(projectRaster(mod17qc,v6,method="bilinear"),v6)
|
|
107
|
mod17qc[mod17qc<0|mod17qc>100]=NA
|
|
108
|
|
|
109
|
## MOD43 via earth engine
|
|
110
|
#mod43=raster("../mod43/3b21aa90cc657523ff31e9559f36fb12.EVI_MEAN.tif",format="GTiff")
|
|
111
|
#mod43=crop(projectRaster(mod43,v6,method="bilinear"),v6)
|
|
112
|
|
|
113
|
mod43qc=raster("../mod43/3b21aa90cc657523ff31e9559f36fb12.Percent_Cloudy.tif",format="GTiff")
|
|
114
|
mod43qc=raster("summary/3b21aa90cc657523ff31e9559f36fb12.Percent_Cloudy.tif",format="GTiff")
|
|
115
|
mod43qc=crop(projectRaster(mod43qc,v6,method="bilinear"),v6)
|
|
116
|
mod43qc[mod43qc<0|mod43qc>100]=NA
|
|
117
|
|
|
118
|
## Summary plot of mod17 and mod43
|
|
119
|
modprod=stack(mod17qc,mod43qc)
|
|
120
|
names(modprod)=c("MOD17","MOD43")
|
|
121
|
|
|
122
|
### worldclim and interpolation
|
|
123
|
m7=raster("summary/h11v08_pred_2_8.nc",varname="ppt")
|
|
124
|
m7=readAll(m7)
|
|
125
|
m7=m7*30
|
|
126
|
pred=stack(subset(wc,2),subset(m7,2))
|
|
127
|
names(pred)=c("WorldClim","MOD35 Precipitation")
|
|
128
|
|
|
129
|
n=100
|
|
130
|
at=seq(0,100,len=n)
|
|
131
|
cols=grey(seq(0,1,len=n))
|
|
132
|
cols=rainbow(n)
|
|
133
|
bgyr=colorRampPalette(c("blue","green","yellow","red"))
|
|
134
|
gbg=colorRampPalette(c(grey(.2),"grey","blue"))
|
|
135
|
|
|
136
|
cols=bgyr(n)
|
|
137
|
|
|
138
|
#levelplot(lulcf,margin=F,layers="LULC")
|
|
139
|
compare=stack(v5,v6)
|
|
140
|
|
|
141
|
m=3
|
|
142
|
mcompare=stack(subset(v5,m),subset(v6,m))
|
|
143
|
names(compare)=paste(rep(c("C5","C6"),each=12),1:12,sep="_")
|
|
144
|
|
|
145
|
mdiff=subset(v5,m)-subset(v6,m)
|
|
146
|
names(mcompare)=c("Collection_5","Collection_6")
|
|
147
|
names(mdiff)=c("Collection_5-Collection_6")
|
|
148
|
|
|
149
|
library(Cairo)
|
|
150
|
(trellis.par.get())
|
|
151
|
trellis.par.set(par.ylab.text=list(cex=5))
|
|
152
|
|
|
153
|
pdf("output/mod35compare.pdf",width=11,height=8.5, fontsize=18)
|
|
154
|
#CairoPNG("output/mod35compare_%d.png",units="in", width=11,height=8.5,pointsize=4000,dpi=1200,antialias="subpixel")
|
|
155
|
|
|
156
|
levelplot(mcompare,col.regions=cols,at=at,margin=F,sub="Cloud Frequency (%) from MOD35 in March")
|
|
157
|
|
|
158
|
## show all months c5 vs c6
|
|
159
|
levelplot(compare,col.regions=cols,at=at,margin=F,sub="C5 and C6 MOD35 Comparison",
|
|
160
|
layout=c(12,2),as.table=T)
|
|
161
|
|
|
162
|
#levelplot(dif,col.regions=bgyr(20),margin=F)
|
|
163
|
levelplot(mdiff,col.regions=bgyr(100),at=seq(mdiff@data@min,mdiff@data@max,len=100),margin=F)
|
|
164
|
|
|
165
|
|
|
166
|
#boxplot(as.matrix(subset(dif,subset=1))~forest,varwidth=T,notch=T);abline(h=0)
|
|
167
|
|
|
168
|
|
|
169
|
levelplot(modprod,main="Missing Data (%) in MOD17 (NPP) and MOD43 (BRDF Reflectance)",
|
|
170
|
sub="Tile H11v08 (Venezuela)",col.regions=cols,at=at)
|
|
171
|
|
|
172
|
levelplot(modprod,main="Missing Data (%) in MOD17 (NPP) and MOD43 (BRDF Reflectance)",
|
|
173
|
sub="Tile H11v08 (Venezuela)",col.regions=cols,at=at,
|
|
174
|
xlim=c(-7300000,-6670000),ylim=c(0,600000))
|
|
175
|
|
|
176
|
levelplot(subset(v5v6compare,1:2),main="Mean Annual Proportion Cloudy Days (%) in Collection 5 and 6 MOD35",
|
|
177
|
sub="Tile H11v08 (Venezuela)",col.regions=cols,at=at,
|
|
178
|
margin=F)
|
|
179
|
|
|
180
|
levelplot(subset(v5v6compare,1:2),main="Mean Annual Proportion Cloudy Days (%) in Collection 5 and 6 MOD35",
|
|
181
|
sub="Tile H11v08 (Venezuela)",col.regions=cols,at=at,
|
|
182
|
xlim=c(-7200000,-6670000),ylim=c(0,400000),margin=F)
|
|
183
|
|
|
184
|
levelplot(subset(v5v6compare,1:2),main="Proportion Cloudy Days (%) in Collection 5 and 6 MOD35",
|
|
185
|
sub="Tile H11v08 (Venezuela)",col.regions=cols,at=at,
|
|
186
|
xlim=c(-7500000,-7200000),ylim=c(700000,1000000),margin=F)
|
|
187
|
|
|
188
|
levelplot(pred,main="Comparison of WorldClim with Cloud-Informed interpolation",
|
|
189
|
sub="Tile H11v08 (Venezuela)",col.regions=grey(seq(0,1,len=n)),at=quantile(as.numeric(as.matrix(pred)),na.rm=T,seq(0,1,len=n)))
|
|
190
|
|
|
191
|
levelplot(pred,main="Comparison of WorldClim with Cloud-Informed interpolation",
|
|
192
|
sub="Tile H11v08 (Venezuela)",col.regions=cols,at=quantile(as.numeric(as.matrix(pred)),na.rm=T,seq(0,1,len=n)))
|
|
193
|
|
|
194
|
levelplot(pred,main="Comparison of WorldClim with Cloud-Informed interpolation",
|
|
195
|
sub="Tile H11v08 (Venezuela)",col.regions=cols,at=seq(0,360,len=n))
|
|
196
|
|
|
197
|
p.strip <- list(cex=1.5, lines=2, fontface='bold')
|
|
198
|
ckey <- list(labels=list(cex=2), height=1)
|
|
199
|
x.scale <- list(cex=1.5, alternating=1)
|
|
200
|
y.scale <- list(cex=1.5, alternating=1)
|
|
201
|
|
|
202
|
levelplot(pred,main="Comparison of WorldClim with Cloud-Informed interpolation",
|
|
203
|
sub="Tile H11v08 (Venezuela)",col.regions=cols,at=c(seq(-2,150,len=n*.75),seq(150.1,360,len=n*.25)),
|
|
204
|
xlim=c(-7700000,-6700000),ylim=c(350000,1100000) , scales=list(x=x.scale,y=y.scale),par.strip.text=p.strip,colorkey=ckey,
|
|
205
|
layout=c(2,1))
|
|
206
|
|
|
207
|
|
|
208
|
dev.off()
|
|
209
|
|
|
210
|
### smoothing plots
|
|
211
|
## explore smoothed version
|
|
212
|
td=subset(v6,m)
|
|
213
|
## build weight matrix
|
|
214
|
s=3
|
|
215
|
w=matrix(1/(s*s),nrow=s,ncol=s)
|
|
216
|
#w[s-1,s-1]=4/12; w
|
|
217
|
td2=focal(td,w=w)
|
|
218
|
td3=stack(td,td2)
|
|
219
|
|
|
220
|
levelplot(td3,col.regions=cols,at=at,margin=F)
|
|
221
|
|
|
222
|
dev.off()
|
|
223
|
plot(stack(difm,lulc))
|
|
224
|
|
|
225
|
### ROI
|
|
226
|
tile_ll=projectExtent(v6, "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
|
|
227
|
|
|
228
|
62,59
|
|
229
|
0,3
|
|
230
|
|
|
231
|
|
|
232
|
|
|
233
|
#### export KML timeseries
|
|
234
|
wd="/Users/adamw/GoogleDrive/Work/presentations/20130426_NASAMeeting/data/"
|
|
235
|
library(plotKML)
|
|
236
|
tile="h11v08"
|
|
237
|
file=paste("summary/MOD35_",tile,".nc",sep="")
|
|
238
|
system(paste("gdalwarp -overwrite -multi -r cubicspline -srcnodata 255 -dstnodata 255 -s_srs '+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs' -t_srs 'EPSG:4326' NETCDF:",file,":PCloud MOD35_",tile,".tif",sep=""))
|
|
239
|
|
|
240
|
v6sp=brick(paste("summary/MOD35_",tile,".tif",sep=""))
|
|
241
|
v6sp=readAll(v6sp)
|
|
242
|
|
|
243
|
v6sp=setZ(v6sp,paste("2011-",1:12,"-15",sep=""))
|
|
244
|
names(v6sp)=month.name
|
|
245
|
|
|
246
|
|
|
247
|
setwd(paste(wd,"/kml/mod35",sep=""))
|
|
248
|
kmlf="mod35.kml"
|
|
249
|
kml_open(file.name=kmlf)
|
|
250
|
|
|
251
|
|
|
252
|
kml_layer.RasterBrick(v6sp,
|
|
253
|
plot.legend = TRUE, dtime = "", tz = "GMT",
|
|
254
|
z.lim = c(0,100),colour_scale = gbg(20),
|
|
255
|
# home_url = get("home_url", envir = plotKML.opts),
|
|
256
|
# metadata = NULL, html.table = NULL,
|
|
257
|
altitudeMode = "clampToGround", balloon = FALSE,folder.name=kmlfl,legend.file="mod35_legend.png"
|
|
258
|
)
|
|
259
|
|
|
260
|
#logo = "yale_logo.png"
|
|
261
|
#kml_screen(image.file = logo, position = "UL", sname = "YALE logo",size=c(.1,.1))
|
|
262
|
kml_close(file.name=kmlf)
|
|
263
|
#kml_compress(file.name=kmlf,files=c(paste(month.name,".png",sep=""),"obj_legend.png"),zip="/usr/bin/zip")
|
|
264
|
|
|
265
|
### worldclim
|
|
266
|
setwd(paste(wd,"/kml/wc",sep=""))
|
|
267
|
kmlf="worldclim.kml"
|
|
268
|
kml_open(file.name=kmlf)
|
|
269
|
|
|
270
|
kml_layer.RasterBrick(wcsp,
|
|
271
|
plot.legend = TRUE, dtime = "", tz = "GMT",
|
|
272
|
z.lim = c(0,650),colour_scale = gbg(100),
|
|
273
|
# home_url = get("home_url", envir = plotKML.opts),
|
|
274
|
# metadata = NULL, html.table = NULL,
|
|
275
|
altitudeMode = "clampToGround", balloon = FALSE
|
|
276
|
)
|
|
277
|
|
|
278
|
kml_close(file.name=kmlf)
|
|
279
|
#kml_compress(file.name=kmlf,files=c(paste("WC_",1:12,".png",sep=""),"obj_legend.png"),zip="/usr/bin/zip")
|
|
280
|
|
|
281
|
setwd(wd)
|
