Перевод долготы-широты в метры и подходы по работе в MapInfo

Эта страница опубликована в основном списке статей сайта
по адресу http://gis-lab.info/qa/dd2meters-mapinfo.html

Описывается решение задачи перевода координат заданных в градусах в координаты в метрах. Приемы статьи могут быть использованы для пересчета координат между различными проекциями, если те могут быть заданы стандартными средствами ГИС MapInfo.

Вступление

Цель данной статьи:

1. показать как выполнить конкретную задачу по вычислению расстояний между точками, заданными в геодезической (географической) системе координат в градусах, что является нормой для приемников сигнала GPS/Глонасс, через перевод координат в спроецированные координаты в метрах.
2. продемонстрировать некоторые типичные подходы при работе в ГИС MapInfo, в том числе и неявные способы, открывающиеся при использовании окна Mapbasic, в который MapInfo выводит практически все команды, вызываемые через стандартный интерфейс.

Статья написана по мотивам обсуждения на форуме.

Постановка задачи

Имеется набор данных, состоящих из координат точек в формате долгота/широта/высота + атрибуты. В рассматриваемом случае атрибуты — время создания точки и некоторые специфические параметры, связанные с работой технического средства. Данные даны в виде плоской таблицы или в виде приводимом к плоской таблице. В нашем случае в формате Excel. Требуется произвести вычисления расстояния между соседними точками и сохранить результат так, что бы его можно было использовать в Excel.

Общее описание алгоритмов

Для решения задачи можно использовать:

• Любую ГИС, позволяющую вычислять координаты точек в различных системах координат (или проекциях);
• Программные средства командной строки, позволяющие производить вычисления над плоскими файлами. Например, утилиту командной строки ogr2ogr или иную, использующую библиотеку PROJ4, или ее аналог;
• Дополнительные модули или формулы в Excel, которые позволят произвести вычисление координат, заданных в градусах, в координаты на плоскости в некоторой проекции. В принципе в тот же Excel можно вставить вызов библиотек PROJ4.

Почему MapInfo + Excel?

В этой статье будет рассмотрен набор манипуляций с данными с использованием ГИС MapInfo. С последующей обработкой результатов в Excel.

Причина выбора этого ПО:

• Предпочтение автора статьи;
• MapInfo хорошо умеет читать данные в формате Excel;
• MapInfo, за счет встроенного интерпретатора языка MapBasic, хорошо автоматизируется;
• Excel идеально подходит для между строчных вычислений, что является трудноразрешимой проблемой для всего, что построено на плоских таблицах с независимыми строками;

Основные приемы использованные в работе

1. Главный прием — 95% процентов команд графического ядра ГИС MapInfo являются командами на языке MapBasic. При стандартном выполнении через меню и «кнопки» эти команды формируются в GUI, выводятся в окно MapBasic в MapInfo и одновременно передаются на выполнение ядру ГИС MapInfo, которое можно считать интерпретатором команд MapBasic.
2. Если выбрать соответствующую команду в окне MapBasic в MapInfo, выделив несколько строк подряд, или просто поместив курсор в конец одной строки, и нажать «Ввод» на клавиатуре, то выбранные команды будут переданы интерпретатору команд MapBasic и они дадут к такой же результат, что и вызов команд через графический интерфейс.
3. Поэтому используя команды от стандартных элементов меню MapInfo, выводимые в MapBasic, как шаблоны для собственных команд, или используя любые команды и операторы языка MapBasic, за исключением операторов проверки условий (например, IF) или циклов (FOR или DO), можно получать большую гибкость при работе с данными в этой ГИС.
4. За пределами этого рассмотрения остается возможность посылки аналогичных команд через OLE интерфейс MapInfo. В этом случае они выполняются так же, как если бы их вводили в окно MapBasic.
5. Так же работа основана на том, что MapInfo хорошо читает многие стандартные форматы табличных данных, и в некоторые из них позволяет сохранять свои таблицы.
6. В свою очередь, Excel легко позволяет выписать функцию на рабочем листе, а потом применить ее к соседним клеткам с сохранением условной или безусловной адресации между ячейками. Например, для того что бы формула, заданная в 3-м столбце строки, выполняющая расчет на основе первых двух ячеек строки, была аналогично использована в следующей строке, достаточно просто растянуть ячейку в 1:3 или С1 следующую ячейку снизу — 2:3 или С2.

