uploadimages=T  # upload all images to imgur.com for easy viewing on github

library(knitr)

```{rm,echo=FALSE}

## subset to 10 stations

dlst2ls=dlst2l[dlst2l$station%in%sample(unique(dlst2l$station),10),]

xyplot (value~month|station, data=dlst2ls, groups=variable,type="l",

## New Solution 2 (also relatively easy)

1. Create 12 monthly LST climatologies if at least 5 or 10 observations are available for each pixel-month.

2. Fill in missing months using temporal interpolation:

   * Fit function through all available months (rather than just +/- 1 month) to capture seasonal variabilty.

   * Use some function that can capture the cyclical nature (sin, spline, polynomial, etc.)

Remember from above that the LST curves are fairly smooth (caveat: we've only looked at locations in the CONUS dataset).  So let's try fitting a simple sin function to estimate missing months:

```{r,echo=T}

### Curve fitting function

### Set parameter starting values and limits

parms=list(phi=5,A=.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<-function(parms,x) {

      fit=sin(parms[["phi"]]+(x*2*pi))*abs(parms[["A"]])

      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<-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=="lst",]

  ## 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="L-BFGS-B",lower=low,upper=high)

  dat$pred=xbar+sin(fit$par[["phi"]]+(dat$time*2*pi))*fit$par[["A"]]

  dat$A=fit$par[["A"]]

  dat$phi=fit$par[["phi"]]

  dat

}

## apply the function

d4=ddply(dlst2ls,.(station),.fun=makepreds)

d4l=melt(d4,id.vars=c("station","month"),measure.vars=c("value","pred"))

```

Let's see what that looks  for the 10 example stations above. The blue line is the observed LST value and the pink line is the predicted LST from the sinusoidal function.

```{r,echo=FALSE}

xyplot (value~month|station, data=d4l, groups=variable,type="l",

        main="LST at station locations by month",

        ylab="LST (monthly)", xlab="LST (2-month Mean)",cex=.5,pch=16,auto.key=list(space="right"),scales=list(x=list(rot=45)))

```

And a histogram of the residuals for these stations:

```{r,echo=FALSE,fig.height=3}

d4$anomoly=d4$value-d4$pred

hist(d4$anomoly,col="grey",main="Histogram of residuals for predictions from sinusoidal model",xlab="Degrees C")

```

Not bad... Now let's drop 25% of the observations from each station and do it again:

```{r,echo=F}

## apply the function

d4a=ddply(dlst2ls,.(station),.fun=makepreds,drop=.25)

d4a$anomoly=d4a$value-d4a$pred

d4al=melt(d4a,id.vars=c("station","month"),measure.vars=c("value","pred"))

xyplot (value~month|station, data=d4al, groups=variable,type="l",

        main="LST at station locations by month",

        ylab="LST (monthly)", xlab="LST (2-month Mean)",cex=.5,pch=16,auto.key=list(space="right"),scales=list(x=list(rot=45)))

```

And the  of residuals for these 10 stations with 25% of the observations removed:

```{r,echo=FALSE,fig.height=3}

hist(d4a$anomoly,col="grey",main="Histogram of residuals for predictions from sinusoidal model",xlab="Degrees C")

```

### Apply this function to the full CONUS dataset described above.

Now look at the distributions of the residuals for the full conus dataset.

```{r,echo=FALSE,fig.height=3}

## apply the function

d5=ddply(dlst2l,.(station),.fun=makepreds,drop=.25)

d5$anomoly=d5$value-d5$pred

hist(d5$anomoly,col="grey",main="Histogram of residuals for predictions from sinusoidal model",xlab="Degrees C")

```

And summarize the residuals into RMSE's by station:

```{r, echo=FALSE,fig.height=3}

d5_rmse=ddply(d5,.(station),.fun=function(dat){

  sqrt(mean( (dat$pred - dat$value)^2, na.rm = TRUE))

})

hist(d5_rmse$V1,col="grey",main="Distribution of station RMSE's for predictions from sinusoidal model",xlab="RMSE")

