issue_comments: 58570858

Ok, I think I got it (for reals this time...)

``` def bcast(spat_only, coord_names): coords = [] for i, n in enumerate(coord_names): if spat_only[n].ndim != len(spat_only.dims): # Needs new axes slices = [np.newaxis] * len(spat_only.dims) slices[i] = slice(None) else: slices = [slice(None)] * len(spat_only.dims) coords.append(spat_only[n].values[slices]) return np.broadcast_arrays(*coords)

def grid_to_points2(grid, points, coord_names): if not coord_names: raise ValueError("No coordinate names provided") spat_dims = {d for n in coord_names for d in grid[n].dims} not_spatial = set(grid.dims) - spat_dims spatial_selection = {n:0 for n in not_spatial} spat_only = grid.isel(**spatial_selection)

coords = bcast(spat_only, coord_names)

kd = KDTree(zip(*[c.ravel() for c in coords]))
_, indx = kd.query(zip(*[points[n].values for n in coord_names]))
indx = np.unravel_index(indx, coords[0].shape)

return xray.concat(
        (grid.isel(**{n:j for n, j in zip(spat_only.dims, i)})
         for i in zip(*indx)),
        dim='station')

```

Needs a lot more tests and comments and such, but I think this works. Best part is that it seems to do a very decent job of keeping memory usage low, and only operates upon the coordinates that I specify. Everything else is left alone. So, I have used this on 4-D data, picking out grid points at specified lat/lon positions, and get back a 3D result (time, level, station). And I have used this on just 2D data, getting back just a 1D result (dimension='station').