issue_comments: 119036789

+1 for this proposal.

I made a slight modification to @WeatherGod's code so that I could use string indices for the "station" coordinate, though I'm sure there is a better way to implement this. Also note the addition of a few list calls for Python 3 compat.

``` python 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(list(zip(*[c.ravel() for c in coords])))
_, indx = kd.query(list(zip(*[points[n].values for n in coord_names])))
indx = np.unravel_index(indx, coords[0].shape)

station_da = xray.DataArray(name='station', dims='station', data=stations.index.values)

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

In [97]: stations = pd.DataFrame({'XLAT':[32.13, 32.43], 'XLONG':[-110.95, -112.02]}, index=['KTUS', 'KPHX']) stations Out[97]: XLAT XLONG KTUS 32.13 -110.95 KPHX 32.43 -112.02

In [98]: station_ds = grid_to_points2(ds, stations, ('XLAT', 'XLONG')) station_ds Out[98]: <xray.Dataset> Dimensions: (Times: 1681, station: 2) Coordinates: * Times (Times) datetime64[ns] 2015-07-02T06:00:00 ... XLAT (station) float32 32.1239 32.4337 * station (station) object 'KTUS' 'KPHX' west_east (station) int64 220 164 XLONG (station) float32 -110.947 -112.012 south_north (station) int64 116 134 Data variables: SWDNBRH (station, Times) float32 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... V10 (station, Times) float32 -2.09897 -1.94047 -1.55494 ... V80 (station, Times) float32 0.0 -1.95921 -1.87583 -1.86289 ... SWDNB (station, Times) float32 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... U10 (station, Times) float32 2.22951 1.89406 1.39955 1.04277 ... SWDDNI (station, Times) float32 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... SWDNBC (station, Times) float32 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... T2 (station, Times) float32 301.419 303.905 304.155 304.296 ... SWDDNIRH (station, Times) float32 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... U80 (station, Times) float32 0.0 1.93936 1.7901 1.63011 1.69481 ...

In [100]: station_ds.sel(station='KTUS')[['U10','V10']] Out[100]: <xray.Dataset> Dimensions: (Times: 1681) Coordinates: west_east int64 220 south_north int64 116 XLONG float32 -110.947 * Times (Times) datetime64[ns] 2015-07-02T06:00:00 ... station object 'KTUS' XLAT float32 32.1239 Data variables: U10 (Times) float32 2.22951 1.89406 1.39955 1.04277 1.16338 ... V10 (Times) float32 -2.09897 -1.94047 -1.55494 -1.34216 ... ```