```

So the vast majority of stations will have a RMSE of <5 using this simple method even if 25% of the data are missing.

# Next steps  

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. 

<img src="http://i.imgur.com/RZbujnI.png" alt="plot of chunk unnamed-chunk-4"/> </p>

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

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

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

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

<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 <- function(parms, x) {

    fit = sin(parms[["phi"]] + (x * 2 * pi)) * abs(parms[["A"]])

    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 <- 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 == "lst", ]

    ## 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 = "L-BFGS-B", 

        lower = low, upper = high)

    dat$pred = xbar + sin(fit$par[["phi"]] + (dat$time * 2 * pi)) * fit$par[["A"]]

    dat$A = fit$par[["A"]]

    dat$phi = fit$par[["phi"]]

    dat

}

## apply the function

d4 = ddply(dlst2ls, .(station), .fun = makepreds)

d4l = melt(d4, id.vars = c("station", "month"), measure.vars = c("value", "pred"))

</code></pre>

<p>Let&#39;s see what that looks  for the 10 example stations above. The blue line is the observed LST value and the pink line is the predicted LST from the sinusoidal function.

<img src="http://i.imgur.com/oNteXlf.png" alt="plot of chunk unnamed-chunk-10"/> </p>

<p>And a histogram of the residuals for these stations:

<img src="http://i.imgur.com/1z4IHQs.png" alt="plot of chunk unnamed-chunk-11"/> </p>

<p>Not bad&hellip; Now let&#39;s drop 25% of the observations from each station and do it again:

<img src="http://i.imgur.com/6zkd5CB.png" alt="plot of chunk unnamed-chunk-12"/> </p>

<p>And the  of residuals for these 10 stations with 25% of the observations removed:

<img src="http://i.imgur.com/c1Wyh9J.png" alt="plot of chunk unnamed-chunk-13"/> </p>

<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.

<img src="http://i.imgur.com/od4iR2S.png" alt="plot of chunk unnamed-chunk-14"/> </p>

<p>And summarize the residuals into RMSE&#39;s by station:

<img src="http://i.imgur.com/KIV2rBk.png" alt="plot of chunk unnamed-chunk-15"/> </p>

<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>

<h1>Next steps</h1>

<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>

![plot of chunk unnamed-chunk-4](http://i.imgur.com/RZbujnI.png) 

![plot of chunk unnamed-chunk-5](http://i.imgur.com/aAmbMaz.png) 

![plot of chunk unnamed-chunk-6](http://i.imgur.com/AeJ0AFq.png) 

![plot of chunk unnamed-chunk-7](http://i.imgur.com/qIlsdWo.png) 

![plot of chunk unnamed-chunk-8](http://i.imgur.com/Ve0mC1C.png) 

## New Solution 2 (also relatively easy)

1. Create 12 monthly LST climatologies if at least 5 or 10 observations are available for each pixel-month.

2. Fill in missing months using temporal interpolation:

   * Fit function through all available months (rather than just +/- 1 month) to capture seasonal variabilty.

   * Use some function that can capture the cyclical nature (sin, spline, polynomial, etc.)

Remember from above that the LST curves are fairly smooth (caveat: we've only looked at locations in the CONUS dataset).  So let's try fitting a simple sin function to estimate missing months:

```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 <- function(parms, x) {

    fit = sin(parms[["phi"]] + (x * 2 * pi)) * abs(parms[["A"]])

    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 <- 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 == "lst", ]

    ## 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 = "L-BFGS-B", 

        lower = low, upper = high)

    dat$pred = xbar + sin(fit$par[["phi"]] + (dat$time * 2 * pi)) * fit$par[["A"]]

    dat$A = fit$par[["A"]]

    dat$phi = fit$par[["phi"]]

    dat

}

## apply the function

d4 = ddply(dlst2ls, .(station), .fun = makepreds)

d4l = melt(d4, id.vars = c("station", "month"), measure.vars = c("value", "pred"))

```

Let's see what that looks  for the 10 example stations above. The blue line is the observed LST value and the pink line is the predicted LST from the sinusoidal function.

![plot of chunk unnamed-chunk-10](http://i.imgur.com/oNteXlf.png) 

And a histogram of the residuals for these stations:

![plot of chunk unnamed-chunk-11](http://i.imgur.com/1z4IHQs.png) 

Not bad... Now let's drop 25% of the observations from each station and do it again:

![plot of chunk unnamed-chunk-12](http://i.imgur.com/6zkd5CB.png) 

And the  of residuals for these 10 stations with 25% of the observations removed:

![plot of chunk unnamed-chunk-13](http://i.imgur.com/c1Wyh9J.png) 

### Apply this function to the full CONUS dataset described above.

Now look at the distributions of the residuals for the full conus dataset.

![plot of chunk unnamed-chunk-14](http://i.imgur.com/od4iR2S.png) 

And summarize the residuals into RMSE's by station:

![plot of chunk unnamed-chunk-15](http://i.imgur.com/KIV2rBk.png) 

So the vast majority of stations will have a RMSE of <5 using this simple method even if 25% of the data are missing.

# Next steps  

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.