Interpolation fails with single-precision data (at least sometimes?)

I'm trying to use reproject with some large images. Even though the FITS files are 12 GB and my system has 64 GB of RAM, I'm running out of memory, so I've been investigating ways to be more memory efficient. I'm currently trying to implement the memmap recommendation [from the docs][memmap].

[memmap]: https://reproject.readthedocs.io/en/stable/celestial.html#very-large-datasets

In so doing, reproject start rejected an innocuous-seeming transformation:

```python
out = np.memmap(filename='temp.dat', mode='w+', shape=shape_out, dtype=input.dtype)
reproject_interp(
(input, wcs),
output_array=out,
output_projection=wcs_out,
shape_out=shape_out,
return_footprint=False,
)
```

with:

> ValueError: An output array of a different type than the input array was specified, which will create an undesired duplicate copy of the input array in memory.

I'm explicitly initializing my arrays to have the same dtype — how is this happening?

The culprit turns out to be:

https://github.com/astropy/reproject/blob/43b0d8a4a5641cfbe6adbc3b1f2d7598f2fd5930/reproject/interpolation/core.py#L73

The input array is float32 data, and on my machine, at least (Linux x86_64 numpy 1.21.2), this line converts it to float64. This is counterintuitive at best and is also going to be doubling my memory consumption.

I'm pretty sure that this is a mistake — unless I'm missing something? I wanted to ask before working on a PR, although I'm pretty sure that there's a one-line way to ensure float-ness without unnecessarily converting.

Start in reproject/interpolation/core.py at the line linked in the issue, then reproduce the reproject_interp call with float32 input and output arrays. Confirm that the interpolation path accepts matching single-precision arrays without converting the input to float64 or creating the reported duplicate copy.

Written by the indexing model from the issue text.