Вычисление длин по дорожному графу, градусы - метры

[stopa85]

Здравствуйте.

Я имею раструвую карту, привязанную к

+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs

стандартной СК QGIS. Я так понимаю, что google maps выдает координаты в этой же системе и привязка была выполнена правильно.

На этот растр наложен выкторный слой (в этой же системе координат) - линии.

Мне нужно вычисление длинн линий в метрах. Как это сделать? Программно или можно чегез ГУИ?

Если я задаю в "свойства проэкта"-> "единицы измерения"-> "метры", то QGIS при измерениях просто меняет единицу измерения, не меняя при этом значения.

[stopa85]

Похоже единственный способ вычислить расстояние - это вычислить ортодромию:
http://gis-lab.info/qa/great-circles.html

И для вычисления длин есть QgsDistanceArea.

Я имею плагин к QGIS, который умеет вычислять длину кратчайшего пути между двумя произвольными точками по некоторой транспортной сети (см img2). Также имеется возможность визуально просматривать и редактировать направление движения по дорогам (см img1)

Транспортная сеть задается в виде векторного слоя:
тип геометрии - линия. Атрибуты - скорость и направление.

Остается разобраться с системами координат, оформить плагин и можно выкладывать на суд общественности.
Если кому интересно - прошу откликнуться)

Рисунок, где видно направление движения по дорогам. Для "соединения" используется притяжка.

img1.jpg (68.9 КБ) 19188 просмотров

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

Конечно интересно! Разработку под QGIS всячески приветствую. Обязательно выкладывайте, если нужно - можно сделать здесь на сайте svn и подключить скрипт автогенерации пакета плагина при обновлении svn, удобно для разработки.

А какую в итоге использовали дорожную сеть? Или пока просто абстрактный набор данных?

В ответ на ваш вопрос, вы не хотите просто спроектировать свои данные? Тогда не нужно будет использовать длины большого круга.

К сожалению пока не знаю как в QGIS, но в Arcview я делаю так:

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

aPrj = aview.getprojection
distancedeg = Feature1.distance(Feature2)
distancemet = Feature1.returnprojected(aPrj).distance(Feature2.returnprojected(aPrj))

возможно найдется что-то соответствующее этой логике, принцип тот же, просто трансформация происходит на лету, в QGIS есть возможность трансформировать на лету, значит расстояния в метрах можно считать и без хранения темы в метрах

[stopa85]

А какую в итоге использовали дорожную сеть? Или пока просто абстрактный набор данных?

Пока абстрактный набор данных. Отсканировал карту своего города (что не законно), привязал к WGS84. Дороги рисую сам, так интереснее, и так можно проверить свои результаты.

В ответ на ваш вопрос, вы не хотите просто спроектировать свои данные? Тогда не нужно будет использовать длины большого круга.

Может быть следует поступить так, подумаю.

P.S. Конечная цель - более или менее вменяемый решатель транспортной задачи (см. viewtopic.php?f=1&t=2938).

[stopa85]

В ответ на ваш вопрос, вы не хотите просто спроектировать свои данные? Тогда не нужно будет использовать длины большого круга.

Спроэктировал. В QGIS это выглядит так:

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

QgsFeature feature;
...
QgsDistanceArea da;
da.setSourceCrs( mQGisIface->mapCanvas()->mapRenderer()->destinationSrs().srsid() );
da.setProjectionsEnabled( true );

double length = da.measure( feature );

После таких преобразований, сопоставления карты и показаний одометра моей машины - я получил правдоподобную дистанцию в метрах.

А почему в метрах? Почему не в дюймах? Это зависит от СК?

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

в принципе единицы измерения могут быть любые, главное что не угловые, как в случае с сфероидальной СК, а линейные

[stopa85]

Вот еще вопросик.

Импортировал дороги и дома из opensreetmap, попробовал свой плагин и вот что видно (смотри рисунок)
Когда я считал как выехать из "дома" на дорогу - пользовался обычной аналитической геометрией, наивно посчитав что погрешность будет мала.

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

Вопрос в том какое преобразование нужно сделать чтобы получить "перпендикуляр" от здания к ближайшей дороге?

[bim2010]

Stopa85 Вы молодец !!! Даже не верится что это диплом.
Метод движущейся прямой.
Я это делал так:

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

