0

I am new to python, and I can't figure out how to find the minimum distance from a given lat/lon point (which is not given from the grid, but selected by me) to a find the closest indices of a lat/lon point on a grid.

Basically , I am reading in an ncfile that contains 2D coordinates:

coords = 'coords.nc'
fh = Dataset(coords,mode='r')
lons = fh.variables['latitudes'][:,:]
lats = fh.variables['longitudes'][:,:]
fh.close()

>>> lons.shape
(94, 83)
>>> lats.shape
(94, 83)

I want to find the indices in the above grid for the nearest lat lon to the below values: 

sel_lat=71.60556
sel_lon=-161.458611

I tried to make lat/lon pairs in order to use the scipy.spatial.distance function, but I still am having problems because I did not set up the input arrays to the format it wants, but I don't understand how to do that:

latLon_pairsGrid = np.vstack(([lats.T],[lons.T])).T

>>> latLon_pairsGrid.shape
(94, 83, 2)

distance.cdist([sel_lat,sel_lon],latLon_pairsGrid,'euclidean')

Any help or hints would be appreciated

pwprnt
  • 521
  • 3
  • 9
  • 27
  • 1
    show us what `lons` and `lats` look like and what are their types. Even better: gives us an example that we can run in our computers and the expected output. – cd98 Oct 12 '16 at 23:10
  • @cd98 I'm not sure what you mean by an example, but maybe I wasn't clear enough in my question. The sel_lat, sel_lon values are given above, and those are the coordinates I need to find the nearest corresponding indices for on the grid from the 'coords.nc' file. The lons.shape = (93,83) and same with the lats.shape. And so that means the latLon_pairsGrid I made is a (93,83,2). I hope this is clear enough, I'm not sure what else to add – pwprnt Oct 12 '16 at 23:28
  • "does not have the format it wants" can you print the error please – Jules G.M. Oct 12 '16 at 23:36
  • also can you try to do `lats.T.as_matrix()` and `lons.T.as_matrix()` to your arrays inside of the np.vstack to convert them to pure numpy arrays instead of pandas stuff – Jules G.M. Oct 12 '16 at 23:38
  • @Julius I'm not using pandas, when I do lats.T.as_matrix() it gives the error ` >>> lats.T.as_matrix() Traceback (most recent call last): File "", line 1, in AttributeError: 'numpy.ndarray' object has no attribute 'as_matrix' ` For the error received in the distance.cdist command in the last line of my question it returns: ` 2027, in cdist raise ValueError('XA must be a 2-dimensional array.') ` – pwprnt Oct 12 '16 at 23:49
  • Maybe this helps: >>> latLon_pairsGrid array([[[ 55.90655136, -168.12886047], [ 55.95423126, -167.83164978], [ 56.00105667, -167.53356934], ..., [ 73.11053467, -133.51173401], [ 73.05083466, -132.90879822], [ 72.98931122, -132.31021118]]], dtype=float32) >>> – pwprnt Oct 12 '16 at 23:51
  • This may help: http://stackoverflow.com/questions/33789379/netcdf-and-python-finding-the-closest-lon-lat-index-given-actual-lon-lat-values/33792778#33792778 – N1B4 Oct 13 '16 at 14:50
  • @pwprnt I think you got an answer already, but by example I mean (hopefully) a reproducible example. Construct (maybe by taking the first few rows of lons and lats) two matrices lons and lats, paste them here and then we can easily run your code in our computers. You'll get much better answers and people that see your question later can experiment by themselves too! – cd98 Oct 17 '16 at 22:52

2 Answers2

2

Checkout the pyresample package. It provides spatial nearest neighbour search using a fast kdtree approach:

import pyresample
import numpy as np

# Define lat-lon grid
lon = np.linspace(30, 40, 100)
lat = np.linspace(10, 20, 100)
lon_grid, lat_grid = np.meshgrid(lon, lat)
grid = pyresample.geometry.GridDefinition(lats=lat_grid, lons=lon_grid)

# Generate some random data on the grid
data_grid = np.random.rand(lon_grid.shape[0], lon_grid.shape[1])

# Define some sample points
my_lons = np.array([34.5, 36.5, 38.5])
my_lats = np.array([12.0, 14.0, 16.0])
swath = pyresample.geometry.SwathDefinition(lons=my_lons, lats=my_lats)

# Determine nearest (w.r.t. great circle distance) neighbour in the grid.
_, _, index_array, distance_array = pyresample.kd_tree.get_neighbour_info(
    source_geo_def=grid, target_geo_def=swath, radius_of_influence=50000,
    neighbours=1)

# get_neighbour_info() returns indices in the flattened lat/lon grid. Compute
# the 2D grid indices:
index_array_2d = np.unravel_index(index_array, grid.shape)

print "Indices of nearest neighbours:", index_array_2d
print "Longitude of nearest neighbours:", lon_grid[index_array_2d]
print "Latitude of nearest neighbours:", lat_grid[index_array_2d]
print "Great Circle Distance:", distance_array

There is also a shorthand method for directly obtaining the data values at the nearest grid points:

data_swath = pyresample.kd_tree.resample_nearest(
    source_geo_def=grid, target_geo_def=swath, data=data_grid,
    radius_of_influence=50000)
print "Data at nearest grid points:", data_swath
sfinkens
  • 1,210
  • 12
  • 15
1

I think I found an answer, but it is a workaround that avoids calculating distance between the chosen lat/lon and the lat/lons on the grid. This doesn't seem completely accurate because I am never calculating distances, just the closest difference between lat/lon values themselves.

I used the answer to the question find (i,j) location of closest (long,lat) values in a 2D array

a = abs(lats-sel_lat)+abs(lons-sel_lon)
i,j = np.unravel_index(a.argmin(),a.shape)

Using those returned indices i,j, I can then find on the grid the coordinates that correspond most closely to my selected lat, lon value:

>>> lats[i,j]
71.490295
>>> lons[i,j]
-161.65045
Community
  • 1
  • 1
pwprnt
  • 521
  • 3
  • 9
  • 27
  • You're minimizing an L1 metric, also called [Manhattan distance](https://en.wikipedia.org/wiki/Taxicab_geometry), which is not totally unreasonable and cheaper to compute than the usual L2 metric `a=sqrt(dx**2+dy**2)`. And like you said, you're computing distance in the lat/lon domain, not geodesic or Cartesian distances, which would require coordinate transformations. If you're going to query multiple points and speed is an issue, definitely check out @Funkensieper's answer or other tree-based approaches. – Brian Hawkins May 09 '17 at 17:24