Project

General

Profile

Download (16.5 KB) Statistics
| Branch: | Revision:
1
Systematic landcover bias in Collection 5 MODIS cloud mask and derived products – a global overview
2
__________
3

    
4
```{r,echo=TRUE,eval=TRUE}
5
 opts_chunk$set(eval=F)
6
```
7

    
8
This document describes the analysis of the Collection 5 MOD35 data.
9

    
10
# Google Earth Engine Processing
11
The following code produces the annual (2009) summaries of cloud frequency from MOD09, MOD35, and MOD11 using the Google Earth Engine 'playground' API [http://ee-api.appspot.com/](http://ee-api.appspot.com/). 
12
```{r, engine='coffee', eval=F}
13
var startdate="2009-01-01"
14
var stopdate="2009-12-31"
15

    
16
// MOD11 MODIS LST
17
var mod11 = ee.ImageCollection("MOD11A2").map(function(img){ 
18
  return img.select(['LST_Day_1km'])});
19
// MOD09 internal cloud flag
20
var mod09 = ee.ImageCollection("MOD09GA").filterDate(new Date(startdate),new Date(stopdate)).map(function(img) {
21
  return img.select(['state_1km']).expression("((b(0)/1024)%2)");
22
});
23
// MOD35 cloud flag
24
var mod35 = ee.ImageCollection("MOD09GA").filterDate(new Date(startdate),new Date(stopdate)).map(function(img) {
25
  return img.select(['state_1km']).expression("((b(0))%4)==1|((b(0))%4)==2");
26
});
27

    
28
//define reducers
29
var COUNT = ee.call("Reducer.count");
30
var MEAN = ee.call("Reducer.mean");
31

    
32
//a few maps of constants
33
c100=ee.Image(100);  //to multiply by 100
34
c02=ee.Image(0.02);  //to scale LST data
35
c272=ee.Image(272.15); // to convert K->C
36

    
37
//calculate mean cloudiness (%), rename, and convert to integer
38
mod09a=mod09.reduce(MEAN).select([0], ['MOD09']).multiply(c100).int8();
39
mod35a=mod35.reduce(MEAN).select([0], ['MOD35']).multiply(c100).int8();
40

    
41
/////////////////////////////////////////////////
42
// Generate the cloud frequency surface:
43
getMiss = function(collection) {
44
  //filter by date
45
i2=collection.filterDate(new Date(startdate),new Date(stopdate));
46
// number of layers in collection
47
i2_n=i2.getInfo().features.length;
48
//get means
49
// image of -1s to convert to % missing
50
c1=ee.Image(-1);
51
// 1 Calculate the number of days with measurements
52
// 2 divide by the total number of layers
53
i2_c=ee.Image(i2_n).float()
54
// 3 Add -1 and multiply by -1 to invert to % cloudy
55
// 4 Rename to "Percent_Cloudy"
56
// 5 multiply by 100 and convert to 8-bit integer to decrease file size
57
i2_miss=i2.reduce(COUNT).divide(i2_c).add(c1).multiply(c1).multiply(c100).int8();
58
return (i2_miss);
59
};
60

    
61
// run the function
62
mod11a=getMiss(mod11).select([0], ['MOD11_LST_PMiss']);
63

    
64
// get long-term mean
65
mod11b=mod11.reduce(MEAN).multiply(c02).subtract(c272).int8().select([0], ['MOD11_LST_MEAN']);
66

    
67
// summary object with all layers
68
summary=mod11a.addBands(mod11b).addBands(mod35a).addBands(mod09a)
69

    
70
var region='[[-180, -60], [-180, 90], [180, 90], [180, -60]]'  //global
71

    
72
// get download link
73
print("All")
74
var path = summary.getDownloadURL({
75
  'scale': 1000,
76
  'crs': 'EPSG:4326',
77
  'region': region
78
});
79
print('https://earthengine.sandbox.google.com' + path);
80
```
81

    
82
# Data Processing
83

    
84
```{r}
85
setwd("~/acrobates/adamw/projects/MOD35C5")
86
library(raster)
87
beginCluster(10)
88
library(rasterVis)
89
library(rgdal)
90
library(plotKML)
91
library(Cairo)
92
library(reshape)
93
library(rgeos)
94
library(splancs)
95

    
96
## get % cloudy
97
mod09=raster("data/MOD09_2009.tif")
98
names(mod09)="C5MOD09CF"
99
NAvalue(mod09)=0
100

    
101
mod35c5=raster("data/MOD35_2009.tif")
102
names(mod35c5)="C5MOD35CF"
103
NAvalue(mod35c5)=0
104

    
105
## mod35C6 annual
106
mod35c6=raster("data/MOD35C6_2009.tif")
107
names(mod35c6)="C6MOD35CF"
108
NAvalue(mod35c6)=255
109

    
110
## landcover
111
lulc=raster("data/MCD12Q1_IGBP_2009_051_wgs84_1km.tif")
112

    
113
#  lulc=ratify(lulc)
114
  data(worldgrids_pal)  #load palette
115
  IGBP=data.frame(ID=0:16,col=worldgrids_pal$IGBP[-c(18,19)],
116
    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"),stringsAsFactors=F)
117
  IGBP$class=rownames(IGBP);rownames(IGBP)=1:nrow(IGBP)
118
  levels(lulc)=list(IGBP)
119
names(lulc)="MCD12Q1"
120

    
121
## MOD17
122
mod17=raster("data/MOD17.tif",format="GTiff")
123
NAvalue(mod17)=65535
124
names(mod17)="MOD17_unscaled"
125

    
126
mod17qc=raster("data/MOD17qc.tif",format="GTiff")
127
NAvalue(mod17qc)=255
128
names(mod17qc)="MOD17CF"
129

    
130
## MOD11 via earth engine
131
mod11=raster("data/MOD11_2009.tif",format="GTiff")
132
names(mod11)="MOD11_unscaled"
133
NAvalue(mod11)=0
134

    
135
mod11qc=raster("data/MOD11qc_2009.tif",format="GTiff")
136
names(mod11qc)="MOD11CF"
137
```
138

    
139
Import the Collection 5 MOD35 processing path:
140
```{r}
141
pp=raster("data/MOD35pp.tif")
142
NAvalue(pp)=255
143
names(pp)="MOD35pp"
144
```
145

    
146
Define transects to illustrate the fine-grain relationship between MOD35 cloud frequency and both landcover and processing path.
147
```{r}
148
r1=Lines(list(
149
  Line(matrix(c(
150
                -61.688,4.098,
151
                -59.251,3.430
152
                ),ncol=2,byrow=T))),"Venezuela")
153
r2=Lines(list(
154
  Line(matrix(c(
155
                133.746,-31.834,
156
                134.226,-32.143
157
                ),ncol=2,byrow=T))),"Australia")
158
r3=Lines(list(
159
  Line(matrix(c(
160
                73.943,27.419,
161
                74.369,26.877
162
                ),ncol=2,byrow=T))),"India")
163
r4=Lines(list(
164
  Line(matrix(c(
165
                33.195,12.512,
166
                33.802,12.894
167
                ),ncol=2,byrow=T))),"Sudan")
168

    
169
trans=SpatialLines(list(r1,r2,r3,r4),CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs "))
170
### write out shapefiles of transects
171
writeOGR(SpatialLinesDataFrame(trans,data=data.frame(ID=names(trans)),match.ID=F),"output",layer="transects",driver="ESRI Shapefile",overwrite=T)
172
```
173

    
174
Buffer transects to identify a small region around each transect for comparison and plotting
175
```{r}
176
transb=gBuffer(trans,byid=T,width=0.4)
177
## make polygons of bounding boxes
178
bb0 <- lapply(slot(transb, "polygons"), bbox)
179
bb1 <- lapply(bb0, bboxx)
180
# turn these into matrices using a helper function in splancs
181
bb2 <- lapply(bb1, function(x) rbind(x, x[1,]))
182
# close the matrix rings by appending the first coordinate
183
rn <- row.names(transb)
184
# get the IDs
185
bb3 <- vector(mode="list", length=length(bb2))
186
# make somewhere to keep the output
187
for (i in seq(along=bb3)) bb3[[i]] <- Polygons(list(Polygon(bb2[[i]])),
188
                   ID=rn[i])
189
# loop over the closed matrix rings, adding the IDs
190
bbs <- SpatialPolygons(bb3, proj4string=CRS(proj4string(transb)))
191
```
192

    
193
Extract the CF and mean values from each raster of interest.
194
```{r}
195
trd1=lapply(1:length(transb),function(x) {
196
  td=crop(mod11,transb[x])
197
  tdd=lapply(list(mod35c5,mod35c6,mod09,mod17,mod17qc,mod11,mod11qc,lulc,pp),function(l) resample(crop(l,transb[x]),td,method="ngb"))
198
  ## normalize MOD11 and MOD17
199
  for(j in which(do.call(c,lapply(tdd,function(i) names(i)))%in%c("MOD11_unscaled","MOD17_unscaled"))){
200
    trange=cellStats(tdd[[j]],range)
201
    tscaled=100*(tdd[[j]]-trange[1])/(trange[2]-trange[1])
202
    tscaled@history=list(range=trange)
203
    names(tscaled)=sub("_unscaled","",names(tdd[[j]]))
204
    tdd=c(tdd,tscaled)
205
  }
206
  return(brick(tdd))
207
})
208
## bind all subregions into single dataframe for plotting
209
trd=do.call(rbind.data.frame,lapply(1:length(trd1),function(i){
210
  d=as.data.frame(as.matrix(trd1[[i]]))
211
  d[,c("x","y")]=coordinates(trd1[[i]])
212
  d$trans=names(trans)[i]
213
  d=melt(d,id.vars=c("trans","x","y"))
214
  return(d)
215
}))
216
transd=do.call(rbind.data.frame,lapply(1:length(trans),function(l) {
217
  td=as.data.frame(extract(trd1[[l]],trans[l],along=T,cellnumbers=F)[[1]])
218
  td$loc=extract(trd1[[l]],trans[l],along=T,cellnumbers=T)[[1]][,1]
219
  td[,c("x","y")]=xyFromCell(trd1[[l]],td$loc)
220
  td$dist=spDistsN1(as.matrix(td[,c("x","y")]), as.matrix(td[1,c("x","y")]),longlat=T)
221
  td$transect=names(trans[l])
222
  td2=melt(td,id.vars=c("loc","x","y","dist","transect"))
223
  td2=td2[order(td2$variable,td2$dist),]
224
  # get per variable ranges to normalize
225
  tr=cast(melt.list(tapply(td2$value,td2$variable,function(x) data.frame(min=min(x,na.rm=T),max=max(x,na.rm=T)))),L1~variable)
226
  td2$min=tr$min[match(td2$variable,tr$L1)]
227
  td2$max=tr$max[match(td2$variable,tr$L1)]
228
  print(paste("Finished ",names(trans[l])))
229
  return(td2)}
230
  ))
231

    
232
transd$type=ifelse(grepl("MOD35|MOD09|CF",transd$variable),"CF","Data")
233
```
234

    
235
Compute difference between MOD09 and MOD35 cloud masks
236
```{r}
237
## comparison of % cloudy days
238
dif_c5_09=raster("data/dif_c5_09.tif",format="GTiff")
239
```
240

    
241
Define a color scheme
242
```{r}
243
n=100
244
at=seq(0,100,len=n)
245
bgyr=colorRampPalette(c("purple","blue","green","yellow","orange","red","red"))
246
bgrayr=colorRampPalette(c("purple","blue","grey","red","red"))
247
cols=bgyr(n)
248
```
249

    
250
Import a global coastline map for overlay
251
```{r} 
252
library(maptools)
253
coast=map2SpatialLines(map("world", interior=FALSE, plot=FALSE),proj4string=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))
254
```
255

    
256
Draw the global cloud frequencies
257
```{r}
258
g1=levelplot(stack(mod35c5,mod09),xlab=" ",scales=list(x=list(draw=F),y=list(alternating=1)),col.regions=cols,at=at)+layer(sp.polygons(bbs[1:4],lwd=2))+layer(sp.lines(coast,lwd=.5))
259

    
260
g2=levelplot(dif_c5_09,col.regions=bgrayr(100),at=seq(-70,70,len=100),margin=F,ylab=" ",colorkey=list("right"))+layer(sp.polygons(bbs[1:4],lwd=2))+layer(sp.lines(coast,lwd=.5))
261
g2$strip=strip.custom(var.name="Difference (C5MOD35-C5MOD09)",style=1,strip.names=T,strip.levels=F)  
262
```
263

    
264
Now illustrate the fine-grain regions
265
```{r}
266
p1=useOuterStrips(levelplot(value~x*y|variable+trans,data=trd[!trd$variable%in%c("MOD17_unscaled","MOD11_unscaled","MCD12Q1","MOD35pp"),],asp=1,scales=list(draw=F,rot=0,relation="free"),
267
                                       at=at,col.regions=cols,maxpixels=7e6,
268
                                       ylab="Latitude",xlab="Longitude"),strip = strip.custom(par.strip.text=list(cex=.7)))+layer(sp.lines(trans,lwd=2))
269

    
270
p2=useOuterStrips(
271
  levelplot(value~x*y|variable+trans,data=trd[trd$variable%in%c("MCD12Q1"),],
272
            asp=1,scales=list(draw=F,rot=0,relation="free"),colorkey=F,
273
            at=c(-1,IGBP$ID),col.regions=IGBP$col,maxpixels=7e7,
274
            legend=list(
275
              right=list(fun=draw.key(list(columns=1,#title="MCD12Q1 \n IGBP Land \n Cover",
276
                           rectangles=list(col=IGBP$col,size=1),
277
                           text=list(as.character(IGBP$ID),at=IGBP$ID-.5))))),
278
            ylab="",xlab=" "),strip = strip.custom(par.strip.text=list(cex=.7)),strip.left=F)+layer(sp.lines(trans,lwd=2))
279
p3=useOuterStrips(
280
  levelplot(value~x*y|variable+trans,data=trd[trd$variable%in%c("MOD35pp"),],
281
            asp=1,scales=list(draw=F,rot=0,relation="free"),colorkey=F,
282
            at=c(-1:4),col.regions=c("blue","cyan","tan","darkgreen"),maxpixels=7e7,
283
            legend=list(
284
              right=list(fun=draw.key(list(columns=1,#title="MOD35 \n Processing \n Path",
285
                           rectangles=list(col=c("blue","cyan","tan","darkgreen"),size=1),
286
                           text=list(c("Water","Coast","Desert","Land")))))),
287
            ylab="",xlab=" "),strip = strip.custom(par.strip.text=list(cex=.7)),strip.left=F)+layer(sp.lines(trans,lwd=2))
288
```
289

    
290
Now draw the profile plots for each transect.
291
```{r}
292
## transects
293
p4=xyplot(value~dist|transect,groups=variable,type=c("smooth","p"),
294
       data=transd,panel=function(...,subscripts=subscripts) {
295
         td=transd[subscripts,]
296
         ## mod09
297
         imod09=td$variable=="C5MOD09CF"
298
         panel.xyplot(td$dist[imod09],td$value[imod09],type=c("p","smooth"),span=0.2,subscripts=1:sum(imod09),col="red",pch=16,cex=.25)
299
         ## mod35C5
300
         imod35=td$variable=="C5MOD35CF"
301
         panel.xyplot(td$dist[imod35],td$value[imod35],type=c("p","smooth"),span=0.09,subscripts=1:sum(imod35),col="blue",pch=16,cex=.25)
302
         ## mod35C6
303
         imod35c6=td$variable=="C6MOD35CF"
304
         panel.xyplot(td$dist[imod35c6],td$value[imod35c6],type=c("p","smooth"),span=0.09,subscripts=1:sum(imod35c6),col="black",pch=16,cex=.25)
305
         ## mod17
306
         imod17=td$variable=="MOD17"
307
         panel.xyplot(td$dist[imod17],100*((td$value[imod17]-td$min[imod17][1])/(td$max[imod17][1]-td$min[imod17][1])),
308
                      type=c("smooth"),span=0.09,subscripts=1:sum(imod17),col="darkgreen",lty=5,pch=1,cex=.25)
309
         imod17qc=td$variable=="MOD17CF"
310
         panel.xyplot(td$dist[imod17qc],td$value[imod17qc],type=c("p","smooth"),span=0.09,subscripts=1:sum(imod17qc),col="darkgreen",pch=16,cex=.25)
311
         ## mod11
312
         imod11=td$variable=="MOD11"
313
         panel.xyplot(td$dist[imod11],100*((td$value[imod11]-td$min[imod11][1])/(td$max[imod11][1]-td$min[imod11][1])),
314
                      type=c("smooth"),span=0.09,subscripts=1:sum(imod17),col="orange",lty="dashed",pch=1,cex=.25)
315
         imod11qc=td$variable=="MOD11CF"
316
         qcspan=ifelse(td$transect[1]=="Australia",0.2,0.05)
317
         panel.xyplot(td$dist[imod11qc],td$value[imod11qc],type=c("p","smooth"),npoints=100,span=qcspan,subscripts=1:sum(imod11qc),col="orange",pch=16,cex=.25)
318
         ## land
319
         path=td[td$variable=="MOD35pp",]
320
         panel.segments(path$dist,-10,c(path$dist[-1],max(path$dist,na.rm=T)),-10,col=c("blue","cyan","tan","darkgreen")[path$value+1],subscripts=1:nrow(path),lwd=10,type="l")
321
         land=td[td$variable=="MCD12Q1",]
322
         panel.segments(land$dist,-20,c(land$dist[-1],max(land$dist,na.rm=T)),-20,col=IGBP$col[land$value+1],subscripts=1:nrow(land),lwd=10,type="l")
323
        },subscripts=T,par.settings = list(grid.pars = list(lineend = "butt")),
324
       scales=list(
325
         x=list(alternating=1,relation="free"),#, lim=c(0,70)),
326
         y=list(at=c(-18,-10,seq(0,100,len=5)),
327
           labels=c("MCD12Q1 IGBP","MOD35 path",seq(0,100,len=5)),
328
           lim=c(-25,100)),
329
         alternating=F),
330
       xlab="Distance Along Transect (km)", ylab="% Missing Data / % of Maximum Value",
331
       legend=list(
332
         bottom=list(fun=draw.key(list( rep=FALSE,columns=1,title=" ",
333
                      lines=list(type=c("b","b","b","b","b","l","b","l"),pch=16,cex=.5,
334
                        lty=c(0,1,1,1,1,5,1,5),
335
                        col=c("transparent","red","blue","black","darkgreen","darkgreen","orange","orange")),
336
                       text=list(
337
                         c("MODIS Products","C5 MOD09 % Cloudy","C5 MOD35 % Cloudy","C6 MOD35 % Cloudy","MOD17 % Missing","MOD17 (scaled)","MOD11 % Missing","MOD11 (scaled)")),
338
                       rectangles=list(border=NA,col=c(NA,"tan","darkgreen")),
339
                       text=list(c("C5 MOD35 Processing Path","Desert","Land")),
340
                       rectangles=list(border=NA,col=c(NA,IGBP$col[sort(unique(transd$value[transd$variable=="MCD12Q1"]+1))])),
341
                       text=list(c("MCD12Q1 IGBP Land Cover",IGBP$class[sort(unique(transd$value[transd$variable=="MCD12Q1"]+1))])))))),
342
  strip = strip.custom(par.strip.text=list(cex=.75)))
343
print(p4)
344
```
345

    
346
Compile the PDF:
347
```{r}
348
CairoPDF("output/mod35compare.pdf",width=11,height=7)
349
### Global Comparison
350
print(g1,position=c(0,.35,1,1),more=T)
351
print(g2,position=c(0,0,1,0.415),more=F)
352
         
353
### MOD35 Desert Processing path
354
levelplot(pp,asp=1,scales=list(draw=T,rot=0),maxpixels=1e6,
355
          at=c(-1:3),col.regions=c("blue","cyan","tan","darkgreen"),margin=F,
356
          colorkey=list(space="bottom",title="MOD35 Processing Path",labels=list(labels=c("Water","Coast","Desert","Land"),at=0:4-.5)))+layer(sp.polygons(bbs,lwd=2))+layer(sp.lines(coast,lwd=.5))
357
### levelplot of regions
358
print(p1,position=c(0,0,.62,1),more=T)
359
print(p2,position=c(0.6,0.21,0.78,0.79),more=T)
360
print(p3,position=c(0.76,0.21,1,0.79))
361
### profile plots
362
print(p4)
363
dev.off()
364
```
365

    
366
Derive summary statistics for manuscript
367
```{r}
368
td=cast(transect+loc+dist~variable,value="value",data=transd)
369
td2=melt.data.frame(td,id.vars=c("transect","dist","loc","MOD35pp","MCD12Q1"))
370

    
371
## function to prettyprint mean/sd's
372
msd= function(x) paste(round(mean(x,na.rm=T),1),"% ±",round(sd(x,na.rm=T),1),sep="")
373

    
374
cast(td2,transect+variable~MOD35pp,value="value",fun=msd)
375
cast(td2,transect+variable~MOD35pp+MCD12Q1,value="value",fun=msd)
376
cast(td2,transect+variable~.,value="value",fun=msd)
377

    
378
cast(td2,transect+variable~.,value="value",fun=msd)
379

    
380
cast(td2,variable~MOD35pp,value="value",fun=msd)
381
cast(td2,variable~.,value="value",fun=msd)
382

    
383
td[td$transect=="Venezuela",]
384
```
385

    
386
Export regional areas as KML for inclusion on website
387
```{r}
388
library(plotKML)
389

    
390
kml_open("output/modiscloud.kml")
391

    
392
readAll(mod35c5)
393

    
394
kml_layer.Raster(mod35c5,
395
     plot.legend = TRUE,raster_name="Collection 5 MOD35 Cloud Frequency",
396
    z.lim = c(0,100),colour_scale = get("colour_scale_numeric", envir = plotKML.opts),
397
#    home_url = get("home_url", envir = plotKML.opts),
398
#    metadata = NULL, html.table = NULL,
399
    altitudeMode = "clampToGround", balloon = FALSE
400
)
401

    
402
system(paste("gdal_translate -of KMLSUPEROVERLAY ",mod35c5@file@name," output/mod35c5.kmz -co FORMAT=JPEG"))
403

    
404
logo = "http://static.tumblr.com/t0afs9f/KWTm94tpm/yale_logo.png"
405
kml_screen(image.file = logo, position = "UL", sname = "YALE logo",size=c(.1,.1))
406
kml_close("modiscloud.kml")
407
kml_compress("modiscloud.kml",files=c(paste(month.name,".png",sep=""),"obj_legend.png"),zip="/usr/bin/zip")
408
```
(43-43/51)