Project

General

Profile

Download (23.1 KB) Statistics
| Branch: | Revision:
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
<meta http-equiv="x-ua-compatible" content="IE=9" >
7

    
8
<title>LST Temporal Aggregation Evaluation</title>
9

    
10
<style type="text/css">
11
body, td {
12
   font-family: sans-serif;
13
   background-color: white;
14
   font-size: 12px;
15
   margin: 8px;
16
}
17

    
18
tt, code, pre {
19
   font-family: 'DejaVu Sans Mono', 'Droid Sans Mono', 'Lucida Console', Consolas, Monaco, monospace;
20
}
21

    
22
h1 { 
23
   font-size:2.2em; 
24
}
25

    
26
h2 { 
27
   font-size:1.8em; 
28
}
29

    
30
h3 { 
31
   font-size:1.4em; 
32
}
33

    
34
h4 { 
35
   font-size:1.0em; 
36
}
37

    
38
h5 { 
39
   font-size:0.9em; 
40
}
41

    
42
h6 { 
43
   font-size:0.8em; 
44
}
45

    
46
a:visited {
47
   color: rgb(50%, 0%, 50%);
48
}
49

    
50
pre {	
51
   margin-top: 0;
52
   max-width: 95%;
53
   border: 1px solid #ccc;
54
   white-space: pre-wrap;
55
}
56

    
57
pre code {
58
   display: block; padding: 0.5em;
59
}
60

    
61
code.r, code.cpp {
62
   background-color: #F8F8F8;
63
}
64

    
65
table, td, th {
66
  border: none;
67
}
68

    
69
blockquote {
70
   color:#666666;
71
   margin:0;
72
   padding-left: 1em;
73
   border-left: 0.5em #EEE solid;
74
}
75

    
76
hr {
77
   height: 0px;
78
   border-bottom: none;
79
   border-top-width: thin;
80
   border-top-style: dotted;
81
   border-top-color: #999999;
82
}
83

    
84
@media print {
85
   * { 
86
      background: transparent !important; 
87
      color: black !important; 
88
      filter:none !important; 
89
      -ms-filter: none !important; 
90
   }
91

    
92
   body { 
93
      font-size:12pt; 
94
      max-width:100%; 
95
   }
96
       
97
   a, a:visited { 
98
      text-decoration: underline; 
99
   }
100

    
101
   hr { 
102
      visibility: hidden;
103
      page-break-before: always;
104
   }
105

    
106
   pre, blockquote { 
107
      padding-right: 1em; 
108
      page-break-inside: avoid; 
109
   }
110

    
111
   tr, img { 
112
      page-break-inside: avoid; 
113
   }
114

    
115
   img { 
116
      max-width: 100% !important; 
117
   }
118

    
119
   @page :left { 
120
      margin: 15mm 20mm 15mm 10mm; 
121
   }
122
     
123
   @page :right { 
124
      margin: 15mm 10mm 15mm 20mm; 
125
   }
126

    
127
   p, h2, h3 { 
128
      orphans: 3; widows: 3; 
129
   }
130

    
131
   h2, h3 { 
132
      page-break-after: avoid; 
133
   }
134
}
135

    
136
</style>
137

    
138
<!-- Styles for R syntax highlighter -->
139
<style type="text/css">
140
   pre .operator,
141
   pre .paren {
142
     color: rgb(104, 118, 135)
143
   }
144

    
145
   pre .literal {
146
     color: rgb(88, 72, 246)
147
   }
148

    
149
   pre .number {
150
     color: rgb(0, 0, 205);
151
   }
152

    
153
   pre .comment {
154
     color: rgb(76, 136, 107);
155
   }
156

    
157
   pre .keyword {
158
     color: rgb(0, 0, 255);
159
   }
160

    
161
   pre .identifier {
162
     color: rgb(0, 0, 0);
163
   }
164

    
165
   pre .string {
166
     color: rgb(3, 106, 7);
167
   }
