|
1
|
<!DOCTYPE html>
|
|
2
|
<!-- saved from url=(0014)about:internet -->
|
|
3
|
<html>
|
|
4
|
<head>
|
|
5
|
<meta http-equiv="Content-Type" content="text/html; charset=utf-8"/>
|
|
6
|
|
|
7
|
<title>Google Earth Engine Processing</title>
|
|
8
|
|
|
9
|
<style type="text/css">
|
|
10
|
body, td {
|
|
11
|
font-family: sans-serif;
|
|
12
|
background-color: white;
|
|
13
|
font-size: 12px;
|
|
14
|
margin: 8px;
|
|
15
|
}
|
|
16
|
|
|
17
|
tt, code, pre {
|
|
18
|
font-family: 'DejaVu Sans Mono', 'Droid Sans Mono', 'Lucida Console', Consolas, Monaco, monospace;
|
|
19
|
}
|
|
20
|
|
|
21
|
h1 {
|
|
22
|
font-size:2.2em;
|
|
23
|
}
|
|
24
|
|
|
25
|
h2 {
|
|
26
|
font-size:1.8em;
|
|
27
|
}
|
|
28
|
|
|
29
|
h3 {
|
|
30
|
font-size:1.4em;
|
|
31
|
}
|
|
32
|
|
|
33
|
h4 {
|
|
34
|
font-size:1.0em;
|
|
35
|
}
|
|
36
|
|
|
37
|
h5 {
|
|
38
|
font-size:0.9em;
|
|
39
|
}
|
|
40
|
|
|
41
|
h6 {
|
|
42
|
font-size:0.8em;
|
|
43
|
}
|
|
44
|
|
|
45
|
a:visited {
|
|
46
|
color: rgb(50%, 0%, 50%);
|
|
47
|
}
|
|
48
|
|
|
49
|
pre {
|
|
50
|
margin-top: 0;
|
|
51
|
max-width: 95%;
|
|
52
|
border: 1px solid #ccc;
|
|
53
|
white-space: pre-wrap;
|
|
54
|
}
|
|
55
|
|
|
56
|
pre code {
|
|
57
|
display: block; padding: 0.5em;
|
|
58
|
}
|
|
59
|
|
|
60
|
code.r, code.cpp {
|
|
61
|
background-color: #F8F8F8;
|
|
62
|
}
|
|
63
|
|
|
64
|
table, td, th {
|
|
65
|
border: none;
|
|
66
|
}
|
|
67
|
|
|
68
|
blockquote {
|
|
69
|
color:#666666;
|
|
70
|
margin:0;
|
|
71
|
padding-left: 1em;
|
|
72
|
border-left: 0.5em #EEE solid;
|
|
73
|
}
|
|
74
|
|
|
75
|
hr {
|
|
76
|
height: 0px;
|
|
77
|
border-bottom: none;
|
|
78
|
border-top-width: thin;
|
|
79
|
border-top-style: dotted;
|
|
80
|
border-top-color: #999999;
|
|
81
|
}
|
|
82
|
|
|
83
|
@media print {
|
|
84
|
* {
|
|
85
|
background: transparent !important;
|
|
86
|
color: black !important;
|
|
87
|
filter:none !important;
|
|
88
|
-ms-filter: none !important;
|
|
89
|
}
|
|
90
|
|
|
91
|
body {
|
|
92
|
font-size:12pt;
|
|
93
|
max-width:100%;
|
|
94
|
}
|
|
95
|
|
|
96
|
a, a:visited {
|
|
97
|
text-decoration: underline;
|
|
98
|
}
|
|
99
|
|
|
100
|
hr {
|
|
101
|
visibility: hidden;
|
|
102
|
page-break-before: always;
|
|
103
|
}
|
|
104
|
|
|
105
|
pre, blockquote {
|
|
106
|
padding-right: 1em;
|
|
107
|
page-break-inside: avoid;
|
|
108
|
}
|
|
109
|
|
|
110
|
tr, img {
|
|
111
|
page-break-inside: avoid;
|
|
112
|
}
|
|
113
|
|
|
114
|
img {
|
|
115
|
max-width: 100% !important;
|
|
116
|
}
|
|
117
|
|
|
118
|
@page :left {
|
|
119
|
margin: 15mm 20mm 15mm 10mm;
|
|
120
|
}
|
|
121
|
|
|
122
|
@page :right {
|
|
123
|
margin: 15mm 10mm 15mm 20mm;
|
|
124
|
}
|
|
125
|
|
|
126
|
p, h2, h3 {
|
|
127
|
orphans: 3; widows: 3;
|
|
128
|
}
|
|
129
|
|
|
130
|
h2, h3 {
|
|
131
|
page-break-after: avoid;
|
|
132
|
}
|
|
133
|
}
|
|
134
|
|
|
135
|
</style>
|
|
136
|
|
|
137
|
<!-- Styles for R syntax highlighter -->
|
|
138
|
<style type="text/css">
|
|
139
|
pre .operator,
|
|
140
|
pre .paren {
|
|
141
|
color: rgb(104, 118, 135)
|
|
142
|
}
|
|
143
|
|
|
144
|
pre .literal {
|
|
145
|
color: rgb(88, 72, 246)
|
|
146
|
}
|
|
147
|
|
|
148
|
pre .number {
|
|
149
|
color: rgb(0, 0, 205);
|
|
150
|
}
|
|
151
|
|
|
152
|
pre .comment {
|
|
153
|
color: rgb(76, 136, 107);
|
|
154
|
}
|
|
155
|
|
|
156
|
pre .keyword {
|
|
157
|
color: rgb(0, 0, 255);
|
|
158
|
}
|
|
159
|
|
|
160
|
pre .identifier {
|
|
161
|
color: rgb(0, 0, 0);
|
|
162
|
}
|
|
163
|
|
|
164
|
pre .string {
|
|
165
|
color: rgb(3, 106, 7);
|
|
166
|
}
|
|
167
|
</style>
|
|
168
|
|
|
169
|
<!-- R syntax highlighter -->
|
|
170
|
<script type="text/javascript">
|
|
171
|
var hljs=new function(){function m(p){return p.replace(/&/gm,"&").replace(/</gm,"<")}function f(r,q,p){return RegExp(q,"m"+(r.cI?"i":"")+(p?"g":""))}function b(r){for(var p=0;p<r.childNodes.length;p++){var q=r.childNodes[p];if(q.nodeName=="CODE"){return q}if(!(q.nodeType==3&&q.nodeValue.match(/\s+/))){break}}}function h(t,s){var p="";for(var r=0;r<t.childNodes.length;r++){if(t.childNodes[r].nodeType==3){var q=t.childNodes[r].nodeValue;if(s){q=q.replace(/\n/g,"")}p+=q}else{if(t.childNodes[r].nodeName=="BR"){p+="\n"}else{p+=h(t.childNodes[r])}}}if(/MSIE [678]/.test(navigator.userAgent)){p=p.replace(/\r/g,"\n")}return p}function a(s){var r=s.className.split(/\s+/);r=r.concat(s.parentNode.className.split(/\s+/));for(var q=0;q<r.length;q++){var p=r[q].replace(/^language-/,"");if(e[p]){return p}}}function c(q){var p=[];(function(s,t){for(var r=0;r<s.childNodes.length;r++){if(s.childNodes[r].nodeType==3){t+=s.childNodes[r].nodeValue.length}else{if(s.childNodes[r].nodeName=="BR"){t+=1}else{if(s.childNodes[r].nodeType==1){p.push({event:"start",offset:t,node:s.childNodes[r]});t=arguments.callee(s.childNodes[r],t);p.push({event:"stop",offset:t,node:s.childNodes[r]})}}}}return t})(q,0);return p}function k(y,w,x){var q=0;var z="";var s=[];function u(){if(y.length&&w.length){if(y[0].offset!=w[0].offset){return(y[0].offset<w[0].offset)?y:w}else{return w[0].event=="start"?y:w}}else{return y.length?y:w}}function t(D){var A="<"+D.nodeName.toLowerCase();for(var B=0;B<D.attributes.length;B++){var C=D.attributes[B];A+=" "+C.nodeName.toLowerCase();if(C.value!==undefined&&C.value!==false&&C.value!==null){A+='="'+m(C.value)+'"'}}return A+">"}while(y.length||w.length){var v=u().splice(0,1)[0];z+=m(x.substr(q,v.offset-q));q=v.offset;if(v.event=="start"){z+=t(v.node);s.push(v.node)}else{if(v.event=="stop"){var p,r=s.length;do{r--;p=s[r];z+=("</"+p.nodeName.toLowerCase()+">")}while(p!=v.node);s.splice(r,1);while(r<s.length){z+=t(s[r]);r++}}}}return z+m(x.substr(q))}function j(){function q(x,y,v){if(x.compiled){return}var u;var s=[];if(x.k){x.lR=f(y,x.l||hljs.IR,true);for(var w in x.k){if(!x.k.hasOwnProperty(w)){continue}if(x.k[w] instanceof Object){u=x.k[w]}else{u=x.k;w="keyword"}for(var r in u){if(!u.hasOwnProperty(r)){continue}x.k[r]=[w,u[r]];s.push(r)}}}if(!v){if(x.bWK){x.b="\\b("+s.join("|")+")\\s"}x.bR=f(y,x.b?x.b:"\\B|\\b");if(!x.e&&!x.eW){x.e="\\B|\\b"}if(x.e){x.eR=f(y,x.e)}}if(x.i){x.iR=f(y,x.i)}if(x.r===undefined){x.r=1}if(!x.c){x.c=[]}x.compiled=true;for(var t=0;t<x.c.length;t++){if(x.c[t]=="self"){x.c[t]=x}q(x.c[t],y,false)}if(x.starts){q(x.starts,y,false)}}for(var p in e){if(!e.hasOwnProperty(p)){continue}q(e[p].dM,e[p],true)}}function d(B,C){if(!j.called){j();j.called=true}function q(r,M){for(var L=0;L<M.c.length;L++){if((M.c[L].bR.exec(r)||[null])[0]==r){return M.c[L]}}}function v(L,r){if(D[L].e&&D[L].eR.test(r)){return 1}if(D[L].eW){var M=v(L-1,r);return M?M+1:0}return 0}function w(r,L){return L.i&&L.iR.test(r)}function K(N,O){var M=[];for(var L=0;L<N.c.length;L++){M.push(N.c[L].b)}var r=D.length-1;do{if(D[r].e){M.push(D[r].e)}r--}while(D[r+1].eW);if(N.i){M.push(N.i)}return f(O,M.join("|"),true)}function p(M,L){var N=D[D.length-1];if(!N.t){N.t=K(N,E)}N.t.lastIndex=L;var r=N.t.exec(M);return r?[M.substr(L,r.index-L),r[0],false]:[M.substr(L),"",true]}function z(N,r){var L=E.cI?r[0].toLowerCase():r[0];var M=N.k[L];if(M&&M instanceof Array){return M}return false}function F(L,P){L=m(L);if(!P.k){return L}var r="";var O=0;P.lR.lastIndex=0;var M=P.lR.exec(L);while(M){r+=L.substr(O,M.index-O);var N=z(P,M);if(N){x+=N[1];r+='<span class="'+N[0]+'">'+M[0]+"</span>"}else{r+=M[0]}O=P.lR.lastIndex;M=P.lR.exec(L)}return r+L.substr(O,L.length-O)}function J(L,M){if(M.sL&&e[M.sL]){var r=d(M.sL,L);x+=r.keyword_count;return r.value}else{return F(L,M)}}function I(M,r){var L=M.cN?'<span class="'+M.cN+'">':"";if(M.rB){y+=L;M.buffer=""}else{if(M.eB){y+=m(r)+L;M.buffer=""}else{y+=L;M.buffer=r}}D.push(M);A+=M.r}function G(N,M,Q){var R=D[D.length-1];if(Q){y+=J(R.buffer+N,R);return false}var P=q(M,R);if(P){y+=J(R.buffer+N,R);I(P,M);return P.rB}var L=v(D.length-1,M);if(L){var O=R.cN?"</span>":"";if(R.rE){y+=J(R.buffer+N,R)+O}else{if(R.eE){y+=J(R.buffer+N,R)+O+m(M)}else{y+=J(R.buffer+N+M,R)+O}}while(L>1){O=D[D.length-2].cN?"</span>":"";y+=O;L--;D.length--}var r=D[D.length-1];D.length--;D[D.length-1].buffer="";if(r.starts){I(r.starts,"")}return R.rE}if(w(M,R)){throw"Illegal"}}var E=e[B];var D=[E.dM];var A=0;var x=0;var y="";try{var s,u=0;E.dM.buffer="";do{s=p(C,u);var t=G(s[0],s[1],s[2]);u+=s[0].length;if(!t){u+=s[1].length}}while(!s[2]);if(D.length>1){throw"Illegal"}return{r:A,keyword_count:x,value:y}}catch(H){if(H=="Illegal"){return{r:0,keyword_count:0,value:m(C)}}else{throw H}}}function g(t){var p={keyword_count:0,r:0,value:m(t)};var r=p;for(var q in e){if(!e.hasOwnProperty(q)){continue}var s=d(q,t);s.language=q;if(s.keyword_count+s.r>r.keyword_count+r.r){r=s}if(s.keyword_count+s.r>p.keyword_count+p.r){r=p;p=s}}if(r.language){p.second_best=r}return p}function i(r,q,p){if(q){r=r.replace(/^((<[^>]+>|\t)+)/gm,function(t,w,v,u){return w.replace(/\t/g,q)})}if(p){r=r.replace(/\n/g,"<br>")}return r}function n(t,w,r){var x=h(t,r);var v=a(t);var y,s;if(v){y=d(v,x)}else{return}var q=c(t);if(q.length){s=document.createElement("pre");s.innerHTML=y.value;y.value=k(q,c(s),x)}y.value=i(y.value,w,r);var u=t.className;if(!u.match("(\\s|^)(language-)?"+v+"(\\s|$)")){u=u?(u+" "+v):v}if(/MSIE [678]/.test(navigator.userAgent)&&t.tagName=="CODE"&&t.parentNode.tagName=="PRE"){s=t.parentNode;var p=document.createElement("div");p.innerHTML="<pre><code>"+y.value+"</code></pre>";t=p.firstChild.firstChild;p.firstChild.cN=s.cN;s.parentNode.replaceChild(p.firstChild,s)}else{t.innerHTML=y.value}t.className=u;t.result={language:v,kw:y.keyword_count,re:y.r};if(y.second_best){t.second_best={language:y.second_best.language,kw:y.second_best.keyword_count,re:y.second_best.r}}}function o(){if(o.called){return}o.called=true;var r=document.getElementsByTagName("pre");for(var p=0;p<r.length;p++){var q=b(r[p]);if(q){n(q,hljs.tabReplace)}}}function l(){if(window.addEventListener){window.addEventListener("DOMContentLoaded",o,false);window.addEventListener("load",o,false)}else{if(window.attachEvent){window.attachEvent("onload",o)}else{window.onload=o}}}var e={};this.LANGUAGES=e;this.highlight=d;this.highlightAuto=g;this.fixMarkup=i;this.highlightBlock=n;this.initHighlighting=o;this.initHighlightingOnLoad=l;this.IR="[a-zA-Z][a-zA-Z0-9_]*";this.UIR="[a-zA-Z_][a-zA-Z0-9_]*";this.NR="\\b\\d+(\\.\\d+)?";this.CNR="\\b(0[xX][a-fA-F0-9]+|(\\d+(\\.\\d*)?|\\.\\d+)([eE][-+]?\\d+)?)";this.BNR="\\b(0b[01]+)";this.RSR="!|!=|!==|%|%=|&|&&|&=|\\*|\\*=|\\+|\\+=|,|\\.|-|-=|/|/=|:|;|<|<<|<<=|<=|=|==|===|>|>=|>>|>>=|>>>|>>>=|\\?|\\[|\\{|\\(|\\^|\\^=|\\||\\|=|\\|\\||~";this.ER="(?![\\s\\S])";this.BE={b:"\\\\.",r:0};this.ASM={cN:"string",b:"'",e:"'",i:"\\n",c:[this.BE],r:0};this.QSM={cN:"string",b:'"',e:'"',i:"\\n",c:[this.BE],r:0};this.CLCM={cN:"comment",b:"//",e:"$"};this.CBLCLM={cN:"comment",b:"/\\*",e:"\\*/"};this.HCM={cN:"comment",b:"#",e:"$"};this.NM={cN:"number",b:this.NR,r:0};this.CNM={cN:"number",b:this.CNR,r:0};this.BNM={cN:"number",b:this.BNR,r:0};this.inherit=function(r,s){var p={};for(var q in r){p[q]=r[q]}if(s){for(var q in s){p[q]=s[q]}}return p}}();hljs.LANGUAGES.cpp=function(){var a={keyword:{"false":1,"int":1,"float":1,"while":1,"private":1,"char":1,"catch":1,"export":1,virtual:1,operator:2,sizeof:2,dynamic_cast:2,typedef:2,const_cast:2,"const":1,struct:1,"for":1,static_cast:2,union:1,namespace:1,unsigned:1,"long":1,"throw":1,"volatile":2,"static":1,"protected":1,bool:1,template:1,mutable:1,"if":1,"public":1,friend:2,"do":1,"return":1,"goto":1,auto:1,"void":2,"enum":1,"else":1,"break":1,"new":1,extern:1,using:1,"true":1,"class":1,asm:1,"case":1,typeid:1,"short":1,reinterpret_cast:2,"default":1,"double":1,register:1,explicit:1,signed:1,typename:1,"try":1,"this":1,"switch":1,"continue":1,wchar_t:1,inline:1,"delete":1,alignof:1,char16_t:1,char32_t:1,constexpr:1,decltype:1,noexcept:1,nullptr:1,static_assert:1,thread_local:1,restrict:1,_Bool:1,complex:1},built_in:{std:1,string:1,cin:1,cout:1,cerr:1,clog:1,stringstream:1,istringstream:1,ostringstream:1,auto_ptr:1,deque:1,list:1,queue:1,stack:1,vector:1,map:1,set:1,bitset:1,multiset:1,multimap:1,unordered_set:1,unordered_map:1,unordered_multiset:1,unordered_multimap:1,array:1,shared_ptr:1}};return{dM:{k:a,i:"</",c:[hljs.CLCM,hljs.CBLCLM,hljs.QSM,{cN:"string",b:"'\\\\?.",e:"'",i:"."},{cN:"number",b:"\\b(\\d+(\\.\\d*)?|\\.\\d+)(u|U|l|L|ul|UL|f|F)"},hljs.CNM,{cN:"preprocessor",b:"#",e:"$"},{cN:"stl_container",b:"\\b(deque|list|queue|stack|vector|map|set|bitset|multiset|multimap|unordered_map|unordered_set|unordered_multiset|unordered_multimap|array)\\s*<",e:">",k:a,r:10,c:["self"]}]}}}();hljs.LANGUAGES.r={dM:{c:[hljs.HCM,{cN:"number",b:"\\b0[xX][0-9a-fA-F]+[Li]?\\b",e:hljs.IMMEDIATE_RE,r:0},{cN:"number",b:"\\b\\d+(?:[eE][+\\-]?\\d*)?L\\b",e:hljs.IMMEDIATE_RE,r:0},{cN:"number",b:"\\b\\d+\\.(?!\\d)(?:i\\b)?",e:hljs.IMMEDIATE_RE,r:1},{cN:"number",b:"\\b\\d+(?:\\.\\d*)?(?:[eE][+\\-]?\\d*)?i?\\b",e:hljs.IMMEDIATE_RE,r:0},{cN:"number",b:"\\.\\d+(?:[eE][+\\-]?\\d*)?i?\\b",e:hljs.IMMEDIATE_RE,r:1},{cN:"keyword",b:"(?:tryCatch|library|setGeneric|setGroupGeneric)\\b",e:hljs.IMMEDIATE_RE,r:10},{cN:"keyword",b:"\\.\\.\\.",e:hljs.IMMEDIATE_RE,r:10},{cN:"keyword",b:"\\.\\.\\d+(?![\\w.])",e:hljs.IMMEDIATE_RE,r:10},{cN:"keyword",b:"\\b(?:function)",e:hljs.IMMEDIATE_RE,r:2},{cN:"keyword",b:"(?:if|in|break|next|repeat|else|for|return|switch|while|try|stop|warning|require|attach|detach|source|setMethod|setClass)\\b",e:hljs.IMMEDIATE_RE,r:1},{cN:"literal",b:"(?:NA|NA_integer_|NA_real_|NA_character_|NA_complex_)\\b",e:hljs.IMMEDIATE_RE,r:10},{cN:"literal",b:"(?:NULL|TRUE|FALSE|T|F|Inf|NaN)\\b",e:hljs.IMMEDIATE_RE,r:1},{cN:"identifier",b:"[a-zA-Z.][a-zA-Z0-9._]*\\b",e:hljs.IMMEDIATE_RE,r:0},{cN:"operator",b:"<\\-(?!\\s*\\d)",e:hljs.IMMEDIATE_RE,r:2},{cN:"operator",b:"\\->|<\\-",e:hljs.IMMEDIATE_RE,r:1},{cN:"operator",b:"%%|~",e:hljs.IMMEDIATE_RE},{cN:"operator",b:">=|<=|==|!=|\\|\\||&&|=|\\+|\\-|\\*|/|\\^|>|<|!|&|\\||\\$|:",e:hljs.IMMEDIATE_RE,r:0},{cN:"operator",b:"%",e:"%",i:"\\n",r:1},{cN:"identifier",b:"`",e:"`",r:0},{cN:"string",b:'"',e:'"',c:[hljs.BE],r:0},{cN:"string",b:"'",e:"'",c:[hljs.BE],r:0},{cN:"paren",b:"[[({\\])}]",e:hljs.IMMEDIATE_RE,r:0}]}};
|
|
172
|
hljs.initHighlightingOnLoad();
|
|
173
|
</script>
|
|
174
|
|
|
175
|
|
|
176
|
|
|
177
|
|
|
178
|
</head>
|
|
179
|
|
|
180
|
<body>
|
|
181
|
<p>Systematic landcover bias in Collection 5 MODIS cloud mask and derived products – a global overview</p>
|
|
182
|
|
|
183
|
<hr/>
|
|
184
|
|
|
185
|
<pre><code class="r">opts_chunk$set(eval = F)
|
|
186
|
</code></pre>
|
|
187
|
|
|
188
|
<p>This document describes the analysis of the Collection 5 MOD35 data.</p>
|
|
189
|
|
|
190
|
<h1>Google Earth Engine Processing</h1>
|
|
191
|
|
|
192
|
<p>The following code produces the annual (2009) summaries of cloud frequency from MOD09, MOD35, and MOD11 using the Google Earth Engine 'playground' API <a href="http://ee-api.appspot.com/">http://ee-api.appspot.com/</a>. </p>
|
|
193
|
|
|
194
|
<pre><code class="coffee">var startdate="2009-01-01"
|
|
195
|
var stopdate="2009-12-31"
|
|
196
|
|
|
197
|
// MOD11 MODIS LST
|
|
198
|
var mod11 = ee.ImageCollection("MOD11A2").map(function(img){
|
|
199
|
return img.select(['LST_Day_1km'])});
|
|
200
|
// MOD09 internal cloud flag
|
|
201
|
var mod09 = ee.ImageCollection("MOD09GA").filterDate(new Date(startdate),new Date(stopdate)).map(function(img) {
|
|
202
|
return img.select(['state_1km']).expression("((b(0)/1024)%2)");
|
|
203
|
});
|
|
204
|
// MOD35 cloud flag
|
|
205
|
var mod35 = ee.ImageCollection("MOD09GA").filterDate(new Date(startdate),new Date(stopdate)).map(function(img) {
|
|
206
|
return img.select(['state_1km']).expression("((b(0))%4)==1|((b(0))%4)==2");
|
|
207
|
});
|
|
208
|
|
|
209
|
//define reducers
|
|
210
|
var COUNT = ee.call("Reducer.count");
|
|
211
|
var MEAN = ee.call("Reducer.mean");
|
|
212
|
|
|
213
|
//a few maps of constants
|
|
214
|
c100=ee.Image(100); //to multiply by 100
|
|
215
|
c02=ee.Image(0.02); //to scale LST data
|
|
216
|
c272=ee.Image(272.15); // to convert K->C
|
|
217
|
|
|
218
|
//calculate mean cloudiness (%), rename, and convert to integer
|
|
219
|
mod09a=mod09.reduce(MEAN).select([0], ['MOD09']).multiply(c100).int8();
|
|
220
|
mod35a=mod35.reduce(MEAN).select([0], ['MOD35']).multiply(c100).int8();
|
|
221
|
|
|
222
|
/////////////////////////////////////////////////
|
|
223
|
// Generate the cloud frequency surface:
|
|
224
|
getMiss = function(collection) {
|
|
225
|
//filter by date
|
|
226
|
i2=collection.filterDate(new Date(startdate),new Date(stopdate));
|
|
227
|
// number of layers in collection
|
|
228
|
i2_n=i2.getInfo().features.length;
|
|
229
|
//get means
|
|
230
|
// image of -1s to convert to % missing
|
|
231
|
c1=ee.Image(-1);
|
|
232
|
// 1 Calculate the number of days with measurements
|
|
233
|
// 2 divide by the total number of layers
|
|
234
|
i2_c=ee.Image(i2_n).float()
|
|
235
|
// 3 Add -1 and multiply by -1 to invert to % cloudy
|
|
236
|
// 4 Rename to "Percent_Cloudy"
|
|
237
|
// 5 multiply by 100 and convert to 8-bit integer to decrease file size
|
|
238
|
i2_miss=i2.reduce(COUNT).divide(i2_c).add(c1).multiply(c1).multiply(c100).int8();
|
|
239
|
return (i2_miss);
|
|
240
|
};
|
|
241
|
|
|
242
|
// run the function
|
|
243
|
mod11a=getMiss(mod11).select([0], ['MOD11_LST_PMiss']);
|
|
244
|
|
|
245
|
// get long-term mean
|
|
246
|
mod11b=mod11.reduce(MEAN).multiply(c02).subtract(c272).int8().select([0], ['MOD11_LST_MEAN']);
|
|
247
|
|
|
248
|
// summary object with all layers
|
|
249
|
summary=mod11a.addBands(mod11b).addBands(mod35a).addBands(mod09a)
|
|
250
|
|
|
251
|
var region='[[-180, -60], [-180, 90], [180, 90], [180, -60]]' //global
|
|
252
|
|
|
253
|
// get download link
|
|
254
|
print("All")
|
|
255
|
var path = summary.getDownloadURL({
|
|
256
|
'scale': 1000,
|
|
257
|
'crs': 'EPSG:4326',
|
|
258
|
'region': region
|
|
259
|
});
|
|
260
|
print('https://earthengine.sandbox.google.com' + path);
|
|
261
|
</code></pre>
|
|
262
|
|
|
263
|
<h1>Data Processing</h1>
|
|
264
|
|
|
265
|
<pre><code class="r">setwd("~/acrobates/adamw/projects/MOD35C5")
|
|
266
|
library(raster)
|
|
267
|
</code></pre>
|
|
268
|
|
|
269
|
<pre><code>## Loading required package: sp
|
|
270
|
</code></pre>
|
|
271
|
|
|
272
|
<pre><code class="r">beginCluster(10)
|
|
273
|
</code></pre>
|
|
274
|
|
|
275
|
<pre><code>## Loading required package: snow
|
|
276
|
</code></pre>
|
|
277
|
|
|
278
|
<pre><code class="r">library(rasterVis)
|
|
279
|
</code></pre>
|
|
280
|
|
|
281
|
<pre><code>## Loading required package: lattice Loading required package: latticeExtra
|
|
282
|
## Loading required package: RColorBrewer Loading required package: hexbin
|
|
283
|
## Loading required package: grid
|
|
284
|
</code></pre>
|
|
285
|
|
|
286
|
<pre><code class="r">library(rgdal)
|
|
287
|
</code></pre>
|
|
288
|
|
|
289
|
<pre><code>## rgdal: version: 0.8-10, (SVN revision 478) Geospatial Data Abstraction
|
|
290
|
## Library extensions to R successfully loaded Loaded GDAL runtime: GDAL
|
|
291
|
## 1.9.2, released 2012/10/08 but rgdal build and GDAL runtime not in sync:
|
|
292
|
## ... consider re-installing rgdal!! Path to GDAL shared files:
|
|
293
|
## /usr/share/gdal/1.9 Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012,
|
|
294
|
## [PJ_VERSION: 480] Path to PROJ.4 shared files: (autodetected)
|
|
295
|
</code></pre>
|
|
296
|
|
|
297
|
<pre><code class="r">library(plotKML)
|
|
298
|
</code></pre>
|
|
299
|
|
|
300
|
<pre><code>## plotKML version 0.3-5 (2013-05-16) URL:
|
|
301
|
## http://plotkml.r-forge.r-project.org/
|
|
302
|
##
|
|
303
|
## Attaching package: 'plotKML'
|
|
304
|
##
|
|
305
|
## The following object is masked from 'package:raster':
|
|
306
|
##
|
|
307
|
## count
|
|
308
|
</code></pre>
|
|
309
|
|
|
310
|
<pre><code class="r">library(Cairo)
|
|
311
|
library(reshape)
|
|
312
|
</code></pre>
|
|
313
|
|
|
314
|
<pre><code>## Loading required package: plyr
|
|
315
|
##
|
|
316
|
## Attaching package: 'plyr'
|
|
317
|
##
|
|
318
|
## The following object is masked from 'package:plotKML':
|
|
319
|
##
|
|
320
|
## count
|
|
321
|
##
|
|
322
|
## The following object is masked from 'package:raster':
|
|
323
|
##
|
|
324
|
## count
|
|
325
|
##
|
|
326
|
## Attaching package: 'reshape'
|
|
327
|
##
|
|
328
|
## The following object is masked from 'package:plyr':
|
|
329
|
##
|
|
330
|
## rename, round_any
|
|
331
|
##
|
|
332
|
## The following object is masked from 'package:raster':
|
|
333
|
##
|
|
334
|
## expand
|
|
335
|
</code></pre>
|
|
336
|
|
|
337
|
<pre><code class="r">library(rgeos)
|
|
338
|
</code></pre>
|
|
339
|
|
|
340
|
<pre><code>## rgeos version: 0.2-19, (SVN revision 394) GEOS runtime version:
|
|
341
|
## 3.3.3-CAPI-1.7.4 Polygon checking: TRUE
|
|
342
|
</code></pre>
|
|
343
|
|
|
344
|
<pre><code class="r">library(splancs)
|
|
345
|
</code></pre>
|
|
346
|
|
|
347
|
<pre><code>## Spatial Point Pattern Analysis Code in S-Plus
|
|
348
|
##
|
|
349
|
## Version 2 - Spatial and Space-Time analysis
|
|
350
|
##
|
|
351
|
## Attaching package: 'splancs'
|
|
352
|
##
|
|
353
|
## The following object is masked from 'package:raster':
|
|
354
|
##
|
|
355
|
## zoom
|
|
356
|
</code></pre>
|
|
357
|
|
|
358
|
<pre><code class="r">
|
|
359
|
## get % cloudy
|
|
360
|
mod09 = raster("data/MOD09_2009.tif")
|
|
361
|
names(mod09) = "C5MOD09CF"
|
|
362
|
NAvalue(mod09) = 0
|
|
363
|
|
|
364
|
mod35c5 = raster("data/MOD35_2009.tif")
|
|
365
|
names(mod35c5) = "C5MOD35CF"
|
|
366
|
NAvalue(mod35c5) = 0
|
|
367
|
|
|
368
|
## mod35C6 annual
|
|
369
|
mod35c6 = raster("data/MOD35C6_2009.tif")
|
|
370
|
names(mod35c6) = "C6MOD35CF"
|
|
371
|
NAvalue(mod35c6) = 255
|
|
372
|
|
|
373
|
## landcover
|
|
374
|
lulc = raster("data/MCD12Q1_IGBP_2009_051_wgs84_1km.tif")
|
|
375
|
|
|
376
|
# lulc=ratify(lulc)
|
|
377
|
data(worldgrids_pal) #load palette
|
|
378
|
IGBP = data.frame(ID = 0:16, col = worldgrids_pal$IGBP[-c(18, 19)], lulc_levels2 = c("Water",
|
|
379
|
"Forest", "Forest", "Forest", "Forest", "Forest", "Shrublands", "Shrublands",
|
|
380
|
"Savannas", "Savannas", "Grasslands", "Permanent wetlands", "Croplands",
|
|
381
|
"Urban and built-up", "Cropland/Natural vegetation mosaic", "Snow and ice",
|
|
382
|
"Barren or sparsely vegetated"), stringsAsFactors = F)
|
|
383
|
IGBP$class = rownames(IGBP)
|
|
384
|
rownames(IGBP) = 1:nrow(IGBP)
|
|
385
|
levels(lulc) = list(IGBP)
|
|
386
|
names(lulc) = "MCD12Q1"
|
|
387
|
|
|
388
|
## MOD17
|
|
389
|
mod17 = raster("data/MOD17.tif", format = "GTiff")
|
|
390
|
NAvalue(mod17) = 65535
|
|
391
|
names(mod17) = "MOD17_unscaled"
|
|
392
|
|
|
393
|
mod17qc = raster("data/MOD17qc.tif", format = "GTiff")
|
|
394
|
NAvalue(mod17qc) = 255
|
|
395
|
names(mod17qc) = "MOD17CF"
|
|
396
|
|
|
397
|
## MOD11 via earth engine
|
|
398
|
mod11 = raster("data/MOD11_2009.tif", format = "GTiff")
|
|
399
|
names(mod11) = "MOD11_unscaled"
|
|
400
|
NAvalue(mod11) = 0
|
|
401
|
|
|
402
|
mod11qc = raster("data/MOD11qc_2009.tif", format = "GTiff")
|
|
403
|
names(mod11qc) = "MOD11CF"
|
|
404
|
</code></pre>
|
|
405
|
|
|
406
|
<p>Import the Collection 5 MOD35 processing path:</p>
|
|
407
|
|
|
408
|
<pre><code class="r">pp = raster("data/MOD35pp.tif")
|
|
409
|
NAvalue(pp) = 255
|
|
410
|
names(pp) = "MOD35pp"
|
|
411
|
</code></pre>
|
|
412
|
|
|
413
|
<p>Define transects to illustrate the fine-grain relationship between MOD35 cloud frequency and both landcover and processing path.</p>
|
|
414
|
|
|
415
|
<pre><code class="r">r1 = Lines(list(Line(matrix(c(-61.688, 4.098, -59.251, 3.43), ncol = 2, byrow = T))),
|
|
416
|
"Venezuela")
|
|
417
|
r2 = Lines(list(Line(matrix(c(133.746, -31.834, 134.226, -32.143), ncol = 2,
|
|
418
|
byrow = T))), "Australia")
|
|
419
|
r3 = Lines(list(Line(matrix(c(73.943, 27.419, 74.369, 26.877), ncol = 2, byrow = T))),
|
|
420
|
"India")
|
|
421
|
r4 = Lines(list(Line(matrix(c(33.195, 12.512, 33.802, 12.894), ncol = 2, byrow = T))),
|
|
422
|
"Sudan")
|
|
423
|
|
|
424
|
trans = SpatialLines(list(r1, r2, r3, r4), CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs "))
|
|
425
|
### write out shapefiles of transects
|
|
426
|
writeOGR(SpatialLinesDataFrame(trans, data = data.frame(ID = names(trans)),
|
|
427
|
match.ID = F), "output", layer = "transects", driver = "ESRI Shapefile",
|
|
428
|
overwrite = T)
|
|
429
|
</code></pre>
|
|
430
|
|
|
431
|
<p>Buffer transects to identify a small region around each transect for comparison and plotting</p>
|
|
432
|
|
|
433
|
<pre><code class="r">transb = gBuffer(trans, byid = T, width = 0.4)
|
|
434
|
## make polygons of bounding boxes
|
|
435
|
bb0 <- lapply(slot(transb, "polygons"), bbox)
|
|
436
|
bb1 <- lapply(bb0, bboxx)
|
|
437
|
# turn these into matrices using a helper function in splancs
|
|
438
|
bb2 <- lapply(bb1, function(x) rbind(x, x[1, ]))
|
|
439
|
# close the matrix rings by appending the first coordinate
|
|
440
|
rn <- row.names(transb)
|
|
441
|
# get the IDs
|
|
442
|
bb3 <- vector(mode = "list", length = length(bb2))
|
|
443
|
# make somewhere to keep the output
|
|
444
|
for (i in seq(along = bb3)) bb3[[i]] <- Polygons(list(Polygon(bb2[[i]])), ID = rn[i])
|
|
445
|
# loop over the closed matrix rings, adding the IDs
|
|
446
|
bbs <- SpatialPolygons(bb3, proj4string = CRS(proj4string(transb)))
|
|
447
|
</code></pre>
|
|
448
|
|
|
449
|
<p>Extract the CF and mean values from each raster of interest.</p>
|
|
450
|
|
|
451
|
<pre><code class="r">trd1 = lapply(1:length(transb), function(x) {
|
|
452
|
td = crop(mod11, transb[x])
|
|
453
|
tdd = lapply(list(mod35c5, mod35c6, mod09, mod17, mod17qc, mod11, mod11qc,
|
|
454
|
lulc, pp), function(l) resample(crop(l, transb[x]), td, method = "ngb"))
|
|
455
|
## normalize MOD11 and MOD17
|
|
456
|
for (j in which(do.call(c, lapply(tdd, function(i) names(i))) %in% c("MOD11_unscaled",
|
|
457
|
"MOD17_unscaled"))) {
|
|
458
|
trange = cellStats(tdd[[j]], range)
|
|
459
|
tscaled = 100 * (tdd[[j]] - trange[1])/(trange[2] - trange[1])
|
|
460
|
tscaled@history = list(range = trange)
|
|
461
|
names(tscaled) = sub("_unscaled", "", names(tdd[[j]]))
|
|
462
|
tdd = c(tdd, tscaled)
|
|
463
|
}
|
|
464
|
return(brick(tdd))
|
|
465
|
})
|
|
466
|
## bind all subregions into single dataframe for plotting
|
|
467
|
trd = do.call(rbind.data.frame, lapply(1:length(trd1), function(i) {
|
|
468
|
d = as.data.frame(as.matrix(trd1[[i]]))
|
|
469
|
d[, c("x", "y")] = coordinates(trd1[[i]])
|
|
470
|
d$trans = names(trans)[i]
|
|
471
|
d = melt(d, id.vars = c("trans", "x", "y"))
|
|
472
|
return(d)
|
|
473
|
}))
|
|
474
|
transd = do.call(rbind.data.frame, lapply(1:length(trans), function(l) {
|
|
475
|
td = as.data.frame(extract(trd1[[l]], trans[l], along = T, cellnumbers = F)[[1]])
|
|
476
|
td$loc = extract(trd1[[l]], trans[l], along = T, cellnumbers = T)[[1]][,
|
|
477
|
1]
|
|
478
|
td[, c("x", "y")] = xyFromCell(trd1[[l]], td$loc)
|
|
479
|
td$dist = spDistsN1(as.matrix(td[, c("x", "y")]), as.matrix(td[1, c("x",
|
|
480
|
"y")]), longlat = T)
|
|
481
|
td$transect = names(trans[l])
|
|
482
|
td2 = melt(td, id.vars = c("loc", "x", "y", "dist", "transect"))
|
|
483
|
td2 = td2[order(td2$variable, td2$dist), ]
|
|
484
|
# get per variable ranges to normalize
|
|
485
|
tr = cast(melt.list(tapply(td2$value, td2$variable, function(x) data.frame(min = min(x,
|
|
486
|
na.rm = T), max = max(x, na.rm = T)))), L1 ~ variable)
|
|
487
|
td2$min = tr$min[match(td2$variable, tr$L1)]
|
|
488
|
td2$max = tr$max[match(td2$variable, tr$L1)]
|
|
489
|
print(paste("Finished ", names(trans[l])))
|
|
490
|
return(td2)
|
|
491
|
}))
|
|
492
|
|
|
493
|
transd$type = ifelse(grepl("MOD35|MOD09|CF", transd$variable), "CF", "Data")
|
|
494
|
</code></pre>
|
|
495
|
|
|
496
|
<p>Compute difference between MOD09 and MOD35 cloud masks</p>
|
|
497
|
|
|
498
|
<pre><code class="r">## comparison of % cloudy days
|
|
499
|
dif_c5_09 = raster("data/dif_c5_09.tif", format = "GTiff")
|
|
500
|
</code></pre>
|
|
501
|
|
|
502
|
<p>Define a color scheme</p>
|
|
503
|
|
|
504
|
<pre><code class="r">n = 100
|
|
505
|
at = seq(0, 100, len = n)
|
|
506
|
bgyr = colorRampPalette(c("purple", "blue", "green", "yellow", "orange", "red",
|
|
507
|
"red"))
|
|
508
|
bgrayr = colorRampPalette(c("purple", "blue", "grey", "red", "red"))
|
|
509
|
cols = bgyr(n)
|
|
510
|
</code></pre>
|
|
511
|
|
|
512
|
<p>Import a global coastline map for overlay</p>
|
|
513
|
|
|
514
|
<pre><code class="r">library(maptools)
|
|
515
|
coast = map2SpatialLines(map("world", interior = FALSE, plot = FALSE), proj4string = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))
|
|
516
|
</code></pre>
|
|
517
|
|
|
518
|
<p>Draw the global cloud frequencies</p>
|
|
519
|
|
|
520
|
<pre><code class="r">g1 = levelplot(stack(mod35c5, mod09), xlab = " ", scales = list(x = list(draw = F),
|
|
521
|
y = list(alternating = 1)), col.regions = cols, at = at) + layer(sp.polygons(bbs[1:4],
|
|
522
|
lwd = 2)) + layer(sp.lines(coast, lwd = 0.5))
|
|
523
|
|
|
524
|
g2 = levelplot(dif_c5_09, col.regions = bgrayr(100), at = seq(-70, 70, len = 100),
|
|
525
|
margin = F, ylab = " ", colorkey = list("right")) + layer(sp.polygons(bbs[1:4],
|
|
526
|
lwd = 2)) + layer(sp.lines(coast, lwd = 0.5))
|
|
527
|
g2$strip = strip.custom(var.name = "Difference (C5MOD35-C5MOD09)", style = 1,
|
|
528
|
strip.names = T, strip.levels = F)
|
|
529
|
</code></pre>
|
|
530
|
|
|
531
|
<p>Now illustrate the fine-grain regions</p>
|
|
532
|
|
|
533
|
<pre><code class="r">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"),
|
|
534
|
at=at,col.regions=cols,maxpixels=7e6,
|
|
535
|
ylab="Latitude",xlab="Longitude"),strip = strip.custom(par.strip.text=list(cex=.7)))+layer(sp.lines(trans,lwd=2))
|
|
536
|
|
|
537
|
p2=useOuterStrips(
|
|
538
|
levelplot(value~x*y|variable+trans,data=trd[trd$variable%in%c("MCD12Q1"),],
|
|
539
|
asp=1,scales=list(draw=F,rot=0,relation="free"),colorkey=F,
|
|
540
|
at=c(-1,IGBP$ID),col.regions=IGBP$col,maxpixels=7e7,
|
|
541
|
legend=list(
|
|
542
|
right=list(fun=draw.key(list(columns=1,#title="MCD12Q1 \n IGBP Land \n Cover",
|
|
543
|
rectangles=list(col=IGBP$col,size=1),
|
|
544
|
text=list(as.character(IGBP$ID),at=IGBP$ID-.5))))),
|
|
545
|
ylab="",xlab=" "),strip = strip.custom(par.strip.text=list(cex=.7)),strip.left=F)+layer(sp.lines(trans,lwd=2))
|
|
546
|
p3=useOuterStrips(
|
|
547
|
levelplot(value~x*y|variable+trans,data=trd[trd$variable%in%c("MOD35pp"),],
|
|
548
|
asp=1,scales=list(draw=F,rot=0,relation="free"),colorkey=F,
|
|
549
|
at=c(-1:4),col.regions=c("blue","cyan","tan","darkgreen"),maxpixels=7e7,
|
|
550
|
legend=list(
|
|
551
|
right=list(fun=draw.key(list(columns=1,#title="MOD35 \n Processing \n Path",
|
|
552
|
rectangles=list(col=c("blue","cyan","tan","darkgreen"),size=1),
|
|
553
|
text=list(c("Water","Coast","Desert","Land")))))),
|
|
554
|
ylab="",xlab=" "),strip = strip.custom(par.strip.text=list(cex=.7)),strip.left=F)+layer(sp.lines(trans,lwd=2))
|
|
555
|
</code></pre>
|
|
556
|
|
|
557
|
<p>Now draw the profile plots for each transect.</p>
|
|
558
|
|
|
559
|
<pre><code class="r">## transects
|
|
560
|
p4=xyplot(value~dist|transect,groups=variable,type=c("smooth","p"),
|
|
561
|
data=transd,panel=function(...,subscripts=subscripts) {
|
|
562
|
td=transd[subscripts,]
|
|
563
|
## mod09
|
|
564
|
imod09=td$variable=="C5MOD09CF"
|
|
565
|
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)
|
|
566
|
## mod35C5
|
|
567
|
imod35=td$variable=="C5MOD35CF"
|
|
568
|
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)
|
|
569
|
## mod35C6
|
|
570
|
imod35c6=td$variable=="C6MOD35CF"
|
|
571
|
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)
|
|
572
|
## mod17
|
|
573
|
imod17=td$variable=="MOD17"
|
|
574
|
panel.xyplot(td$dist[imod17],100*((td$value[imod17]-td$min[imod17][1])/(td$max[imod17][1]-td$min[imod17][1])),
|
|
575
|
type=c("smooth"),span=0.09,subscripts=1:sum(imod17),col="darkgreen",lty=5,pch=1,cex=.25)
|
|
576
|
imod17qc=td$variable=="MOD17CF"
|
|
577
|
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)
|
|
578
|
## mod11
|
|
579
|
imod11=td$variable=="MOD11"
|
|
580
|
panel.xyplot(td$dist[imod11],100*((td$value[imod11]-td$min[imod11][1])/(td$max[imod11][1]-td$min[imod11][1])),
|
|
581
|
type=c("smooth"),span=0.09,subscripts=1:sum(imod17),col="orange",lty="dashed",pch=1,cex=.25)
|
|
582
|
imod11qc=td$variable=="MOD11CF"
|
|
583
|
qcspan=ifelse(td$transect[1]=="Australia",0.2,0.05)
|
|
584
|
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)
|
|
585
|
## land
|
|
586
|
path=td[td$variable=="MOD35pp",]
|
|
587
|
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")
|
|
588
|
land=td[td$variable=="MCD12Q1",]
|
|
589
|
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")
|
|
590
|
},subscripts=T,par.settings = list(grid.pars = list(lineend = "butt")),
|
|
591
|
scales=list(
|
|
592
|
x=list(alternating=1,relation="free"),#, lim=c(0,70)),
|
|
593
|
y=list(at=c(-18,-10,seq(0,100,len=5)),
|
|
594
|
labels=c("MCD12Q1 IGBP","MOD35 path",seq(0,100,len=5)),
|
|
595
|
lim=c(-25,100)),
|
|
596
|
alternating=F),
|
|
597
|
xlab="Distance Along Transect (km)", ylab="% Missing Data / % of Maximum Value",
|
|
598
|
legend=list(
|
|
599
|
bottom=list(fun=draw.key(list( rep=FALSE,columns=1,title=" ",
|
|
600
|
lines=list(type=c("b","b","b","b","b","l","b","l"),pch=16,cex=.5,
|
|
601
|
lty=c(0,1,1,1,1,5,1,5),
|
|
602
|
col=c("transparent","red","blue","black","darkgreen","darkgreen","orange","orange")),
|
|
603
|
text=list(
|
|
604
|
c("MODIS Products","C5 MOD09 % Cloudy","C5 MOD35 % Cloudy","C6 MOD35 % Cloudy","MOD17 % Missing","MOD17 (scaled)","MOD11 % Missing","MOD11 (scaled)")),
|
|
605
|
rectangles=list(border=NA,col=c(NA,"tan","darkgreen")),
|
|
606
|
text=list(c("C5 MOD35 Processing Path","Desert","Land")),
|
|
607
|
rectangles=list(border=NA,col=c(NA,IGBP$col[sort(unique(transd$value[transd$variable=="MCD12Q1"]+1))])),
|
|
608
|
text=list(c("MCD12Q1 IGBP Land Cover",IGBP$class[sort(unique(transd$value[transd$variable=="MCD12Q1"]+1))])))))),
|
|
609
|
strip = strip.custom(par.strip.text=list(cex=.75)))
|
|
610
|
print(p4)
|
|
611
|
</code></pre>
|
|
612
|
|
|
613
|
<p>Compile the PDF:</p>
|
|
614
|
|
|
615
|
<pre><code class="r">CairoPDF("output/mod35compare.pdf", width = 11, height = 7)
|
|
616
|
### Global Comparison
|
|
617
|
print(g1, position = c(0, 0.35, 1, 1), more = T)
|
|
618
|
print(g2, position = c(0, 0, 1, 0.415), more = F)
|
|
619
|
|
|
620
|
### MOD35 Desert Processing path
|
|
621
|
levelplot(pp, asp = 1, scales = list(draw = T, rot = 0), maxpixels = 1e+06,
|
|
622
|
at = c(-1:3), col.regions = c("blue", "cyan", "tan", "darkgreen"), margin = F,
|
|
623
|
colorkey = list(space = "bottom", title = "MOD35 Processing Path", labels = list(labels = c("Water",
|
|
624
|
"Coast", "Desert", "Land"), at = 0:4 - 0.5))) + layer(sp.polygons(bbs,
|
|
625
|
lwd = 2)) + layer(sp.lines(coast, lwd = 0.5))
|
|
626
|
### levelplot of regions
|
|
627
|
print(p1, position = c(0, 0, 0.62, 1), more = T)
|
|
628
|
print(p2, position = c(0.6, 0.21, 0.78, 0.79), more = T)
|
|
629
|
print(p3, position = c(0.76, 0.21, 1, 0.79))
|
|
630
|
### profile plots
|
|
631
|
print(p4)
|
|
632
|
dev.off()
|
|
633
|
</code></pre>
|
|
634
|
|
|
635
|
<p>Derive summary statistics for manuscript</p>
|
|
636
|
|
|
637
|
<pre><code class="r">td = cast(transect + loc + dist ~ variable, value = "value", data = transd)
|
|
638
|
td2 = melt.data.frame(td, id.vars = c("transect", "dist", "loc", "MOD35pp",
|
|
639
|
"MCD12Q1"))
|
|
640
|
|
|
641
|
## function to prettyprint mean/sd's
|
|
642
|
msd = function(x) paste(round(mean(x, na.rm = T), 1), "% ±", round(sd(x, na.rm = T),
|
|
643
|
1), sep = "")
|
|
644
|
|
|
645
|
cast(td2, transect + variable ~ MOD35pp, value = "value", fun = msd)
|
|
646
|
cast(td2, transect + variable ~ MOD35pp + MCD12Q1, value = "value", fun = msd)
|
|
647
|
cast(td2, transect + variable ~ ., value = "value", fun = msd)
|
|
648
|
|
|
649
|
cast(td2, transect + variable ~ ., value = "value", fun = msd)
|
|
650
|
|
|
651
|
cast(td2, variable ~ MOD35pp, value = "value", fun = msd)
|
|
652
|
cast(td2, variable ~ ., value = "value", fun = msd)
|
|
653
|
|
|
654
|
td[td$transect == "Venezuela", ]
|
|
655
|
</code></pre>
|
|
656
|
|
|
657
|
<p>Export regional areas as KML for inclusion on website</p>
|
|
658
|
|
|
659
|
<pre><code class="r">library(plotKML)
|
|
660
|
|
|
661
|
kml_open("output/modiscloud.kml")
|
|
662
|
|
|
663
|
readAll(mod35c5)
|
|
664
|
|
|
665
|
kml_layer.Raster(mod35c5,
|
|
666
|
plot.legend = TRUE,raster_name="Collection 5 MOD35 Cloud Frequency",
|
|
667
|
z.lim = c(0,100),colour_scale = get("colour_scale_numeric", envir = plotKML.opts),
|
|
668
|
# home_url = get("home_url", envir = plotKML.opts),
|
|
669
|
# metadata = NULL, html.table = NULL,
|
|
670
|
altitudeMode = "clampToGround", balloon = FALSE
|
|
671
|
)
|
|
672
|
|
|
673
|
system(paste("gdal_translate -of KMLSUPEROVERLAY ",mod35c5@file@name," output/mod35c5.kmz -co FORMAT=JPEG"))
|
|
674
|
|
|
675
|
logo = "http://static.tumblr.com/t0afs9f/KWTm94tpm/yale_logo.png"
|
|
676
|
kml_screen(image.file = logo, position = "UL", sname = "YALE logo",size=c(.1,.1))
|
|
677
|
kml_close("modiscloud.kml")
|
|
678
|
kml_compress("modiscloud.kml",files=c(paste(month.name,".png",sep=""),"obj_legend.png"),zip="/usr/bin/zip")
|
|
679
|
</code></pre>
|
|
680
|
|
|
681
|
</body>
|
|
682
|
|
|
683
|
</html>
|
|
684
|
|
