python equivalent of NCL "getind_latlon2d " · python-questions

Stream: python-questions

Topic: python equivalent of NCL "getind_latlon2d "

Stephen Yeager (Apr 02 2020 at 23:18):

Can anyone point me in the right python direction for coding the equivalent of the following NCL line:
ij = getind_latlon2d(lat2d,lon2d,lat1d,lon1d)
which returns an array of i,j indices of the lat2d and lon2d arrays (think TLAT, TLONG from POP) closest to the locations given by the lat1d and lon1d arrays ? How do you accomplish this using xarray or numpy methods?

Brian Bonnlander (Apr 02 2020 at 23:48):

https://stackoverflow.com/questions/10818546/finding-index-of-nearest-point-in-numpy-arrays-of-x-and-y-coordinates

Brian Bonnlander (Apr 02 2020 at 23:52):

If you come up with something you think others would benefit from, I can help package it.

Brian Bonnlander (Apr 02 2020 at 23:57):

This one might be better, since it references xarray: https://stackoverflow.com/questions/58758480/xarray-select-nearest-lat-lon-with-multi-dimension-coordinates

Brian Bonnlander (Apr 03 2020 at 00:04):

Maybe not better, but another data point. I like the first discussion better.

Fred Castruccio (Apr 06 2020 at 17:15):

@Stephen Yeager I usually use a Kdtree for that.

Something along those lines:

import numpy as np
from scipy import spatial

lon = np.linspace(0,360,361)
lat = np.linspace(-90,90,181)
lon2d, lat2d = np.meshgrid(lon, lat)
ii, jj = np.meshgrid(range(len(lon)), range(len(lat)))

grd = list( zip(np.ravel(lon2d), np.ravel(lat2d)) )
idx = list( zip(np.ravel(ii), np.ravel(jj)) )

pts = ([60.4, 45.1], [275.3, -34.8])

distance, index = spatial.KDTree(grd).query(pts)

pts_coord = np.asarray(grd)[index]
pts_idx = np.asarray(idx)[index]

Riley Brady (Apr 07 2020 at 20:15):

@Stephen Yeager , here's the simple function I use. If you pass in xarray objects it'll return xarray objects.

def find_indices(xgrid, ygrid, xpoint, ypoint):
    """Returns the i, j index for a latitude/longitude point on a grid.

    .. note::
        Longitude and latitude points (``xpoint``/``ypoint``) should be in the same
        range as the grid itself (e.g., if the longitude grid is 0-360, should be
        200 instead of -160).

    Args:
        xgrid (array_like): Longitude meshgrid (shape `M`, `N`)
        ygrid (array_like): Latitude meshgrid (shape `M`, `N`)
        xpoint (int or double): Longitude of point searching for on grid.
        ypoint (int or double): Latitude of point searching for on grid.

    Returns:
        i, j (int):
            Keys for the inputted grid that lead to the lat/lon point the user is
            seeking.

    Examples:
        >>> import numpy as np
        >>> x = np.linspace(0, 360, 37)
        >>> y = np.linspace(-90, 90, 19)
        >>> xx, yy = np.meshgrid(x, y)
        >>> xp = 20
        >>> yp = -20
        >>> i, j = find_indices(xx, yy, xp, yp)
        >>> print(xx[i, j])
        20.0
        >>> print(yy[i, j])
        -20.0
    """
    dx = xgrid - xpoint
    dy = ygrid - ypoint
    reduced_grid = abs(dx) + abs(dy)
    min_ix = np.nanargmin(reduced_grid)
    i, j = np.unravel_index(min_ix, reduced_grid.shape)
    return i, j

Last updated: Jan 30 2022 at 12:01 UTC