title: Applying a function along an axis of a dask array
body: <p>I'm analyzing ocean temperature data from a climate model simulation where the 4D data arrays (time, depth, latitude, longitude; denoted <code>dask_array</code> below) typically have a shape of (6000, 31, 189, 192) and a size of ~25GB (hence my desire to use dask; I've been getting memory errors trying to process these arrays using numpy).</p>
  
  <p>I need to fit a cubic polynomial along the time axis at each level / latitude / longitude point and store the resulting 4 coefficients. I've therefore set <code>chunksize=(6000, 1, 1, 1)</code> so I have a separate chunk for each grid point.</p>
  
  <p>This is my function for getting the coefficients of the cubic polynomial (the <code>time_axis</code> axis values are a global 1D numpy array defined elsewhere):</p>
  
  <pre><code>def my_polyfit(data):    
      return numpy.polyfit(data.squeeze(), time_axis, 3)
  </code></pre>
  
  <p>(So in this case, <code>numpy.polyfit</code> returns a list of length 4) </p>
  
  <p>and this is the command I thought I'd need to apply it to each chunk:</p>
  
  <pre><code>dask_array.map_blocks(my_polyfit, chunks=(4, 1, 1, 1), drop_axis=0, new_axis=0).compute()
  </code></pre>
  
  <p>Whereby the time axis is now gone (hence <code>drop_axis=0</code>) and there's a new coefficient axis in it's place (of length 4).</p>
  
  <p>When I run this command I get <code>IndexError: tuple index out of range</code>, so I'm wondering where/how I've misunderstood the use of <code>map_blocks</code>?</p>