Порядок обработки данных

1. Открываем таблицу Excel. Приводим ее к виду набора однородных данных – заголовки должны быть у всех столбцов, имеющих данные, и располагаться в первой строке. Столбцы без данных должны быть удалены. 2. Сохраним файл с простым именем, например D20120503.xls. С вот таким заголовком:

Время	Высота над УМ	Датчик скорости	Координата-1	Координата-2	Обороты двигателя	Конечная
2012-05-03 07:31:40	521	18	51.822845	107.679138	1482	мелькомбинат
2012-05-03 07:32:27	525	22	51.81982	107.682434	1208.4

3. Загрузим данные в MapInfo (галочку – создать копию ставим обязательно):

Примечание: для обработки данных использована ГИС MapInfo версии 9.5, однако описанные алгоритмы работают в версиях от 7.8 (а должны и ранее). Исключения могут быть связаны только с форматом файла Excel, в котором созданы исходные данные, и который умеете читать "напрямую" программа Мапифно.

4. Используем заголовки для наименования столбцов:

5. Подтвердим выбранные самой MapInfo типы данных:

6. Откроем таблицу для просмотра в виде списка, т.к. пока у нас есть только атрибутика, но нет геометрии:

7. Выполним фильтрацию данных, т.к. в некоторых строках есть некорректные данные, т.е. точки с отрицательной высотой над уровнем моря, отрицательными координатами или отрицательной скоростью:

Такие данные не позволили бы нам произвести расчет в Excel даже если бы мы использовали функции для пересчета сферических координат, видимо полученных от устройства GPS, в плановые метровые. Для этого создадим и выполним SQL запрос:

8. Получим вот такой результат:

9. Удалим из таблицы все строки, содержащие некорректные данные, нажав Del. После этого сохраним таблицу и упакуем ее.

10. Теперь вычислим несколько параметров наших данных. Опять же через SQL запрос:

11. И сразу еще один запрос на данных предыдущего запроса:

12. В результат получим:

Все эти данные нам понадобятся, что бы создать геоинформацию к нашей Excel таблице.

На основании выполненных вычислений, мы получили пределы, в которые вписываются все координаты измерений. Используя явное указание границ (см. ниже) при задании системы координат (проекции) таблицы MapInfo, мы существенно повышаем точность измерения. Вычислив среднее значение "Координата_2", мы получили возможность выбрать центральный меридиан проекции "Поперечная Меркатора" (или проекция "Гаусса-Крюгера"). Это сделано для того, что бы уменьшить различие между спроецированными расстояниями (переведенными на плоскость спроецированной поверхности) и расстояниями, вычисленными с учетом кривизны Земли, т.к. в дальнейшем при переходе от координат в градусах к координатам в метрах, для вычисления расстояний между точками мы будем использовать обычную формулу Пифагора для 3-х мерного пространства.

13. Откроем окно "Управление параметрами и полями таблицы" в MapInfo добавим два поля для плановых (плоских, метровых ) координат – X,Y:

14. В этом же окне присвоим таблице систему координат - установим проекцию «Долгота/Широта WGS84»:

15. В окне MapBasic должны появиться две команды:

Alter Table "D20120503" ( add X Float,Y Float ) Interactive

Create Map For D20120503 CoordSys Earth Projection 1, 104

16. Для нас важна вторая строка, которая и создает систему координат (СК) таблицы. Нас она не устраивает, т.к. создает СК, которая охватывает весь мир, что плохо для точности наших данных. Обсуждение проблемы смотрим вот здесь Изменим строку команды, что она приобрела вот такой вид – добавим ей границы координат, полученные в п.12, с некоторым запасом:

Create Map For D20120503 CoordSys Earth Projection 1, 104 Bounds (107.5, 51.7) (107.8, 52)

