@shoyer I saw where you were going with this and here is the result I got. The trick was to preseed the dataset with an array of dtype 'object' ``` MasterShape=SimDataSet['CoorSpace'].shape; MasterShape

dealing with the non broadcasting function

SimDataSet['vector']=(['x', 'y', 'z', 't'], np.zeros(MasterShape, dtype=object)) for i in Iter.product(SimDataSet.indexes.values(), repeat=1): #print(i, vector_funcN(i)) SimDataSet['vector'].loc[dict(x=i[0], y=i[1], z=i[2], t=i[3])]=vector_funcN(*i)

full code import numpy as np import xarray as xr from sympy import * init_printing() import itertools as Iter

scale_func = lambda x, y, z, t: np.cos(1x+2y+3z-4t) scale_func(1,1,1,1)

x, y, z, t=symbols('x, y, z, t') xDirVec=Matrix([1,0,0]) vector_func=cos(1x+2y+3z-4t)*xDirVec vector_func

vector_funcN=lambdify((x, y, z, t), vector_func, dummify=False)

data point test

vector_funcN(1,1,1,1)

DomainSpaceTimeSize = 5 # using cartesian 4D SpaceTimeDensity = [10, 5] # 100 divisions in space 5 in time

x_coord = np.linspace(-DomainSpaceTimeSize, +DomainSpaceTimeSize, SpaceTimeDensity[0]) y_coord = np.linspace(-DomainSpaceTimeSize, +DomainSpaceTimeSize, SpaceTimeDensity[0]) z_coord = np.linspace(-DomainSpaceTimeSize, +DomainSpaceTimeSize, SpaceTimeDensity[0]) t_coord = np.linspace(0, +DomainSpaceTimeSize, SpaceTimeDensity[1])

dx=xr.DataArray(x_coord, dims='x') dy=xr.DataArray(y_coord, dims='y') dz=xr.DataArray(z_coord, dims='z') dt=xr.DataArray(t_coord, dims='t')

CoorEnterFunc = lambda x, y, z, t: 1+0x+0y+0z+0t

SimDataSet = xr.Dataset({'CoorSpace':(['x', 'y', 'z', 't'], CoorData)}, coords={'x':x_coord, 'y':y_coord, 'z':z_coord, 't':t_coord}) MasterShape=SimDataSet['CoorSpace'].shape; MasterShape

SimDataSet['Scaler']=(['x', 'y', 'z', 't'],scale_func(*np.meshgrid(x_coord, y_coord, z_coord, t_coord)) )

dealing with the non broadcasting function

SimDataSet['vector']=(['x', 'y', 'z', 't'], np.zeros(MasterShape, dtype=object)) for i in Iter.product(SimDataSet.indexes.values(), repeat=1): #print(i, vector_funcN(i)) SimDataSet['vector'].loc[dict(x=i[0], y=i[1], z=i[2], t=i[3])]=vector_funcN(*i)

SimDataSet.info() ```

I think there is definitely room for improvement. But I am now going to start implementing said EM example and y'alls help should get me to step 4. I step 5 should be interesting (I am not jinxing myself here) and will let y'all know when I get started on step 6 (I have done it before but very ad-hoc to say the least ) . I know holoviews works with xarray but it is high time that yt starts really taking in data from sources other than ancient obscure astronomy Fortran programs.