/climate/research/LST_evaluation.html - Annotate - Environment and organisms - NCEAS Projects

root/climate/research/LST_evaluation.html @ f9695052

-e1d
+<!DOCTYPE html>
 <!-- saved from url=(0014)about:internet -->
 <html>
 <head>
 <meta http-equiv="Content-Type" content="text/html; charset=utf-8"/>
 <meta http-equiv="x-ua-compatible" content="IE=9" >
 <title>LST Temporal Aggregation Evaluation</title>
 <style type="text/css">
 body, td {
    font-family: sans-serif;
    background-color: white;
    font-size: 12px;
    margin: 8px;
+}
 tt, code, pre {
    font-family: 'DejaVu Sans Mono', 'Droid Sans Mono', 'Lucida Console', Consolas, Monaco, monospace;
+}
 h1 {
    font-size:2.2em;
+}
 h2 {
    font-size:1.8em;
+}
 h3 {
    font-size:1.4em;
+}
 h4 {
    font-size:1.0em;
+}
 h5 {
    font-size:0.9em;
+}
 h6 {
    font-size:0.8em;
+}
 a:visited {
    color: rgb(50%, 0%, 50%);
+}
 pre {
    margin-top: 0;
    max-width: 95%;
    border: 1px solid #ccc;
    white-space: pre-wrap;
+}
 pre code {
    display: block; padding: 0.5em;
+}
 code.r, code.cpp {
    background-color: #F8F8F8;
+}
 table, td, th {
   border: none;
+}
 blockquote {
    color:#666666;
    margin:0;
    padding-left: 1em;
    border-left: 0.5em #EEE solid;
+}
 hr {
    height: 0px;
    border-bottom: none;
    border-top-width: thin;
    border-top-style: dotted;
    border-top-color: #999999;
+}
 @media print {
    * {
       background: transparent !important;
       color: black !important;
       filter:none !important;
       -ms-filter: none !important;
+   }
    body {
       font-size:12pt;
       max-width:100%;
+   }
    a, a:visited {
       text-decoration: underline;
+   }
    hr {
       visibility: hidden;
       page-break-before: always;
+   }
    pre, blockquote {
       padding-right: 1em;
       page-break-inside: avoid;
+   }
    tr, img {
       page-break-inside: avoid;
+   }
    img {
       max-width: 100% !important;
+   }
    @page :left {
       margin: 15mm 20mm 15mm 10mm;
+   }
    @page :right {
       margin: 15mm 10mm 15mm 20mm;
+   }
    p, h2, h3 {
       orphans: 3; widows: 3;
+   }
    h2, h3 {
       page-break-after: avoid;
+   }
+}
 </style>
 <!-- Styles for R syntax highlighter -->
 <style type="text/css">
    pre .operator,
    pre .paren {
      color: rgb(104, 118, 135)
+   }
    pre .literal {
      color: rgb(88, 72, 246)
+   }
    pre .number {
      color: rgb(0, 0, 205);
+   }
    pre .comment {
      color: rgb(76, 136, 107);
+   }
    pre .keyword {
      color: rgb(0, 0, 255);
+   }
    pre .identifier {
      color: rgb(0, 0, 0);
+   }
    pre .string {
      color: rgb(3, 106, 7);
+   }
 </style>
 <!-- R syntax highlighter -->
 <script type="text/javascript">
 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}]}};
 hljs.initHighlightingOnLoad();
 </script>
 </head>
 <body>
 <h1>LST Temporal Aggregation Evaluation</h1>
