Feature request: vector cross product

[smartass101]

xarray currently has the xarray.dot() function for calculating arbitrary dot products which is indeed very handy.
Sometimes, especially for physical applications I also need a vector cross product. I' wondering whether you would be interested in having xarray.cross as a wrapper of numpy.cross. I currently use the following implementation:

def cross(a, b, spatial_dim, output_dtype=None):
    """xarray-compatible cross product
    
    Compatible with dask, parallelization uses a.dtype as output_dtype
    """
    # TODO find spatial dim default by looking for unique 3(or 2)-valued dim?
    for d in (a, b):
        if spatial_dim not in d.dims:
            raise ValueError('dimension {} not in {}'.format(spatial_dim, d))
        if d.sizes[spatial_dim] != 3:  #TODO handle 2-valued cases
            raise ValueError('dimension {} has not length 3 in {}'.format(d))
        
    if output_dtype is None: 
        output_dtype = a.dtype  # TODO some better way to determine default?
    c = xr.apply_ufunc(np.cross, a, b,
                       input_core_dims=[[spatial_dim], [spatial_dim]], 
                       output_core_dims=[[spatial_dim]], 
                       dask='parallelized', output_dtypes=[output_dtype]
                      )
    return c

Example usage

import numpy as np
import xarray as xr
a = xr.DataArray(np.empty((10, 3)), dims=['line', 'cartesian'])
b = xr.full_like(a, 1)
c = cross(a, b, 'cartesian')

Main question

Do you want such a function (and possibly associated DataArray.cross methods) in the xarray namespace, or should it be in some other package? I didn't find a package which would be a good fit as this is close to core numpy functionality and isn't as domain specific as some geo packages. I'm not aware of some "xrphysics" package.

I could make a PR if you'd want to have it in xarray directly.

Output of xr.show_versions()

# Paste the output here xr.show_versions() here INSTALLED VERSIONS ------------------ commit: None python: 3.7.3 (default, Mar 27 2019, 22:11:17) [GCC 7.3.0] python-bits: 64 OS: Linux OS-release: 4.9.0-9-amd64 machine: x86_64 processor: byteorder: little LC_ALL: None LANG: en_US.UTF-8 LOCALE: en_US.UTF-8 libhdf5: 1.10.4 libnetcdf: 4.6.1

xarray: 0.12.3
pandas: 0.24.2
numpy: 1.16.4
scipy: 1.3.0
netCDF4: 1.4.2
pydap: None
h5netcdf: 0.7.4
h5py: 2.9.0
Nio: None
zarr: None
cftime: 1.0.3.4
nc_time_axis: None
PseudoNetCDF: None
rasterio: None
cfgrib: None
iris: None
bottleneck: 1.2.1
dask: 2.1.0
distributed: 2.1.0
matplotlib: 3.1.0
cartopy: None
seaborn: 0.9.0
numbagg: None
setuptools: 41.0.1
pip: 19.1.1
conda: 4.7.11
pytest: 5.0.1
IPython: 7.6.1
sphinx: 2.1.2

[dschwoerer]

Very useful 👍
I would add:

    try:
        c.attrs["units"] = a.attrs["units"] + '*' + b.attrs["units"]
    except KeyError:
        pass

to preserve units - but I am not sure that is in scope for xarray.

[keewis]

it is not, but we have been working on unit aware arrays with pint. Once that is done, unit propagation should work automatically.