После этого вставим эту строку в окно MapBasic, встанем на ее конец (или выделим ее целиком) и нажмем «Ввод». Далее для любой строки команд, необходимость выделения вставленной команды и ее ввода будет подразумеваться.

17. Теперь создадим точки на месте координат:

18. Кнопку «Проекция» - не трогаем. Тогда проекция вычисления координат и проекция таблицы – совпадут по умолчанию.

19. В результате получим вот такую карту:

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

20. Пришло время пересчитать координаты точек. Но для этого надо установить прямоугольную систему координат. Сделаем это на примере проекции UTM. Причем опять с использованием шаблона из стандартных проекций, которые потом подправим «под себя». Для этого сделаем самое простое – установим проекцию окна карты:

21. Выберем самую первую – UTM-1N:

22. В окне карты получится ужас, но выставленные параметры нам нужны только как образец, который мы исправим под свои нужды.

23. В окне MapBasic’а будет выведена следующая команда, соответствующая установке для текущего окна карт выбранной нами проекции:

Set Map XY Units "m" CoordSys Earth Projection 8, 104, "m", -177, 0, 0.9996, 500000, 0

24. Заменим эту проекцию на нужную нам, взяв новый центральный меридиан из среднего значения координат таблицы (п.13):

Set Map XY Units "m" CoordSys Earth Projection 8, 104, "m", 107.63, 51, 1, 250000, 0 Bounds (0,0) (500000, 500000)

26. Обращаю внимание, что мы в параметрах проекции заменили их практически все – оставив только первые «8, 104, "m",» (что означает проекция «поперечная Меркатора», на датуме WGS84, единица измерения координат «метры»):

a. Затем задали – центральный меридиан 107.63;

b. 0 (экватор) заменили на минимальную южную границу = 51 (градус северной широты);

c. 0.9996 – стандартный масштабный коэффициент для проекции UTM заменили на 1, т.к. наши измерения предполагаются возле центрального меридиана (ЦМ) проекции, а стандартный масштабный коэффициент для UTM уменьшает расстояния возле ЦМ.;

d. Уменьшили ложное восточное смещение – заменили 500000 на 250000 (что то же для наших данных). Принципиального значения это не имеет. Могли и 0 поставить, округления во внутренних расчетах MapInfo для таких величин будут незначительны. Вот в случаях, когда используется проекция СК-42 или СК-95, то «добавочные» миллионы, задающие зону, могут оказать неприятное воздействие на MapInfo, да и на Excel - то же.;

e. 0 – ложное смещение на север, оставили 0, т.к. уже сдвинули начало системы координат задав смещение в п. 26.b;

f. Задали границы карты - Bounds (0,0) (500000, 500000). Это существенно для Mapinfo только при вычислении координат. Для окна карты не имеет никакого значения. Границы выставили из предположения, что все наши координаты впишутся в квадрат 500 на 500 км.

27. Проверим, что наши параметры вернули карте нормальный вид, такой как она имела изначально – дадим команду показать все слои карты в нашем окне:

28. Откроем таблицу в виде списка, что бы видеть результаты вычислений.

29. Теперь установим систему координат, которую MapInfo , будет использовать для расчета координат. Для этого введем следующую строку в окно MapBasic, используя ранее полученную строку как шаблон:

Set CoordSys Earth Projection 8, 104, "m", 107.63, 51, 1, 250000, 0 Bounds (0, 0) (500000, 500000)

30. Пришло время вычислить координаты. Для этого используем команду:

Update D20120503 set col8=centroidx(obj),col9=centroidy(obj)

31. COL8 и COL9 – это псевдонимы для 8 и 9 колонки, соответственно. Вместо них обычно используют настоящие названия колонок, но так проще и короче.

32. Наша открытая списком таблица обновится и приобретет вот такой вид:

33. Обязательно сохраним ее.

34. Теперь нам надо вернуть данные в Excel. Проще всего это сделать, сохранив данные в тестовый формат, т.к. экспорта в Excel в MapInfo – нет. По крайней мере, в моей версии.