-a3363b
+<h3>Adam M. Wilson and Giuseppe Amatulli</h3>
-e1d
+<p>A short script to explore the implications of a 3-month moving window Land Surface Temperature Climatology.</p>
 <h1>Data</h1>
 <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>
 <table><thead>
 <tr>
 <th align="left">station</th>
 <th align="left">month</th>
 <th align="right">TMax</th>
 <th align="right">lst</th>
 <th align="right">lstp1</th>
 <th align="right">lstm1</th>
 <th align="left">tile</th>
 </tr>
 </thead><tbody>
 <tr>
 <td align="left">USC00100528</td>
 <td align="left">January</td>
 <td align="right">0.9592</td>
 <td align="right">1.880</td>
 <td align="right">3.406</td>
 <td align="right">1.303</td>
 <td align="left">50.0_-115.0</td>
 </tr>
 <tr>
 <td align="left">USC00100667</td>
 <td align="left">January</td>
 <td align="right">2.5101</td>
 <td align="right">1.630</td>
 <td align="right">2.264</td>
 <td align="right">0.340</td>
 <td align="left">50.0_-115.0</td>
 </tr>
 <tr>
 <td align="left">USC00101079</td>
 <td align="left">January</td>
 <td align="right">0.8731</td>
 <td align="right">1.598</td>
 <td align="right">3.976</td>
 <td align="right">3.330</td>
 <td align="left">50.0_-115.0</td>
 </tr>
 <tr>
 <td align="left">USC00101363</td>
 <td align="left">January</td>
 <td align="right">1.5249</td>
 <td align="right">2.682</td>
 <td align="right">3.616</td>
 <td align="right">1.030</td>
 <td align="left">50.0_-115.0</td>
 </tr>
 <tr>
 <td align="left">USC00101956</td>
 <td align="left">January</td>
 <td align="right">2.6857</td>
 <td align="right">2.354</td>
 <td align="right">4.175</td>
 <td align="right">3.123</td>
 <td align="left">50.0_-115.0</td>
 </tr>
 <tr>
 <td align="left">USC00102159</td>
 <td align="left">January</td>
 <td align="right">2.6410</td>
 <td align="right">4.978</td>
 <td align="right">5.748</td>
 <td align="right">3.740</td>
 <td align="left">50.0_-115.0</td>
 </tr>
 </tbody></table>
 <h1>TMax~LST relationship seasonal variability</h1>
 <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).
-e95ca
+<img src="http://i.imgur.com/KYSRkMV.png" alt="plot of chunk unnamed-chunk-4"/> </p>
-e1d
+            Adam M. Wilson
 <h1>Comparison of monthly means with 3-month moving window</h1>
 <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>
-e95ca
+<p><img src="http://i.imgur.com/kInQcBr.png" alt="plot of chunk unnamed-chunk-5"/> </p>
-e1d
+            Adam M. Wilson
 <h1>Comparison of monthly means with 2-month moving window that does not include the month</h1>
 <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>
-e95ca
+<p><img src="http://i.imgur.com/8AfNmiq.png" alt="plot of chunk unnamed-chunk-6"/> </p>
-e1d
+            Adam M. Wilson
 <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>
 <pre><code class="r">histogram(dlst2$lst - dlst2$lst2m, col = &quot;grey&quot;, xlab = &quot;Anomolies (1 month - 2 month means)&quot;)
 </code></pre>
-e95ca
+<p><img src="http://i.imgur.com/7E9SHnc.png" alt="plot of chunk unnamed-chunk-7"/> </p>
-e1d
+            Adam M. Wilson
-            5600791d
+<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>
-e1d
+            Adam M. Wilson
 <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>
-e95ca
+<p><img src="http://i.imgur.com/wOsUVGT.png" alt="plot of chunk unnamed-chunk-8"/> </p>
-a3363b
+            Adam M. Wilson
 <h1>Processing Options</h1>
 <h2>Solution 1 (relatively easy)</h2>
 <ol>
 <li>Spatial interpolation: fill in LST when another observation is available within  some small distance (3? 5? 7? kilometers)</li>
 <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>
 </ol>
 <p>This is more robust than the simple mean of 3 months.  For example, imagine using a simple mean in the following case:  </p>
 <ul>
 <li>month -1: no observation </li>
 <li>month  0: no observation<br/></li>
 <li>month +1: observation in the end of the month</li>
 </ul>
 <p>Would lead to a prediction = month +1.</p>
 <p>Or, alternatively, one in which:</p>
 <ul>
 <li>month -1: has data for the full month</li>
 <li>month 0: has no data</li>
 <li>month +1: has only data near end of month</li>
 </ul>
 <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>
