|
1
|
## explore the MOD35 data downloaded and gridded by the DAAC
|
|
2
|
setwd("~/acrobates/adamw/projects/interp/data/modis/mod35")
|
|
3
|
|
|
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
|
projection(v5)="+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs"
|
|
32
|
v6=stack(brick("summary/MOD35_h11v08.nc",varname="PCloud"))
|
|
33
|
projection(v6)="+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs"
|
|
34
|
v6=setZ(v6,as.Date(paste("2011-",1:12,"-15",sep="")))
|
|
35
|
names(v6)=month.name
|
|
36
|
|
|
37
|
## generate means
|
|
38
|
v6m=mean(v6)
|
|
39
|
v5m=mean(v5)
|
|
40
|
|
|
41
|
|
|
42
|
## landcover
|
|
43
|
lulc=raster("~/acrobatesroot/jetzlab/Data/environ/global/landcover/MODIS/MCD12Q1_IGBP_2005_v51.tif")
|
|
44
|
projection(lulc)="+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs"
|
|
45
|
lulc=crop(lulc,v6)
|
|
46
|
|
|
47
|
Mode <- function(x,na.rm=T) { #get MODE
|
|
48
|
x=na.omit(x)
|
|
49
|
ux <- unique(x)
|
|
50
|
ux[which.max(tabulate(match(x, ux)))]
|
|
51
|
}
|
|
52
|
## aggregate to 1km resolution
|
|
53
|
lulc2=aggregate(lulc,2,fun=function(x,na.rm=T) Mode(x))
|
|
54
|
## convert to factor table
|
|
55
|
lulcf=lulc2
|
|
56
|
ratify(lulcf)
|
|
57
|
levels(lulcf)[[1]]
|
|
58
|
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")
|
|
59
|
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")
|
|
60
|
|
|
61
|
levels(lulcf)=list(data.frame(ID=0:16,LULC=lulc_levels,LULC2=lulc_levels2))
|
|
62
|
|
|
63
|
|
|
64
|
### load WORLDCLIM elevation
|
|
65
|
dem=raster("../../tiles/h11v08/dem_h11v08.tif",format="GTiff")
|
|
66
|
projection(dem)="+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs"
|
|
67
|
|
|
68
|
dif=v6-v5
|
|
69
|
names(dif)=month.name
|
|
70
|
|
|
71
|
difm=v6m-v5m
|
|
72
|
v5v6compare=stack(v5m,v6m,difm)
|
|
73
|
names(v5v6compare)=c("Collection 5","Collection 6","Difference (C6-C5)")
|
|
74
|
|
|
75
|
tile=extent(v6)
|
|
76
|
|
|
77
|
### compare differences between v5 and v6 by landcover type
|
|
78
|
lulcm=as.matrix(lulc)
|
|
79
|
forest=lulcm>=1&lulcm<=5
|
|
80
|
|
|
81
|
|
|
82
|
#####################################
|
|
83
|
### compare MOD43 and MOD17 products
|
|
84
|
|
|
85
|
## MOD17
|
|
86
|
mod17=raster("../MOD17/Npp_1km_C5.1_mean_00_to_06.tif",format="GTiff")
|
|
87
|
mod17=crop(projectRaster(mod17,v6,method="bilinear"),v6)
|
|
88
|
NAvalue(mod17)=32767
|
|
89
|
|
|
90
|
mod17qc=raster("../MOD17/Npp_QC_1km_C5.1_mean_00_to_06.tif",format="GTiff")
|
|
91
|
mod17qc=crop(projectRaster(mod17qc,v6,method="bilinear"),v6)
|
|
92
|
mod17qc[mod17qc<0|mod17qc>100]=NA
|
|
93
|
|
|
94
|
## MOD43 via earth engine
|
|
95
|
mod43=raster("../mod43/3b21aa90cc657523ff31e9559f36fb12.EVI_MEAN.tif",format="GTiff")
|
|
96
|
mod43=crop(projectRaster(mod43,v6,method="bilinear"),v6)
|
|
97
|
|
|
98
|
mod43qc=raster("../mod43/3b21aa90cc657523ff31e9559f36fb12.Percent_Cloudy.tif",format="GTiff")
|
|
99
|
mod43qc=crop(projectRaster(mod43qc,v6,method="bilinear"),v6)
|
|
100
|
mod43qc[mod43qc<0|mod43qc>100]=NA
|
|
101
|
|
|
102
|
## Summary plot of mod17 and mod43
|
|
103
|
modprod=stack(mod17qc,mod43qc)
|
|
104
|
names(modprod)=c("MOD17","MOD43")
|
|
105
|
|
|
106
|
###
|
|
107
|
|
|
108
|
n=100
|
|
109
|
at=seq(0,100,len=n)
|
|
110
|
cols=grey(seq(0,1,len=n))
|
|
111
|
cols=rainbow(n)
|
|
112
|
bgyr=colorRampPalette(c("blue","green","yellow","red"))
|
|
113
|
cols=bgyr(n)
|
|
114
|
|
|
115
|
#levelplot(lulcf,margin=F,layers="LULC")
|
|
116
|
|
|
117
|
m=3
|
|
118
|
mcompare=stack(subset(v5,m),subset(v6,m))
|
|
119
|
|
|
120
|
mdiff=subset(v5,m)-subset(v6,m)
|
|
121
|
names(mcompare)=c("Collection_5","Collection_6")
|
|
122
|
names(mdiff)=c("Collection_5-Collection_6")
|
|
123
|
|
|
124
|
|
|
125
|
CairoPDF("output/mod35compare.pdf",width=11,height=8.5)
|
|
126
|
#CairoPNG("output/mod35compare_%d.png",units="in", width=11,height=8.5,pointsize=4000,dpi=1200,antialias="subpixel")
|
|
127
|
|
|
128
|
levelplot(mcompare,col.regions=cols,at=at,margin=F,sub="Frequency of MOD35 Clouds in March")
|
|
129
|
#levelplot(dif,col.regions=bgyr(20),margin=F)
|
|
130
|
levelplot(mdiff,col.regions=bgyr(100),at=seq(mdiff@data@min,mdiff@data@max,len=100),margin=F)
|
|
131
|
|
|
132
|
|
|
133
|
boxplot(as.matrix(subset(dif,subset=1))~forest,varwidth=T,notch=T);abline(h=0)
|
|
134
|
|
|
135
|
|
|
136
|
levelplot(modprod,main="Missing Data (%) in MOD17 (NPP) and MOD43 (BRDF Reflectance)",
|
|
137
|
sub="Tile H11v08 (Venezuela)",col.regions=cols,at=at)
|
|
138
|
|
|
139
|
levelplot(modprod,main="Missing Data (%) in MOD17 (NPP) and MOD43 (BRDF Reflectance)",
|
|
140
|
sub="Tile H11v08 (Venezuela)",col.regions=cols,at=at,
|
|
141
|
xlim=c(-7300000,-6670000),ylim=c(0,600000))
|
|
142
|
|
|
143
|
levelplot(v5m,main="Missing Data (%) in MOD17 (NPP) and MOD43 (BRDF Reflectance)",
|
|
144
|
sub="Tile H11v08 (Venezuela)",col.regions=cols,at=at,
|
|
145
|
xlim=c(-7200000,-6670000),ylim=c(0,400000),margin=F)
|
|
146
|
|
|
147
|
|
|
148
|
levelplot(subset(v5v6compare,1:2),main="Proportion Cloudy Days (%) in Collection 5 and 6 MOD35",
|
|
149
|
sub="Tile H11v08 (Venezuela)",col.regions=cols,at=at,
|
|
150
|
margin=F)
|
|
151
|
|
|
152
|
levelplot(subset(v5v6compare,1:2),main="Proportion Cloudy Days (%) in Collection 5 and 6 MOD35",
|
|
153
|
sub="Tile H11v08 (Venezuela)",col.regions=cols,at=at,
|
|
154
|
xlim=c(-7200000,-6670000),ylim=c(0,400000),margin=F)
|
|
155
|
|
|
156
|
levelplot(subset(v5v6compare,1:2),main="Proportion Cloudy Days (%) in Collection 5 and 6 MOD35",
|
|
157
|
sub="Tile H11v08 (Venezuela)",col.regions=cols,at=at,
|
|
158
|
xlim=c(-7500000,-7200000),ylim=c(700000,1000000),margin=F)
|
|
159
|
|
|
160
|
|
|
161
|
dev.off()
|
|
162
|
|
|
163
|
### smoothing plots
|
|
164
|
## explore smoothed version
|
|
165
|
td=subset(v6,m)
|
|
166
|
## build weight matrix
|
|
167
|
s=3
|
|
168
|
w=matrix(1/(s*s),nrow=s,ncol=s)
|
|
169
|
#w[s-1,s-1]=4/12; w
|
|
170
|
td2=focal(td,w=w)
|
|
171
|
td3=stack(td,td2)
|
|
172
|
|
|
173
|
levelplot(td3,col.regions=cols,at=at,margin=F)
|
|
174
|
|
|
175
|
dev.off()
|
|
176
|
plot(stack(difm,lulc))
|
|
177
|
|
|
178
|
### ROI
|
|
179
|
tile_ll=projectExtent(v6, "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
|
|
180
|
|
|
181
|
62,59
|
|
182
|
0,3
|
|
183
|
|
|
184
|
|
|
185
|
|
|
186
|
#### export KML timeseries
|
|
187
|
library(plotKML)
|
|
188
|
tile="h11v08"
|
|
189
|
file=paste("summary/MOD35_",tile,".nc",sep="")
|
|
190
|
system(paste("gdalwarp -overwrite -multi -ot INT16 -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=""))
|
|
191
|
|
|
192
|
v6sp=brick(paste("MOD35_",tile,".tif",sep=""))
|
|
193
|
v6sp=readAll(v6sp)
|
|
194
|
|
|
195
|
## wasn't working with line below, perhaps Z should just be text? not date?
|
|
196
|
v6sp=setZ(v6sp,as.Date(paste("2011-",1:12,"-15",sep="")))
|
|
197
|
names(v6sp)=month.name
|
|
198
|
|
|
199
|
kml_open("output/mod35.kml")
|
|
200
|
|
|
201
|
|
|
202
|
kml_layer.RasterBrick(v6sp,
|
|
203
|
plot.legend = TRUE, dtime = "", tz = "GMT",
|
|
204
|
z.lim = c(0,100),colour_scale = get("colour_scale_numeric", envir = plotKML.opts))
|
|
205
|
# home_url = get("home_url", envir = plotKML.opts),
|
|
206
|
# metadata = NULL, html.table = NULL,
|
|
207
|
# altitudeMode = "clampToGround", balloon = FALSE,
|
|
208
|
)
|
|
209
|
|
|
210
|
logo = "http://static.tumblr.com/t0afs9f/KWTm94tpm/yale_logo.png"
|
|
211
|
kml_screen(image.file = logo, position = "UL", sname = "YALE logo",size=c(.1,.1))
|
|
212
|
kml_close("mod35.kml")
|
|
213
|
kml_compress("mod35.kml",files=c(paste(month.name,".png",sep=""),"obj_legend.png"),zip="/usr/bin/zip")
|
