See #8402 for more discussion. Bottom line is that numpy 2 changes the rules for casting between two inputs. Due to this and xarray's preference for promoting python scalars to 0d arrays (scalar arrays), xarray objects are being upcast to higher data types when they previously didn't.

I'm mainly opening this PR for further and more detailed discussion.

CC @dcherian

• [ ] Closes #8402
• [ ] Tests added
• [ ] User visible changes (including notable bug fixes) are documented in whats-new.rst
• [ ] New functions/methods are listed in api.rst

What happened?

I'm testing my code with numpy 2.0 and current main xarray and dask and ran into a change that I guess is expected given the way xarray does things, but want to make sure as it could be unexpected for many users.

Doing DataArray.where with an integer array less than 64-bits and an integer as the new value will upcast the array to 64-bit integers (python's int). With old versions of numpy this would preserve the dtype of the array. As far as I can tell the relevant xarray code hasn't changed so this seems to be more about numpy making things more consistent.

The main problem seems to come down to:

https://github.com/pydata/xarray/blob/d933578ebdc4105a456bada4864f8ffffd7a2ced/xarray/core/duck_array_ops.py#L218

As this converts my scalar input int to a numpy array. If it didn't do this array conversion then numpy works as expected. See the MCVE for the xarray specific example, but here's the numpy equivalent:

```python import numpy as np

a = np.zeros((2, 2), dtype=np.uint16)

what I'm intending to do with my xarray data_arr.where(cond, 2)

np.where(a != 0, a, 2).dtype

dtype('uint16')

equivalent to what xarray does:

np.where(a != 0, a, np.asarray(2)).dtype

dtype('int64')

workaround, cast my scalar to a specific numpy type

np.where(a != 0, a, np.asarray(np.uint16(2))).dtype

dtype('uint16')

```

From a numpy point of view, the second where call makes sense that 2 arrays should be upcast to the same dtype so they can be combined. But from an xarray user point of view, I'm entering a scalar so I expect it to be the same as the first where call above.

What did you expect to happen?

See above.

Minimal Complete Verifiable Example

```Python import xarray as xr import numpy as np

data_arr = xr.DataArray(np.array([1, 2], dtype=np.uint16)) print(data_arr.where(data_arr == 2, 3).dtype)

int64

```

MVCE confirmation

• [X] Minimal example — the example is as focused as reasonably possible to demonstrate the underlying issue in xarray.
• [X] Complete example — the example is self-contained, including all data and the text of any traceback.
• [X] Verifiable example — the example copy & pastes into an IPython prompt or Binder notebook, returning the result.
• [X] New issue — a search of GitHub Issues suggests this is not a duplicate.
• [X] Recent environment — the issue occurs with the latest version of xarray and its dependencies.

Relevant log output

No response

Anything else we need to know?

Numpy 1.x preserves the dtype.

```python In [1]: import numpy as np

In [2]: np.asarray(2).dtype Out[2]: dtype('int64')

In [3]: a = np.zeros((2, 2), dtype=np.uint16)

In [4]: np.where(a != 0, a, np.asarray(2)).dtype Out[4]: dtype('uint16')

In [5]: np.where(a != 0, a, np.asarray(np.uint16(2))).dtype Out[5]: dtype('uint16') ```

Environment

I'm helping @joleenf debug an issue where some old code that uses xr.open_dataset no longer works since the introduction of engines or at least as far as we can tell. The main issue is that she's using code that assumes the NetCDF4 C library was compiled with HDF4 support (ex. conda-forge builds with this functionality enabled). So in this case netCDF4.Dataset("my_file.hdf") can actually read the HDF4 file through the NetCDF4 C library.

However, with xr.open_dataset("my_file.hdf") will fail because xarray (or rather the netcdf4 engine) doesn't know that it could potentially read HDF4 files. This PR adds the .hdf extension to the 'netcdf4' engine to allow this to be automatic without needing engine='netcdf4' to be specified.

What do people think? I didn't want to put any more work into this until others weighed in.

• [ ] Closes #xxxx
• [ ] Tests added
• [ ] User visible changes (including notable bug fixes) are documented in whats-new.rst
• [ ] New functions/methods are listed in api.rst

I've turned on testing my Satpy package with unstable or pre-releases of some of our dependencies including numpy and xarray. I've found one error so far where in previous versions of xarray it was possible to assign to the numpy array taken from a DataArray.

MCVE Code Sample

```python import numpy as np import dask.array as da import xarray as xr

data = np.arange(15, 301, 15).reshape(2, 10) data_arr = xr.DataArray(data, dims=('y', 'x'), attrs={'test': 'test'}) data_arr = data_arr.copy() data_arr = data_arr.expand_dims('bands') data_arr['bands'] = ['L'] n_arr = np.asarray(data_arr.data) n_arr[n_arr == 45] = 5

```

Which results in:

```

ValueError Traceback (most recent call last) <ipython-input-12-90dae37dd808> in <module> ----> 1 n_arr = np.asarray(data_arr.data); n_arr[n_arr == 45] = 5

ValueError: assignment destination is read-only ```

Expected Output

A writable array. No error.

Problem Description

If this is expected new behavior then so be it, but wanted to check with the xarray devs before I tried to work around it.

Output of xr.show_versions()

Problem description

This issue is to start the discussion for a feature that would be helpful to a lot of people. It may not necessarily be best to put it in xarray, but let's figure that out. I'll try to describe things below to the best of my knowledge. I'm typically thinking of raster/image data when it comes to this stuff, but it could probably be used for GIS-like point data.

Geographic data can be projected (uniform grid) or unprojected (nonuniform). Unprojected data typically has longitude and latitude values specified per-pixel. I don't think I've ever seen non-uniform data in a projected space. Projected data can be specified by a CRS (PROJ.4), a number of pixels (shape), and extents/bbox in CRS units (xmin, ymin, xmax, ymax). This could also be specified in different ways like origin (X, Y) and pixel size. Seeing as xarray already computes all coords data it makes sense for extents and array shape to be used. With this information provided in an xarray object any library could check for these properties and know where to place the data on a map.

So the question is: Should these properties be standardized in xarray Dataset/DataArray objects and how?

Related libraries and developers

• pyresample (me, @mraspaud, @pnuu)
• verde and gmt-python (@leouieda)
• metpy (@dopplershift)
• geo-xarray (@andrewdhicks)
• rasterio
• cartopy

I know @WeatherGod also showed interest on gitter.

Complications and things to consider

1. Other related coordinate systems like ECEF where coordinates are specified in three dimensions (X, Y, Z). Very useful for calculations like nearest neighbor of lon/lat points or for comparisons between two projected coordinate systems.
2. Specifying what coords arrays are the CRS coordinates or geographic coordinates in general.
3. If xarray should include these properties, where is the line drawn for what functionality xarray supports? Resampling/gridding, etc?
4. How is the CRS object represented? PROJ.4 string, PROJ.4 dict, existing libraries CRS object, new CRS object, pyproj.Proj object?
5. Affine versus geotransforms instead of extents: https://github.com/mapbox/rasterio/blob/master/docs/topics/migrating-to-v1.rst#affineaffine-vs-gdal-style-geotransforms
6. Similar to 4, I never mentioned "rotation" parameters which some users may want and are specified in the affine/geotransform.
7. Dynamically generated extents/affine objects so that slicing operations don't have to be handled specially.
8. Center of pixel coordinates versus outer edge of pixel coordinates.

Code Sample, a copy-pastable example if possible

I'm working on a proof-of-concept for the geoxarray project where I'd like to store coordinate reference system (CRS) information in the coordinates of a DataArray or Dataset object. I'd like to avoid subclassing objects and instead depend completely on xarray accessors to implement any utilities I need.

I'm having trouble deciding what the best place is for this CRS information so that it benefits the user; .coords made the most sense. My hope was that adding two DataArrays together with two different crs coordinates would cause an error, but found out that since crs is not a dimension it doesn't get treated the same way; even when changing join method to 'exact'.

```python from pyproj import CRS import xarray as xr import dask.array as da

crs1 = CRS.from_string('+proj=lcc +datum=WGS84 +lon_0=-95 +lat_0=25 +lat_1=25') crs2 = CRS.from_string('+proj=lcc +datum=WGS84 +lon_0=-95 +lat_0=35 +lat_1=35')

a = xr.DataArray(da.zeros((5, 5), chunks=2), dims=('y', 'x'), coords={'y': da.arange(1, 6, chunks=3), 'x': da.arange(2, 7, chunks=3), 'crs': crs1, 'test': 1, 'test2': 2})

b = xr.DataArray(da.zeros((5, 5), chunks=2), dims=('y', 'x'), coords={'y': da.arange(1, 6, chunks=3), 'x': da.arange(2, 7, chunks=3), 'crs': crs2, 'test': 2, 'test2': 2})

a + b

Results in:

<xarray.DataArray 'zeros-e5723e7f9121b7ac546f61c19dabe786' (y: 5, x: 5)>

dask.array<shape=(5, 5), dtype=float64, chunksize=(2, 2)>

Coordinates:

* y (y) int64 1 2 3 4 5

* x (x) int64 2 3 4 5 6

test2 int64 2

```

In the above code I was hoping that because the crs coordinates are different (lat_0 and lat_1 are different and crs1 != crs2) that I could get it to raise an exception.

Any ideas for how I might be able to accomplish something like this? I'm not an expert on xarray/pandas indexes, but could this be another possible solution?

Edit: xr.merge with compat='no_conflicts' does detect this difference.