feature request: provide exisiting output array to `ndinterp.affine`

With the recent activity around `ndinterp.affine` and many eyes on that functionality I will open an issue for the following feature request:

It would be great if one could provide an existing dask array as an output array through an optional argument.

I will describe the use case briefly.

Given a number of tiles (2D or 3D) each with a corresponding affine transformations that map these tiles into a joint coordinate system, it would be great if one could provide a pre-allocated output array that represents this joint coordinate system. In the output array only the chunks that map into the source array would be processed. Existing data in each output chunk would be overwritten, but only if the output pixel maps to a valid input pixel.

An application for this would for example be implementing something roughly equivalent to "fusion" step in Big Stitcher or other panorama stitching software.

I had discussed this use case before @GenevieveBuckley via email but at the time I didn't have time to look into details.
I believe @m-albert has worked on similar applications.

Start by reading the implementation and usage of `ndinterp.affine`, then trace how affine-transformed tiles produce the output dask array. Define the optional existing-output-array behavior, including chunk overwrites only where input pixels are valid, and verify it against a 2D or 3D tile-stitching case.

Written by the indexing model from the issue text.