| 1 |
555815c9
|
Adam M. Wilson
|
<!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>
|