SEEK XCODEKART+substr(STXYFILE+space(12),1,12)+STR(STABSX,10,3)+STR(STABSY,10,3)
IF FOUND()
*** ЗАХВАТ КЛИЕНТА ВЫПОЛНЕН
    SELECT MAPY
    YYRECNO=KOD
    SELECT MAPX
    XXRECNO=KOD
    * ОПРЕДЕЛЕНИЕ КООРДИНАТ ОБЛАСТИ АНАЛИЗА
    QXBMPFILE=STXYFILE   
    TIK7=2
    DO WHILE TIK7<13
      X=SUBSTR(QXBMPFILE,TIK7,1)
      IF X="T"
         TTYBMPY=VAL(ALLTRIM(SUBSTR(QXBMPFILE,TIK7+1,12)))
         TTXBMPX=VAL(ALLTRIM(SUBSTR(QXBMPFILE,2,TIK7)))
         EXIT
      ENDIF
      TIK7=TIK7+1  
    ENDDO
    QXOTFILE="L"+ALLTRIM(STR(TTXBMPX-200,10,0))+"T"+ALLTRIM(STR(TTYBMPY-200,10,0))
    QXDOFILE="L"+ALLTRIM(STR(TTXBMPX+200,10,0))+"T"+ALLTRIM(STR(TTYBMPY+200,10,0))

    TIKX3=(TTXBMPX-XXRECNO)/100
    poutn3=(TTYBMPY-YYRECNO)/100
  
    XBMPX3=TIKX3*XMAHTAB+STABSX*XMAHTAB/100
    YBMPY3=poutn3*DELITY+STABSY*XMAHTAB/100
    KORXKLI=YBMPY3
    KORYKLI=XBMPX3
***  АНАЛИЗ ЧАСТИ КАРТЫ С КООРДИНАТАМИ  ОТ "+qXOTFILE+"   ДО "+qXDOFILE
     * ЦИКЛ ПО КВАДРАТАМ
     XOT1=TTXBMPX-200   
     YOT1=TTYBMPY-200
     XDO2=TTXBMPX+200
     YDO2=TTYBMPY+200
     * эти координаты абсолютные
     KORMINX=0
 	 KORMINY=0
	 MINDLINA=999999    
	 YYYTEK2=YOT1
	 DO whil YYYTEK2<YDO2+100
		XXXTEK2=XOT1
		DO whil XXXTEK2<XDO2+100
 	                    SKVADRAT="L"+ALLTRIM(STR(XXXTEK2,10,0))+"T"+ALLTRIM(STR(YYYTEK2,10,0))
	                    * ЦИКЛ ПО БАЗЕ STREETMA ПО ТЕКУЩЕМУ КВАДРАТУ
	                    SELECT STREETMA
	                    *INDEX ON CODEKARTA+BMPFILE1+STR(BMPX1,10,3)+STR(BMPY1,10,3) TAG STREETMA2
	                     SET ORDER TO  2
	                     SEEK  XCODEKART+SUBSTR(SKVADRAT+space(12),1,12)
	                     IF FOUND()          
  			                       TIKX=(XXXTEK2-XXRECNO)/100
			                       poutn=(YYYTEK2-YYRECNO)/100

    	                                * ЕСЛИ В КВАДРАТЕ НЕСКОЛОКО ОТРЕЗКОВ
	                            DO WHILE .NOT.EOF()
			                        IF CODEKARTA==XCODEKART.AND.BMPFILE1=SKVADRAT
                               * ПОИК МИНИМАЛЬНОГО РАССТОЯНИЯ ОТ ТОЧКИ ДО ОТРЕЗКА
                               * МЕТОД ДВИЖУЩЕЙСЯ ПРЯМОЙ
                                                       XBMPX1=BMPX1
 	                                       YBMPY1=BMPY1
	                                       XBMPX1=TIKX*XMAHTAB+XBMPX1*XMAHTAB/100
	                                       YBMPY1=poutn*DELITY+YBMPY1*XMAHTAB/100
	                                       XBMPX2=BMPX2
	                                       YBMPY2=BMPY2
	                                       XBMPX2=TIKX*XMAHTAB+XBMPX2*XMAHTAB/100
	                                       YBMPY2=poutn*DELITY+YBMPY2*XMAHTAB/100
 
                                       = SEEKLine(YBMPY1,XBMPX1,YBMPY2,XBMPX2,KORXKLI,KORYKLI)
	
     		                         SELECT STREETMA             
		   	   ENDIF
           		IF EOF()
	             EXIT
		        ELSE
	             SKIP
	            ENDIF
        	ENDDO
	      ENDIF
	     XXXTEK2=XXXTEK2+100    
	     IF XXXTEK2>XDO2
	       EXIT
	     ENDIF
             ENDDO  
	  YYYTEK2=YYYTEK2+100
	  IF YYYTEK2>YDO2
	     EXIT
	  ENDIF
    ENDDO