35. Выберем вот такие параметры экспорта:

36. Сознательно добавим букву «а» к имени, предложенному Mapinfo, т.к. ниже мы загрузим экспорт обратно, а если имя текстового файла совпадет с именем таблицы, MapInfo не даст нам ее открыть.

37. Откроем экспортированную таблицу:

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

Create Map For D20120503a CoordSys Earth Projection 8, 104, "m", 107.63, 51, 1, 250000, 0 Bounds (0,0) (500000, 500000)

39. То, что это удалось можно увидеть, добавив новую таблицу в окно карты:

До момента создания пространственной привязки таблица не появляется в списке слоев.

40. Создадим точки и в этой таблице, только изменим их форму и цвет:

41. И естественно выберем новые расчетные столбцы в качестве координат создаваемых точек. После нажатия на «OK» - окно карты измениться. В нем появятся новые геокодированные точки:

42. Максимально приблизим одну из точек, что бы стало заметно различие между старыми (фиолетовыми) и новыми (зелеными) точками:

Видим, что оно составило 4 мм. Это плата за пересчет и округление. На мой взгляд, в большей степени за округление, с которым Mapinfo вывел данные в текстовый файл – всего 2 символа после запятой (точки). Эту точность можно принудительно увеличивать с помощью функции Format$, но это в этом случае не играет большой роли.

43. Все. С MapInfo закончили. Надо возвращаться в Excel. Пункты с 37-го по 42-ой можно опустить. Они были необходимы только для проверки корректности выполненных нами манипуляций с координатными системами.

44. Открываем результат экспорта как текстовый файл в Excel.

45. Получим вот такую таблицу в Excel и добавим в нее две колонки «расстояние» (между соседними точками) и «путь» (сумму расстояний):

46. Заполним ячейки в новых колонках формулой:

a. ячейка J2 ‑ «0»

b. ячейка K2 – «0»

c. В ячейке J3 напишем формулу - «=КОРЕНЬ((H2-H3)*(H2-H3)+(I2-I3)*(I2-I3)+(B2-B3)*(B2-B3))» (Теорема Пифагора для 3-х мерного пространства)

d. в ячейке K3 напишем «=K2+J3» - это суммарный путь, пройденный до этой точки.

e. Теперь «растянем» формулу из ячеек J3:K3 на все оставшееся заполненное пространство листа Excel. Существуют и более правильные способы задания формулы на столбцы, но они слегка длиннее и мудренее, чем просто «взять за угол и растянуть».

47. Мы получили искомый результат :) Все строки содержат расстояние от точки до предшествующей ей, а так же пройденный путь.

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

49. Теперь, по уму, надо бы построить линию по этим точкам, что бы убедиться, что расстояние мы измерили без погрешности.

50. Для этого воспользуемся программой перевода последовательности точек в полилинию (или полигон):

Примечание: приведенная программа является собственной разработкой (с)2011. В новые версии MapInfo входит улучшенный аналог, который позволяет задавать построение множества геобъектов из наборов точек, в то время как эта программа создает единственный объект на основании входного списка. Более функциональные программы используют для определения порядка создания вершин два индекса - первый задает объекты, второй, порядок вершин в объекте. Приведенные в статье программа не использует индексов, а создает вершины в том же порядке как они записаны в исходных файлах. Кроме того возможно создание геообъектов через программы библиотеки GDAL/ORG, как это описано в статье Конвертация данных из CSV в SHP и обратно с OGR 51. Вот что получим, измерив длину новой линии:

Или с учетом кривизны Земли:

52. Тут надо помнить, что вычисления в MapInfo производятся только по двумерным координатам при геообъекте, высоту h она считать не умеет и для узлов не хранит. Результат вычислений в Exec = 226 781 m.

53. Если я ничего не перепутал то, то что казалось набором последовательных точек, оказалось набором вершин от разных проходов по одному пути, пройденных в разное время:

Ключевые слова: MapInfo, MapBasic команды, Select SQL, долгота/широта, преобразование координат, пользовательские проекции