168
</style>
169

    
170
<!-- R syntax highlighter -->
171
<script type="text/javascript">
172
var hljs=new function(){function m(p){return p.replace(/&/gm,"&amp;").replace(/</gm,"&lt;")}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}]}};
173
hljs.initHighlightingOnLoad();
174
</script>
175

    
176

    
177

    
178

    
179
</head>
180

    
181
<body>
182
<h1>LST Temporal Aggregation Evaluation</h1>
183

    
184
<h3>Adam M. Wilson and Giuseppe Amatulli</h3>
185

    
186
<p>A short script to explore the implications of a 3-month moving window Land Surface Temperature Climatology.</p>
187

    
188
<h1>Data</h1>
189

    
190
<p>The following data were extracted from the shapefiles used to fit the tiled interpolations that Alberto has already run. They are limited to the CONUS region and include mean monthly maximum temperature observed at the station locations, as well as the LST from the month before, the month of, and the month following.  Here&#39;s what the data look like:</p>
191

    
192
<table><thead>
193
<tr>
194
<th align="left">station</th>
195
<th align="left">month</th>
196
<th align="right">TMax</th>
197
<th align="right">lst</th>
198
<th align="right">lstp1</th>
199
<th align="right">lstm1</th>
200
<th align="left">tile</th>
201
</tr>
202
</thead><tbody>
203
<tr>
204
<td align="left">USC00100528</td>
205
<td align="left">January</td>
206
<td align="right">0.9592</td>
207
<td align="right">1.880</td>
208
<td align="right">3.406</td>
209
<td align="right">1.303</td>
210
<td align="left">50.0_-115.0</td>
211
</tr>
212
<tr>
213
<td align="left">USC00100667</td>
214
<td align="left">January</td>
215
<td align="right">2.5101</td>
216
<td align="right">1.630</td>
217
<td align="right">2.264</td>
218
<td align="right">0.340</td>
219
<td align="left">50.0_-115.0</td>
220
</tr>
221
<tr>
222
<td align="left">USC00101079</td>
223
<td align="left">January</td>
224
<td align="right">0.8731</td>
225
<td align="right">1.598</td>
226
<td align="right">3.976</td>
227
<td align="right">3.330</td>
228
<td align="left">50.0_-115.0</td>
229
</tr>
230
<tr>
231
<td align="left">USC00101363</td>
232
<td align="left">January</td>
233
<td align="right">1.5249</td>
234
<td align="right">2.682</td>
235
<td align="right">3.616</td>
236
<td align="right">1.030</td>
237
<td align="left">50.0_-115.0</td>
238
</tr>
239
<tr>
240
<td align="left">USC00101956</td>
241
<td align="left">January</td>
242
<td align="right">2.6857</td>
243
<td align="right">2.354</td>
244
<td align="right">4.175</td>
245
<td align="right">3.123</td>
246
<td align="left">50.0_-115.0</td>
247
</tr>
248
<tr>
249
<td align="left">USC00102159</td>
250
<td align="left">January</td>
251
<td align="right">2.6410</td>
252
<td align="right">4.978</td>
253
<td align="right">5.748</td>
254
<td align="right">3.740</td>
255
<td align="left">50.0_-115.0</td>
256
</tr>
257
</tbody></table>
258

    
259
<h1>TMax~LST relationship seasonal variability</h1>
260

    
261
<p>First let&#39;s explore the variability of the TMax~LST relationship by month, grouped by tile.  In this plot grey lines indicate a Tmax~LST regression within each tile (stations may be present in multiple tiles). Variability among grey lines represents spatial variability in the fitted lapse rate and the heavy black line is overall mean relationship (by month).
262
<img src="http://i.imgur.com/KYSRkMV.png" alt="plot of chunk unnamed-chunk-4"/> </p>
263

    
264
<h1>Comparison of monthly means with 3-month moving window</h1>
265

    
266
<p>Here we compare the monthly means with a three month moving window (e.g. January mean LST with December-February mean LST).  Note that the relationship is very good (R<sup>2</sup> &gt;0.95) but slightly weaker in winter months, likely due primarily to seasonal minimums.  Heavy black line is 1:1, and thin line is the fitted regression.</p>
267

    
268
<p><img src="http://i.imgur.com/kInQcBr.png" alt="plot of chunk unnamed-chunk-5"/> </p>
269

    
270
<h1>Comparison of monthly means with 2-month moving window that does not include the month</h1>
271

    
272
<p>Here we compare the monthly means with a two month moving window that does not include the month of interest (e.g. January mean LST with December and February mean LST, but not including the January LST).  Heavy black line is 1:1, and thin line is the fitted regression.</p>
273

    
274
<p><img src="http://i.imgur.com/8AfNmiq.png" alt="plot of chunk unnamed-chunk-6"/> </p>
275

    
276
<p>Now let&#39;s look at the differences between the 1-month and 2-month LST values.  These represent a measure of how wrong we would be if we only had data from the two surrounding months and not the month in question.  </p>
277

    
278
<pre><code class="r">histogram(dlst2$lst - dlst2$lst2m, col = &quot;grey&quot;, xlab = &quot;Anomolies (1 month - 2 month means)&quot;)
279
</code></pre>
280

    
281
<p><img src="http://i.imgur.com/7E9SHnc.png" alt="plot of chunk unnamed-chunk-7"/> </p>
282

    
283
<p>The 95th quantile of the absolute value is only 4.2, so the differences are quite small. From this analysis, it appears that broadening the temporal window will maintain a relatively consistent estimate of LST (R<sup>2</sup> ranged from 0.88-0.9) even when using only data from the surrounding months.</p>
284

    
285
<p>Let&#39;s see how the seasonal cycle is represented by these proxies for a few randomly selected stations.  Here the red line is the observed TMax data, the heavy black line represents the mean LST in that month, and the green line is a three-month moving window, while the purple line is the 2-month window (not including the month of interest).</p>
286

    
287
<p><img src="http://i.imgur.com/wOsUVGT.png" alt="plot of chunk unnamed-chunk-8"/> </p>
288

    
289
<h1>Processing Options</h1>
290

    
291
<h2>Solution 1 (relatively easy)</h2>
292

    
293
<ol>
294
<li>Spatial interpolation: fill in LST when another observation is available within  some small distance (3? 5? 7? kilometers)</li>
295
<li>Temporal interpolation: Estimate the monthly value at the 15th day by linear interpolation (in the temporal domain) considering the temporally closest observations in the previous and following month.</li>
296
</ol>
297

    
298
<p>This is more robust than the simple mean of 3 months.  For example, imagine using a simple mean in the following case:  </p>
299

    
300
<ul>
301
<li>month -1: no observation </li>
302
<li>month  0: no observation<br/></li>
303
<li>month +1: observation in the end of the month</li>
304
</ul>
305

    
306
<p>Would lead to a prediction = month +1.</p>
307

    
308
<p>Or, alternatively, one in which:</p>
309

    
310
<ul>
311
<li>month -1: has data for the full month</li>
312
<li>month 0: has no data</li>
313
<li>month +1: has only data near end of month</li>
314
</ul>
315

    
316
<p>In this case a mean would be &#39;pulled&#39; towards the mean value of month-1.  Using a fixed number of observations in t-1 and t+1 and a linear interpolation that accounts for the time of those observations would prevent </p>
317

    
318
<h2>New Solution 2 (also relatively easy)</h2>
319

    
320
<ol>
321
<li>Create 12 monthly LST climatologies if at least 5 or 10 observations are available for each pixel-month.</li>
322
<li>Fill in missing months using temporal interpolation:
323

    
324
<ul>
325
<li>Fit function through all available months (rather than just +/- 1 month) to capture seasonal variabilty.</li>
326
<li>Use some function that can capture the cyclical nature (sin, spline, polynomial, etc.)</li>
327
</ul></li>
328
</ol>
329

    
330
<p>Remember from above that the LST curves are fairly smooth (caveat: we&#39;ve only looked at locations in the CONUS dataset).  So let&#39;s try fitting a simple sin function to estimate missing months:</p>
331

    
332
<pre><code class="r">### Curve fitting function Set parameter starting values and limits
333
parms = list(phi = 5, A = 0.05)  #starting parameters
334
low = list(phi = -999, A = 0)  #lower bounds - constrain s3 to be positive 
335
high = list(phi = 999, A = 999)  #lower bounds - constrain s3 to be positive
336

    
337
## define a simple model estimating a sin function with variable shift and
338
## amplitude
339
fit.model &lt;- function(parms, x) {
340
    fit = sin(parms[[&quot;phi&quot;]] + (x * 2 * pi)) * abs(parms[[&quot;A&quot;]])
341
    fit[fit == Inf] = .Machine$double.xmax  #if Inf, replace with maximum value
342
    fit[fit == -Inf] = .Machine$double.xmin  #if Inf, replace with maximum value
343
    list(fit = fit)
344
}
345
## function to minimize the RMSE of predictions
346
minimize.model &lt;- function(parms, x, y) {
347
    ss = sum((y - fit.model(parms, x)$fit)^2)
348
    ss = min(ss, .Machine$double.xmax)  #if Inf, replace with maximum value
349
    ss = max(ss, .Machine$double.xmin)  #if -Inf, replace with maximum value
350
    list(ss = ss)
351
}
352

    
353
## Process station data and make predictions for each month
354
makepreds = function(dat, drop = NULL) {
355
    dat = na.omit(dat)
356
    dat = dat[dat$variable == &quot;lst&quot;, ]
357
    ## drop
358
    if (nrow(dat) == 0) {
359
        return(NULL)
360
    }
361
    ## drop some proportion of data if drop is specified
362
    if (!is.null(drop)) 
363
        dat = dat[sample(1:nrow(dat), round(nrow(dat) * (1 - drop))), ]
364
    dat = dat[order(dat$month), ]
365
    dat$time = (1:nrow(dat))/nrow(dat)
366
    xbar = mean(dat$value, na.rm = T)
367
    fit = optim(unlist(parms), minimize.model, x = dat$time, y = dat$value - 
368
        xbar, control = list(trace = 0, maxit = 10000), method = &quot;L-BFGS-B&quot;, 
369
        lower = low, upper = high)
370

    
371
    dat$pred = xbar + sin(fit$par[[&quot;phi&quot;]] + (dat$time * 2 * pi)) * fit$par[[&quot;A&quot;]]
372
    dat$A = fit$par[[&quot;A&quot;]]
373
    dat$phi = fit$par[[&quot;phi&quot;]]
374
    dat
375
}
376
## apply the function
377
d4 = ddply(dlst2ls, .(station), .fun = makepreds)
378
d4l = melt(d4, id.vars = c(&quot;station&quot;, &quot;month&quot;), measure.vars = c(&quot;value&quot;, &quot;pred&quot;))
379
</code></pre>
380

    
381
<p>Let&#39;s see what that looks  for the 10 example stations above. The red line is the observed LST value and the black line is the predicted LST from the sinusoidal function.
382
<img src="http://i.imgur.com/WYgS7qx.png" alt="plot of chunk unnamed-chunk-10"/> </p>
383

    
384
<p>And a histogram of the residuals for these stations:
385
<img src="http://i.imgur.com/qRekbih.png" alt="plot of chunk unnamed-chunk-11"/> </p>
386

    
387
<p>Not bad&hellip; Now let&#39;s drop 25% of the observations from each station and do it again:
388
<img src="http://i.imgur.com/n6ZP6yB.png" alt="plot of chunk unnamed-chunk-12"/> </p>
389

    
390
<p>And the  of residuals for these 10 stations with 25% of the observations removed:
391
<img src="http://i.imgur.com/8Z1Niwd.png" alt="plot of chunk unnamed-chunk-13"/> </p>
392

    
393
<h3>Apply this function to the full CONUS dataset described above.</h3>
394

    
395
<p>Now look at the distributions of the residuals for the full conus dataset.
396
<img src="http://i.imgur.com/iSdwpOe.png" alt="plot of chunk unnamed-chunk-14"/> </p>
397

    
398
<p>And summarize the residuals into RMSE&#39;s by station:
399
<img src="http://i.imgur.com/M29nxV1.png" alt="plot of chunk unnamed-chunk-15"/> </p>
400

    
401
<p>So the vast majority of stations will have a RMSE of &lt;5 using this simple method even if 25% of the data are missing.</p>
402

    
403
<h1>Next steps</h1>
404

    
405
<p>We should pick 3-4 problematic tiles and 1 tile with little missing data and explore how these various options differ in infilling LST.  It would also be ideal if we could complete the interpolation using these different datasets and compare the resulting estimates of TMax over these tiles. </p>
406

    
407
<h1>Remaining Questions</h1>
408

    
409
<ol>
410
<li>How does this variability compare to spatial variability (e.g. which is worse, spatial or temporal smoothing)?</li>
411
<li>How do the anomolies from different temporal windows vary spatially?<br/></li>
412
</ol>
413

    
414
</body>
415

    
416
</html>
417

    
(6-6/16)