IF MINDLINA#999999
  * KORMINX=XTEK  
   * KORMINY=YTEK
   * ПРЕОБРАЗОВАТЬ ИЗ АБСОЛЮТЫХ КООРДИНАТ В КООРДИНАТЫ ЭКРАНА
...
...

Т.е. двигаюсь по линиям транспортного графа в пределах квадрата экрана, нахожу наименьшее (дискретно)
расстояние и создаю в базе дополнительный отдельный слой к транспортному графу.

Абсолютные координаты - упрощенно говоря это метры, т.е. проекция не географическая.
И в последнюю очередь преобразую метры в координаты экрана.
Конечно запись дополнительного слоя увеличивает размер базы, но зато не надо тратить время на расчет
при решении траспортной задачи.
Вопрос: что Вы используете в качестве координат дома? Из рисунка не понятно.
Я использую слой подписи зданий. Не знаю правильно ли? Сомневаюсь.
Если Вам нужны данные для тестовой задачи чтобы не тратить время на создание данных по своему городу
рекомендую:
Data from the City of Albuquerque
http://www.cabq.gov/gis/download.html

Или попроще Берлин:
http://download.geofabrik.de/osm/europe ... in.shp.zip

Успеха!

[stopa85]

Вопрос: что Вы используете в качестве координат дома? Из рисунка не понятно.

У меня просто: кликнул мышкой - получил точку А. Кликнул мышкой получил точку Б. (нужной найти кратчайшее растояние между ними)

В качестве координаты дома можно, например, взять среднюю арифметическую координату всех вершин полигона или центр минимального включающего прямоугольника. А, кстати, Ваша идея - использовать слой "подписи" не лишена здравого смысла:)

Сначала я строю в оперативной памяти компьютера граф дорог.
Затем чтобы понять, какой дороге принадлежит точка А ищю ближайшее ребро по формуле:
Ax+By+c= "расстояние", где Ax+Bx+C=0 уравнение прямой образованное конкретным участком дороги и проверяю, принадлеит ли пересечение этой прямой и перпендикуляром опущеным на нее из точки А прямоугольнику образованному координатами ребра.

А беда в том, что формулы "прямых, расстояний и перпендикуляров ..." распространяются только на эвклидовы системы координат (на "план-схему"), но не на сферические и идеально бы было получить все эти данные именно спомощью API QGIS.

P.S. Я ж математик, мой диплом одни формулы + маленькая демонстрационная программа. То что я интегрирую свой решатель в QGIS - это моя инициатива, ради удовольствия да портфолио.

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

вернусь к предыдущему вопросу-предложению, так почему вы не хотите спроектировать данные? линии станут перпендикулярнее, чем были.

[stopa85]

вернусь к предыдущему вопросу-предложению, так почему вы не хотите спроектировать данные? линии станут перпендикулярнее, чем были.

Мы легких путей не ищем:) Я хочу получить более или менее универсальное решение - должен же быть способ сделать это и для сферической системы координат. Ну а если весь проэкт в "план-схеме" то тут, конечно, все вопросы отпыдают.

Можно спроэктировать данные на какую-либо геометрическую СК. Но как ее выбрать для общего случая? Пусть это делатет пользователь?

P.S. Взгляните на класс который я использовал для получения длинн дорог в метрах: http://doc.qgis.org/stable/classQgsDistanceArea.html
Там нет функции позволяющей найти расстояние от произвольной точки до прямой (отрезка). Нет функции "опустить перпенидикуляр". И вполне возможно что и быть не должно...

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

Я хочу получить более или менее универсальное решение - должен же быть способ сделать это и для сферической системы координат.

Я не математик, но мне кажется он есть. Просто нужно измерять на сфере не только расстояния, но и углы. В одной из статей на сайте, когда я занимался для своих задач сферической геометрией был пример расчета. Что думаете?

