find indices of lat lon point on a grid using python
<p>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.</p>

<p>Basically , I am reading in an ncfile that contains 2D coordinates:</p>

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

&gt;&gt;&gt; lons.shape
(94, 83)
&gt;&gt;&gt; lats.shape
(94, 83)
</code></pre>

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

<pre><code>sel_lat=71.60556
sel_lon=-161.458611
</code></pre>

<p>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:</p>

<pre><code>latLon_pairsGrid = np.vstack(([lats.T],[lons.T])).T

&gt;&gt;&gt; latLon_pairsGrid.shape
(94, 83, 2)

distance.cdist([sel_lat,sel_lon],latLon_pairsGrid,'euclidean')
</code></pre>

<p>Any help or hints would be appreciated</p>

<p>Checkout the <a href="http://pyresample.readthedocs.io" rel="nofollow">pyresample</a> package. It provides spatial nearest neighbour search using a fast kdtree approach:</p>

<pre><code>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
</code></pre>

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

<pre><code>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
</code></pre>

python coordinates netcdf