I'm investigating an error in my application that appeared after upgrading from 1.3.10 to 1.4.1. A call to reproject with a single band, 2D masked array as source now returns a 3D masked array when it previously returned a 2D array.
I suspect this change in behavior was introduced with this change, but it is unclear to me whether this is an intentional result or an unintended regression.

I can demonstrate this change using one of the tests added in the PR (test_reproject_masked_masked_output) by adding a simple assertion check on the number of dimensions of inp and out:

def test_reproject_masked_masked_output(test3d, count_nonzero, path_rgb_byte_tif):
    with rasterio.open(path_rgb_byte_tif) as src:
        if test3d:
            inp = src.read(masked=True)
        else:
            inp = src.read(1, masked=True)
    out = np.ma.masked_array(np.empty(inp.shape, dtype=np.uint8))
    out, _ = reproject(
        inp,
        out,
        src_transform=src.transform,
        src_crs=src.crs,
        dst_transform=DST_TRANSFORM,
        dst_crs="EPSG:3857",
    )
    # new assert
    assert out.ndim == inp.ndim
    assert np.count_nonzero(out[out != np.ma.masked]) == count_nonzero

When test3d = False, the test fails with:

>       assert out.ndim == inp.ndim
E       assert 3 == 2

For comparison, this assertion check passes when running this test on 1.3.10.

Is this new behavior intended?