<p><code>forecast::na.interp</code> is a good approach. From the <a href="https://cran.r-project.org/web/packages/forecast/forecast.pdf" rel="nofollow">documentation</a></p>

<blockquote>
  <p>Uses linear interpolation for non-seasonal series and a periodic stl decomposition with seasonal series to replace missing values.</p>
</blockquote>

<pre><code>library(forecast)
fit &lt;- na.interp(myzoo)
fit[10]  # 32.5, vs. 31.0 actual and 32.0 from Rob Hyndman's answer
</code></pre>

<p><a href="http://arxiv.org/pdf/1510.03924.pdf" rel="nofollow">This paper</a> evaluates several interpolation methods against real time series, and finds that <code>na.interp</code> is both accurate and efficient:</p>

<blockquote>
  <p>From the R implementations tested in this paper, na.interp from the forecast package and na.StructTS from the zoo package showed the best overall results.</p>
  
  <p>The na.interp function is also not that much slower than
  na.approx [the fastest method], so the loess decomposition seems not to be very demanding in terms of computing time.</p>
</blockquote>

<p>Also worth noting that Rob Hyndman wrote the <code>forecast</code> package, and included <code>na.interp</code> after providing his answer to this question. It's likely that <code>na.interp</code> is an improvement upon this approach, even though it performed worse in this instance (probably due to specifying the period in <code>StructTS</code>, where <code>na.interp</code> figures it out).</p>
Interpolate missing values in a time series with a seasonal cycle
<p>I have a time series for which I want to intelligently interpolate the missing values. The value at a particular time is influenced by a multi-day trend, as well as its position in the daily cycle. </p>

<p>Here is an example in which the tenth observation is missing from <code>myzoo</code></p>

<pre><code>start &lt;- as.POSIXct("2010-01-01") 
freq &lt;- as.difftime(6, units = "hours") 
dayvals &lt;- (1:4)*10 
timevals &lt;- c(3, 1, 2, 4) 
index &lt;- seq(from = start, by = freq, length.out = 16)
obs &lt;- (rep(dayvals, each = 4) + rep(timevals, times = 4))
myzoo &lt;- zoo(obs, index)
myzoo[10] &lt;- NA
</code></pre>

<p>If I had to implement this, I'd use some kind of weighted mean of close times on nearby days, or add a value for the day to a function line fitted to the larger trend, but I hope there already exist some package or functions that apply to this situation?</p>

<p>EDIT: Modified the code slightly to clarify my problem. There are <code>na.*</code> methods that interpolate from nearest neighbors, but in this case they do not recognize that the missing value is at the time that is the lowest value of the day. Maybe the solution is to reshape the data to wide format and then interpolate, but I wouldn't like to completely disregard the contiguous values from the same day. It is worth noting that <code>diff(myzoo, lag = 4)</code> returns a vector of 10's. The solution may lie with some combination of <code>reshape</code>, <code>na.spline</code>, and <code>diff.inv</code>, but I just can't figure it out.</p>

<p>Here are three approaches that don't work:
<img src="https://i.stack.imgur.com/xm05G.jpg" alt="enter image description here"></p>

<p>EDIT2. Image produced using the following code.</p>

<pre><code>myzoo &lt;- zoo(obs, index)
myzoo[10] &lt;- NA # knock out the missing point
plot(myzoo, type="o", pch=16) # plot solid line
points(na.approx(myzoo)[10], col = "red")
points(na.locf(myzoo)[10], col = "blue")
points(na.spline(myzoo)[10], col = "green")
myzoo[10] &lt;- 31 # replace the missing point
lines(myzoo, type = "o", lty=3, pch=16) # dashed line over the gap
legend(x = "topleft", 
       legend = c("na.spline", "na.locf", "na.approx"), 
       col=c("green","blue","red"), pch = 1)
</code></pre>