NetCDF And Python: Finding The Closest Lon/lat Index Given Actual Lon/lat Values

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.findXYindexFrom

Solution 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]

Solution 2:

You can write a fairly simple algorithm using Scipy's cdist function. You just need to compute distances from the target lat/lon coordinates (lat_value, lon_value) to the set of coordinates in the data set. Locate minimum distance and return its associated lat_index and lon_index (which might be helped with Numpy's argmin).