2016-10-12 1 views
0

Я новичок в python, и не могу понять, как найти минимальное расстояние от заданной точки lat/lon (которая не указана из сетки, но выбранный мной) для поиска ближайших индексов точки lat/lon на сетке.найти индексы lat lon point на сетке с использованием python

В основном, я читаю в ncfile, который содержит 2D координаты:

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) 

Я хочу, чтобы найти индексы в приведенной выше сетки для ближайшего Lat Lon к ниже значений:

sel_lat=71.60556 
sel_lon=-161.458611 

Я попытался сделать пары lat/lon, чтобы использовать функцию scipy.spatial.distance, но у меня все еще возникают проблемы, потому что я не настраивал входные массивы в желаемый формат, но я не понимаю, как сделать это:

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

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

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

Любая помощь или советы будут оценены

+1

показать нам, что 'lons' и 'lats' выглядят и каковы их типы. Еще лучше: дает нам пример, который мы можем запустить на наших компьютерах и ожидаемый результат. – cd98

+0

@ cd98 Я не уверен, что вы подразумеваете на примере, но, возможно, я не был достаточно ясен в своем вопросе. Значения sel_lat, sel_lon приведены выше, и это координаты, которые мне нужно найти для ближайших соответствующих индексов для сетки из файла coords.nc. Lons.shape = (93,83) и тот же с lats.shape. И это значит, что latLon_pairsGrid я сделал (93,83,2). Я надеюсь, что это достаточно ясно, я не уверен, что еще добавить – pwprnt

+0

«не имеет формата, который он хочет», можете ли вы распечатать ошибку, пожалуйста, – Julius

ответ

1

Я думаю, что я нашел ответ, но это решение, которое позволяет избежать вычисления расстояния между выбранной широты/долготы и широты/Лон на сетке. Это не кажется абсолютно точным, потому что я никогда не вычисляю расстояния, а самое близкое различие между значениями lat/lon.

Я использовал ответ на вопрос 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) 

Используя эти возвращаемые индексы I, J, я могу потом найти на сетке координат, которые соответствуют наиболее близко к моему выбранной широте, значению Lon:

>>> lats[i,j] 
71.490295 
>>> lons[i,j] 
-161.65045 
+0

Вы минимизируете метрику L1, также называемую [Манхэттенское расстояние] (https://en.wikipedia.org/wiki/Taxicab_geometry), что не совсем необоснованно и дешевле вычислять, чем обычная метрика L2 'a = sqrt (дх ** 2 + ау ** 2) '. И, как вы сказали, вы вычисляете расстояние в области lat/lon, а не на геодезическом или декартовом расстояниях, что потребует преобразования координат. Если вы собираетесь запрашивать несколько точек, а скорость - проблема, обязательно проверьте ответ @ Funkensieper или другие древовидные подходы. –

1

Оформить заказ pyresample упаковка. Это обеспечивает пространственный ближайший сосед поиск, используя быстрый kdtree подход:

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 

Там также является обобщающим методом прямого получения значений данных в ближайших точках сетки:

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