NetCDF and Python: Finding the closest lon/lat index given actual lon/lat values

Question

I'd like to be able to find the lon/lat coordinate indices of the closest location to a lon/lat tuple. This is already available in the Java API as GridCoordSystem.findXYindexFromLatLon(), but I haven't found anything comparable in the Python API. What I'm hoping to find in the API, or write myself and contribute to the API (if useful), is something like this:

def get_indices(netcdf_dataset, lon_value, lat_value):
    '''
    :param netcdf_dataset an open NetCDF data set object
    :param lon_value a longitude value, in degrees (-180...180)
    :param lat_value a latitude value, in degrees (-90...90)
    :return indices into the lon and lat coordinate variables corresponding to the closest point to the lon/lat value arguments
    '''

    # some (trigonometry?) code here...

    return lon_index, lat_index

Maybe this isn't as complicated as I assume it is, and I can get away with just using the closest neighbor?

Thanks in advance for any comments or suggestions.

Answer 1

This is what I use for a regular lat/lon grid of decimal degrees:

def geo_idx(dd, dd_array):
   """
     search for nearest decimal degree in an array of decimal degrees and return the index.
     np.argmin returns the indices of minium value along an axis.
     so subtract dd from all values in dd_array, take absolute value and find index of minium.
    """
   geo_idx = (np.abs(dd_array - dd)).argmin()
   return geo_idx

Called like:

  in_lat = 44.67
  in_lon = -79.25
  nci = netCDF4.Dataset(infile)
  lats = nci.variables['lat'][:]
  lons = nci.variables['lon'][:]

  lat_idx = geo_idx(in_lat, lats)
  lon_idx = geo_idx(in_lon, lons)

To test:

  print lats[lat_idx]
  print lons[lon_idx]