Выборка данных для построение разреза

[magik]

Здравствуйте уважаемое сообщество. Вопрос мой не совсем про ГИС, но возможно вы сможете мне помочь.

Если коротко: Есть два массива с координатами центров сетки, нужен алгоритм (если он уже реализован на Python было бы идеально) который находил бы гридпоинт с координатами ближайшими к заданным пользователем. Наверняка задача тривиальна, но похоже что именно по этому вопросу меня в гугле забанили )

Если долго и нудно:

Исходные данные: Имеется трехмерный массив данных (x,y, глубина). Это гридированная информация из океанологической модели. Также имеются два двумерных массива с координатами широт и долгот центров ячеек (гридпоинтов) сетки.

Задача: Пользователь вводит координаты двух точек и в результате получает разрез где по x расстояние между точками, по y глубина и по z некий параметр (температура или соленость).

Мое решение: По координатам двух точек при помощи вот этой команды из пакета PyNGL http://www.pyngl.ucar.edu/Functions/Ngl.gc_interp.shtml я могу получить координаты промежуточных точек. Далее собственно нужно найти гридпоинты ближайшие к полученным мною промежуточным и крайним точкам. Сетка достаточно плотная, так что мне не нужно интерполировать значения в координаты промежуточных точек, достаточно просто взять значения из ближайших гридпоинтов.

Я написал велосипед который просто последовательно выбирает из массивов все гридпоинты с координатами в определенном окне и последовательно уменьшает это окно до тех пор пока не останется только один гридпоинт, но это очень долго и не очень правильно поскольку окно одновременно уменьшается на одинаковое значение для широт и долгот. Кому интересно может взглянуть (функция pointfind):
https://bitbucket.org/koldunovn/komod/src/tip/ut.py

Было бы здорово если вы сможете указать на алгоритм который может быстро и эффективно это делать.
Заранее спасибо.

[trir]

А ты знаешь что такое Триангуляция Делоне?
http://ru.wikipedia.org/wiki/%D0%A2%D1% ... 0%BD%D0%B5
Самый простой способ загрузить твои точки в Civil 3D/Land Desktop, построить по ним
поверхность и уже там строить разрезы или я чего то не понимаю?

[magik]

А ты знаешь что такое Триангуляция Делоне?
http://ru.wikipedia.org/wiki/%D0%A2%D1% ... 0%BD%D0%B5
Самый простой способ загрузить твои точки в Civil 3D/Land Desktop, построить по ним
поверхность и уже там строить разрезы или я чего то не понимаю?

Я не очень понимаю как мне может помочь триангуляция Делоне, не объясните поподробнее?
Вариант со сторонним софтом не подходит, мне нужно реализовать все на питоне.

[Максим Дубинин]

Voltron, не напоминает ли тебе задача ту, которую мы решали в экстрактилке данных, когда на входе точка, а нужно найти ячейку? Если да, давай поделимся кодом.

[Voltron]

Что-то похожее есть... Мы использовали аффиную трансформацию, формулы пересчета приведены здесь.

[Voltron]

Можно еще в этой теме посмотреть viewtopic.php?f=27&t=5470. Возможно, что-то подойдет

[magik]

Спасибо большое за хинты, я попробую, если что то получится отпишусь )

[trir]

Цифровая модель рельефа.
ЦМР загружена в бд MyCMR. Таблицы mypoints, mylinesid, my3dfaceid.

Функция, которая возвращает разрез по линии:
DELIMITER $$

DROP PROCEDURE IF EXISTS `GiveProfAllPoints` $$
CREATE DEFINER=`root`@`localhost` PROCEDURE `GiveProfAllPoints`(in E1 double, in N1 double,
in E2 double, in N2 double)
BEGIN
declare tk, tb double;
set tk = (N2 - N1)/(E2 - E1); set tb = N1 - tk*E1;
Select lineid, NewN, NewE, (StH + l*EndH)/(1+l) as NewH, r3
from (Select lineid, NewN, NewE, r1/r2 as l, r3, StH, EndH
from (Select lineid, NewN, NewE, StN, StE, StH, EndN, EndE, EndH, (Pow((Pow((StN-NewN), 2) + Pow((StE-NewE),2)),0.5)) as r1,
Pow((Pow((EndN-NewN),2) + Pow((EndE-NewE),2)),0.5) as r2, Pow((Pow((N1-NewN),2) + Pow((E1-NewE),2)),0.5) as r3
from (SELECT mylinekoordview.lineid, ((b - tb)/(tk - k)) as NewE,
(k*((b - tb)/(tk - k))+b) as NewN, StN, StE, StH, EndN, EndE, EndH
FROM mylinekoordview left join dlina on mylinekoordview.lineid=dlina.lineid) as NewTbl
where ((NewE between StE and EndE) or (NewE between EndE and StE))
and ((NewN between StN and EndN) or (NewN between EndN and StN))
and ((NewE between E1 and E2) or (NewE between E2 and E1))
and ((NewN between N1 and N2) or (NewN between N2 and N1))
order by r3) as OldTbl) as NextTbl;
END $$

DELIMITER ;

[magik]

Спасибо всем за помощь, я решил в итоге проблему средствами питона. Я просто нахожу расстояние от заданной точки до всех центров гридпоинтов и выбираю наименьшее

Код: Выделить всё

import os
import matplotlib.pyplot as plt
import numpy
import Ngl

def pointfind2(plat, plon, lat, lon, pdif=1):
	""" return indeces and values of the grid point closest to lat/lon point
	the same as pointfind but could be faster
	Usage: 
		pointfind(plat, plon, lat, lon, pdif = 0.5)
	Input:
		plat - latitude of the point
		plon - longitude of the point
		lat  - 2d array of latutudes
		lon  - 2d array of longitudes
		pdif  - we don't need it but leave it to have the same input as pointfind
	Output:
		indeces and values of the points that fulfil conditions 
		
		
	
	"""

	dist_min = 1000000.
	
	
	for i in range(lon.shape[0]):
		for j in range(lon.shape[1]):
			dist = Ngl.gc_dist(plat,plon,lat[i,j],lon[i,j])
			if dist_min > dist:
				dist_min = dist
				i_min = i
				j_min = j
				lat_min = lat[i,j]
				lon_min = lon[i,j]
	
	print(i_min,j_min,lat_min,lon_min)
	gg1 = i_min, j_min
	
	return(gg1, lat_min, lon_min)