Actually there's another conversion when you reuse an xarray dimension coordinate in array-like computations:

```python ds = xr.Dataset(coords={"x": np.array([1.2, 1.3, 1.4], dtype=np.float16)})

coordinate data is a wrapper around a pandas.Index object

(it keeps track of the original array dtype)

ds.variables["x"]._data

PandasIndexingAdapter(array=Float64Index([1.2001953125, 1.2998046875, 1.400390625], dtype='float64', name='x'), dtype=dtype('float16'))

This coerces the pandas.Index back as a numpy array

np.asarray(ds.x)

array([1.2, 1.3, 1.4], dtype=float16)

which is equivalent to

ds.variables["x"]._data.array()

array([1.2, 1.3, 1.4], dtype=float16)

```

The round-trip conversion preserves the original dtype so different execution times may be expected.

I can't tell much why the results are different (how much are they different?), but I wouldn't be surprised if it's caused by rounding errors accumulated through the computation of a complex formula like haversine.

The differences are larger than I would expect (order of 0.1 in some variables) but could be related to the fact that, when using different precisions, the closest grid points to the target point could change. This would eventually lead to a different value of the variable extracted from the original dataset.

Unfortunately I didn't have time to verify if it was the case, but I think this is the only valid explanation because the variables of the dataset are untouched.

It is still puzzling because, as the target points have a precision of e.g. (45.820497820, 13.003510004), I would expect the cast of the dataset coordinates from e.g. (45.8, 13.0) to preserve the 0 (45.800000000, 13.00000000), so that the closest point should not change.

Anyway, I think we're getting off-topic, thanks for the help :)