issue_comments: 58550403

Starting using the above snippet for more datasets, some with interdependent coordinates and some without (so the coordinates would be 1-d). I think I have generalized it significantly...

``` def grid_to_points(grid, points, coord_names): not_spatial = set(grid.dims) - set(coord_names) spatial_selection = {n:0 for n in not_spatial} spat_only = grid.isel(*spatial_selection) coords = [] for i, n in enumerate(spat_only.dims): 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]) coords = [c.flatten() for c in np.broadcast_arrays(coords)]

kd = KDTree(zip(*coords))
_, indx = kd.query(zip(*[points[n].values for n in spat_only.dims]))
indx = np.unravel_index(indx, spat_only.shape)

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

```

I can still imagine some situations where this won't work, such as a requested set of dimensions that are a mix of dependent and independent variables. Currently, if the dimensions are independent, then the number of dimensions of each one is assumed to be 1 and np.newaxis is used for the others. Meanwhile, if the dimensions are dependent, then the number of dimensions for each one is assumed to be the same as the number of dependent variables and is merely flattened (the broadcast is essentially no-op).

I should also note that this is technically not restricted to spatial coordinates even though the code says so. Just anything that can be represented in euclidean space.