Пусть это делатет пользователь?

На первом этапе я бы сделал именно так. Сильно упрощает. Но универсальный подход - еще лучше.

Взгляните на класс который я использовал для получения длинн дорог в метрах: http://doc.qgis.org/stable/classQgsDistanceArea.html

Посмотрел, вроде нет. Но почему бы не добавить? Если вы все равно это делаете для себя и функция достаточно универсальная, то, если хотите, можем результат предложить для включения в ядро QGIS.

Кстати, я правильно понял, что как раз вычисления длин больших кругов присутствуют, причем на эллипсоиде (Винценти), а не просто на сфере.

[stopa85]

Кстати, я правильно понял, что как раз вычисления длин больших кругов присутствуют, причем на эллипсоиде (Винценти), а не просто на сфере.

Да присутствует. И что интересно вычисление расстояния в метрах - это итерационный процесс... тут без глубокого знания теории точности+скорости вычислений не добиться, да и нужно ли? ведь если использовать "план-схему" все и так в ажуре.

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

P.S. Я откладываю решение этой проблемы до лучших времен. Если кто-то подскажет мне решение, буду очень признателен. В настоящее время считаю важным зафиксировать результат, так что я пошел "оформлять код".

[bim2010]

Предлагаю посмотреть как делают другие
Здесь вариант решения схожий с моим и это не план-схема
http://openlayers.org/pipermail/users/2 ... 05233.html

А здесь решение наиболее Вам близкое Quantum Navigator
http://blog.qgis.org/?q=node/73
http://mapserver.sk/~wonder/qnavigator/
Если не лень напишите Ваши впечатления.

Если говорить о теме решения транспортной задачи то здесь замечательная библиотека:
http://www.boost.org/doc/libs/1_38_0/li ... tents.html

Все сводится к необходимости изучения pgrouting и сопутсвующих библиотек.
Кроме того нашел ряд литературы и статей по теме решения транспортной задачи.
Там есть ответы на мои вопросы по TSP и VRP.

Успеха!

[stopa85]

http://openlayers.org/pipermail/users/2 ... 05233.html

Тут, как я понял, используют postgresql для поиска кратчайшего маршрута. И все задачи берет на себя он.
Для Веб приложения это естественный плюс (ну не на php же маршруты считать), для всех остальных приложений + только в том, что не нужна куча ОЗУ, для того чтобы хранить очень большой граф дорог в ОЗУ.

В моей реализации для компьютера с ОЗУ 512МБ. можно хранить и обрабатывать матрицу смежности порядка 5.000.000 ребер (сколько вершин не важно), что в конечном счете выльется в решение VRP для 1000 клиенто-складо-машин.
Много это или мало даже не знаю... в любом случае можно прикрутить "подкачку из БД"
да и не проблема это - найти кратчайший маршрут между двух вершин... потому и сделал сам

А здесь решение наиболее Вам близкое Quantum Navigator
http://blog.qgis.org/?q=node/73
http://mapserver.sk/~wonder/qnavigator/
Если не лень напишите Ваши впечатления.

Действительно близкое, надо скачать и изучить. Только вот цели у нас с автором разные (я так понял автор собирается запихнуть это в GPS-навигатор, простите меня за мой английский)...
Кратчайший маршрут - для меня промежуточный результат.

Если говорить о теме решения транспортной задачи то здесь замечательная библиотека:
http://www.boost.org/doc/libs/1_38_0/li ... tents.html

Ага! библиотека замечательная.

Все сводится к необходимости изучения pgrouting и сопутствующих библиотек.

Фиг его знает

Рекомендую ознакомиться еще с работой: http://www.emis.de/journals/FPM/ps/k03/k031/k03114.pdf

Смысл работы таков. Решая задачи доставки товара точно в срок у нас может возникнуть такая ситуация:

Есть маршрут из А в Б стоимостью 10, продолжительностью 15 - пусть этот маршрут оптимален по стоимости.
Есть маршрут из А в Б стоимостью 15, продолжительностью 10 - пусть этот маршрут оптимален по времени.
Очевидно, что маршрут стоимостью 20 и продолжительностью 20 нас не интересует, но маршрут стоимостью 12 продолжительностью 12 интересует.

Вот там представлен и доказан алгоритм находящий множество интересующих нас маршрутов.