-d01d
+<h2>New Solution 2 (also relatively easy)</h2>
 <ol>
 <li>Create 12 monthly LST climatologies if at least 5 or 10 observations are available for each pixel-month.</li>
 <li>Fill in missing months using temporal interpolation:
 <ul>
 <li>Fit function through all available months (rather than just +/- 1 month) to capture seasonal variabilty.</li>
 <li>Use some function that can capture the cyclical nature (sin, spline, polynomial, etc.)</li>
 </ul></li>
 </ol>
 <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>
 <pre><code class="r">### Curve fitting function Set parameter starting values and limits
 parms = list(phi = 5, A = 0.05)  #starting parameters
 low = list(phi = -999, A = 0)  #lower bounds - constrain s3 to be positive
 high = list(phi = 999, A = 999)  #lower bounds - constrain s3 to be positive
 ## define a simple model estimating a sin function with variable shift and
 ## amplitude
 fit.model &lt;- function(parms, x) {
     fit = sin(parms[[&quot;phi&quot;]] + (x * 2 * pi)) * abs(parms[[&quot;A&quot;]])
     fit[fit == Inf] = .Machine$double.xmax  #if Inf, replace with maximum value
     fit[fit == -Inf] = .Machine$double.xmin  #if Inf, replace with maximum value
     list(fit = fit)
+}
 ## function to minimize the RMSE of predictions
 minimize.model &lt;- function(parms, x, y) {
     ss = sum((y - fit.model(parms, x)$fit)^2)
     ss = min(ss, .Machine$double.xmax)  #if Inf, replace with maximum value
     ss = max(ss, .Machine$double.xmin)  #if -Inf, replace with maximum value
     list(ss = ss)
+}
 ## Process station data and make predictions for each month
 makepreds = function(dat, drop = NULL) {
     dat = na.omit(dat)
     dat = dat[dat$variable == &quot;lst&quot;, ]
     ## drop
     if (nrow(dat) == 0) {
         return(NULL)
+    }
     ## drop some proportion of data if drop is specified
     if (!is.null(drop))
         dat = dat[sample(1:nrow(dat), round(nrow(dat) * (1 - drop))), ]
     dat = dat[order(dat$month), ]
     dat$time = (1:nrow(dat))/nrow(dat)
     xbar = mean(dat$value, na.rm = T)
     fit = optim(unlist(parms), minimize.model, x = dat$time, y = dat$value -
         xbar, control = list(trace = 0, maxit = 10000), method = &quot;L-BFGS-B&quot;,
         lower = low, upper = high)
     dat$pred = xbar + sin(fit$par[[&quot;phi&quot;]] + (dat$time * 2 * pi)) * fit$par[[&quot;A&quot;]]
     dat$A = fit$par[[&quot;A&quot;]]
     dat$phi = fit$par[[&quot;phi&quot;]]
     dat
+}
 ## apply the function
 d4 = ddply(dlst2ls, .(station), .fun = makepreds)
 d4l = melt(d4, id.vars = c(&quot;station&quot;, &quot;month&quot;), measure.vars = c(&quot;value&quot;, &quot;pred&quot;))
 </code></pre>
-e95ca
+<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.
 <img src="http://i.imgur.com/WYgS7qx.png" alt="plot of chunk unnamed-chunk-10"/> </p>
-d01d
+            Adam M. Wilson
 <p>And a histogram of the residuals for these stations:
-e95ca
+<img src="http://i.imgur.com/qRekbih.png" alt="plot of chunk unnamed-chunk-11"/> </p>
-d01d
+            Adam M. Wilson
 <p>Not bad&hellip; Now let&#39;s drop 25% of the observations from each station and do it again:
-e95ca
+<img src="http://i.imgur.com/n6ZP6yB.png" alt="plot of chunk unnamed-chunk-12"/> </p>
-d01d
+            Adam M. Wilson
 <p>And the  of residuals for these 10 stations with 25% of the observations removed:
-e95ca
+<img src="http://i.imgur.com/8Z1Niwd.png" alt="plot of chunk unnamed-chunk-13"/> </p>
-d01d
+            Adam M. Wilson
 <h3>Apply this function to the full CONUS dataset described above.</h3>
 <p>Now look at the distributions of the residuals for the full conus dataset.
-e95ca
+<img src="http://i.imgur.com/iSdwpOe.png" alt="plot of chunk unnamed-chunk-14"/> </p>
-d01d
+            Adam M. Wilson
 <p>And summarize the residuals into RMSE&#39;s by station:
-e95ca
+<img src="http://i.imgur.com/M29nxV1.png" alt="plot of chunk unnamed-chunk-15"/> </p>
-a3363b
+            Adam M. Wilson
-d01d
+<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>
-a3363b
+            Adam M. Wilson
-d01d
+<h1>Next steps</h1>
-a3363b
+            Adam M. Wilson
-d01d
+<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>
-e1d
+            Adam M. Wilson
 <h1>Remaining Questions</h1>
 <ol>
 <li>How does this variability compare to spatial variability (e.g. which is worse, spatial or temporal smoothing)?</li>
-a3363b
+<li>How do the anomolies from different temporal windows vary spatially?<br/></li>
-e1d
+</ol>
 </body>
 </html>

Project

General

Profile