WWW.DISSERS.RU

БЕСПЛАТНАЯ ЭЛЕКТРОННАЯ БИБЛИОТЕКА

   Добро пожаловать!

МИНИСТЕРСТВО ОБРАЗОВАНИЯ РОССИЙСКОЙ ФЕДЕРАЦИИ КРАСНОЯРСКИЙ ГОСУДАРСТВЕННЫЙ ТЕХНИЧЕСКИЙ УНИВЕРСИТЕТ АЛГОРИТМЫ И СТРУКТУРЫ ДАННЫХ ГЕОИНФОРМАЦИОННЫХ СИСТЕМ Методические указания Красноярск 2003 3

УДК 528.9(07) М21 М21 Алгоритмы и структуры данных геоинформационных систем: Методические указания для студентов специальности 071903 – «Геоинформационные системы» / Сост. И.В. Варфоломеев, И.Г. Ермакова, А.С. Савельев. Красноярск: КГТУ, 2003, 34 с.

Печатается по решению редакционно-издательского совета университета © КГТУ, 2003 Печатается в авторской редакции 4 Общие сведения Методические указания подготовлены в соответствии с рабочей про граммой по курсу «Проектирование геоинформационных систем» для сту дентов Красноярского Государственного технического университета специ альности 071903 – «Геоинформационные системы».

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

В этих методических указаниях авторы преследуют две цели. Во первых, студенты должны не только уметь работать с современным про граммным обеспечением ГИС, но и понимать, как внутри системы выполня ется та или иная операция. Во-вторых, проектирование ГИС вовсе не ограни чено использованием существующего коммерческого программного обеспе чения.

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

1. Структуры пространственных данных ГИС 1.1. Хранение растровых данных В геоинформационных системах широко распространена растровая мо дель данных. Растры применяются для хранения и обработки данных дистан ционного зондирования, для представления цифровых моделей рельефа, при визуализации геоданных и т.д. Существует множество вариантов кодирова ния растровых структур. Некоторые из них более экономно расходуют па мять, другие позволяют получать более быстрые алгоритмы. Растровая мо дель соответствует двумерному ячеистому изображению, которое хранится в памяти компьютера в виде одномерной последовательности значений. Рас тровые изображения обычно разлагаются по строке сверху – слева. Далее бу дут описаны другие способы эффективного представления растров.

В некоторых форматах графических файлов используется сжатие изо бражения, основанное на замене длительных последовательностей повто ряющихся значений парой <значение, количество повторов> (рис. 1-а). Гео графические данные обычно автокоррелированны. В растровой модели это означает, что соседние ячейки имеют большую вероятность быть одинако выми, чем разобщенные. При обычном порядке сканирования в конце каж дой строки происходит скачок на начало следующей строки. Предложим простое изменение порядка сканирования. Нечетные строки будем кодиро вать слева направо, а четные – в обратном направлении (рис. 1-б). Направле ние сканирования напоминает движение быка, вспахивающего поле. Отсюда название этого способа сканирования – Boustrophedon (греч. – бык, вспахи вающий поле). Теперь при переходе к новой строке первая ячейка является смежной последней ячейке старой строки. Так в линейном разложении рас тра сохраняется автокорреляция и повышается эффективность кодирования.

A A A A A A A A A B B B а) A B B B A A B B б) A A B B A A A B A A A B AAAA BBBA AABB BAAA 16 байт 4A 3B 3A 3B 3A 8 байт AAAA ABBB AABB AAAB 16 байт 5A 3B 2A 2B 3A 1B 12 байт Рис. 1. Порядки сканирования растров, их линейное разложение и сжатие: а) обычный порядок;

б) Boustrophedon.

Порядок сканирования Мортона (назван по имени Гая Мортона, впервые использовавшего этот способ в Canada GIS) основан на иерархическом раз биении карты. В предыдущих способах сканирования учитывалась автокор реляция значений ячеек только по одному направлению (по строке). Геогра фические объекты образуют на растровом изображении пятна. В порядке Мортона предпринимается попытка сканирования ячеек таким образом, что бы охватить линией обхода эти двумерные пятна.

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

На следующем уровне матрица размера 4 x 4 складывается из четырех мат риц размера 2 x 2, расположенных в таком же порядке, как ячейки матрицы 2 x 2 (рис. 2). Аналогично формируется линия сканирования любой матрицы порядка 2n. Матрица формируется уровень за уровнем, повторяя один и тот же шаблон размера 2 x 2. При сканировании растра по Мортону линия скани рования представляет собой фрактал. Недостатки сканирования по Мортону очевидны. Во-первых, присутствуют скачки, например, от ячейки 7 к ячейке 8. Во-вторых, таким способом можно кодировать только растры размера, кратного двум.

2 3 10 11 14 15 42 43 46 47 58 59 62 8 9 12 0 1 40 41 44 45 56 57 60 2 3 6 34 35 38 39 50 51 54 0 1 4 32 33 36 37 48 49 52 10 11 14 15 26 27 29 A A A A 8 9 12 13 24 25 28 A B B B 2 3 6 7 18 19 22 A A B B 0 1 4 5 16 17 20 A A A B AAAAABBBABAABBAA 5A 3B 1A 1B 2A 2B 2A Рис. 2. Порядок сканирования растра по Мортону.

Рассмотрим следующий способ сканирования, в котором отсутствуют скачки между ячейками. На рис. 3 ячейки растра сканируется по линии Пеа но. Имеется базовый П-образный шаблон, который поворачивается от уровня к уровню так, чтобы обеспечить непрерывность линии сканирования.

При работе с растровыми данными важной является задача определения местоположения ячейки в последовательном файле по растровым координа там и наоборот. Для обычного порядка сканирования и для Boustrophedon – сканирования получение такого отображения не составляет труда. При ска нировании по Мортону задача усложняется.

Рассмотрим пример. Пусть требуется по растровым координатам ячейки A(2, 3) определить ее номер в последовательности Мортона. Для этого пред ставим координаты A в двоичной системе счисления и на их основе сформи руем число N так, что координаты столбца ячейки A задают нечетные биты N, а координаты строки – четные биты. Получившееся число N=( 1 1 0 1 )2= соответствует позиции ячейки A в последовательности Мортона.

15 12 11 10 21 22 25 26 37 38 41 14 13 8 20 23 24 27 36 39 40 1 2 6 19 18 29 28 35 34 45 0 3 4 16 17 30 31 32 33 46 15 12 11 10 53 52 51 A A A A 14 13 8 9 54 55 50 A B B B 1 2 7 6 57 56 61 A A B B 0 3 4 5 58 59 60 A A A B AAAAABBBBBAAABAA 5A 5B 3A 1B 2A Рис. 3. Порядок сканирования растра Пеано.

Обратная задача решается похожим способом. Пусть ячейка B записана в десятой позиции последовательности Мортона. Представим ее номер в двоичной системе счисления и разделим четные и нечетные биты между рас тровыми координатами столбца и строки ячейки: N=10=(1010)2. Получим растровые координаты ячейки A(3, 0).

1.2. Иерархические структуры данных Рассмотренные выше порядки сканирования растровых изображений дают незначительные различия в компрессии данных. Основное преимуще ство Мортон-сканирования и других иерархических структур данных заклю чается в более быстром доступе к данным. Информация распределена по кар те неравномерно. Увеличение разрешения растрового изображения приводит к увеличению размеров файлов, а уменьшение – к потере информации. Далее пойдет речь об адаптивных методах представления растровых данных с раз ной плотностью информации.

На рис. 4 изображена растровая матрица размера 16 x 16, в которой со держатся 255 значений “A” и одно “B”. Индексируем растр следующим спо собом. Разделим матрицу на четыре подматрицы размера 8 x 8 и нумеруем их 0, 1, 2, 3 в порядке Мортона. Назовем подматрицу гомогенной, если в ней со держатся одинаковые значения. Будем рекурсивно разбивать негомогенные подматрицы до тех пор, пока не достигнем гомогенности всех подматриц.

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

Идея выделения гомогенных блоков растра тождественна кодированию растра по Мортону. Гомогенный блок растра размера m x m при сканирова нии по Мортону соответствует коду .

Рис. 4. Разбиение растра на гомогенные блоки.

Соответствие двумерных растровых координат ячейки и адреса ячейки в последовательном файле похоже на аналогичное преобразование при коди ровании растра по Мортону. Единственное отличие в том, что используется система счисления с основанием четыре. В примере на рис. 4 ячейка “B” имеет код 0311. В двоичной системе счисления N=0311=(00110101)2. Разде лим биты между растровыми координатами и выясним, что ячейка лежит в четвертой строке и седьмом столбце.

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

Позиция: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 Содержание: 2 6 A A A A A A 10 A 14 A A A B A A Уровень: 0 1 1 1 1 2 2 2 2 3 3 3 3 4 4 4 Рис. 5. Кодирование квадродерева.

Как уже отмечалось выше, основное преимущество иерархической орга низации данных в ГИС заключается в пространственном упорядочении ин формации и более быстром ее поиске. Поэтому рассматриваются две задачи ГИС, связанные с индексированием квадродеревьями: первая – поиск всех частей карты с заданным значением и вторая – определение содержимого не которой ячейки.

Обозначим n – число уровней квадродерева (тогда размер растра 2n*2n) и через m – число листьев в дереве. Чтобы найти части карты с некоторым значением B, необходимо проверить каждый лист дерева, что потребует m шагов. Определение значения ячейки происходит путем спуска по квадроде реву до тех пор, пока не будет получен гомогенный блок. В худшем случае, когда ячейка находится на самой вершине дерева (как, например, ячейка B на рис. 4), поиск займет n шагов. Сравним теперь трудоемкости обеих задач на квадродереве с трудоемкостями этих задач при различных вариантах скани рования растра.

Таблица 1.

Трудоемкость алгоритмов при различной организации растров.

Структура данных Поиск частей с задан- Определение значения ным значением ячейки Квадродерево m n Обычный порядок 4 n2 * 1 ** Boustrophedon m *** m **** Мортон m *** m **** Прим. * – проверяется каждая ячейка матрицы;

** – непосредственное вычисление позиции ячейки;

*** – число цепочек примерно соответствует числу листьев;

**** – проверяется каждая цепочка.

Существуют различные модификации квадродеревьев, позволяющие, например, эффективно индексировать трехмерные данные (при этом куб ре курсивно делится на восемь частей). При кодировании глобальных данных в проекции Меркатора, представленных в растровой форме, существует про блема различий в форме и размерах ячеек, что приводит к отклонениям в мо дели. Проблема может быть решена путем представления данных в иерархи ческой форме. Для этого строится глобальная тесселяция: земная поверх ность проецируется на октаэдр, содержащий восемь пронумерованных тре угольников. Далее каждый треугольник делится на четыре треугольника со единением отрезками середин его сторон. Получившаяся модификация квад родерева позволяет получать разрешение 1 метр при уровне вложенности де рева, равном 20.

1.3. Алгоритмы на квадродеревьях Иерархическая организация данных позволяет получать быстрые спосо бы доступа к пространственным данным. Рассмотрим теперь некоторые ал горитмы ГИС на квадродеревьях: вычисление площади, оверлейный алго ритм и алгоритмы определения смежности ячеек.

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

Вычислим, например, на карте 1 (рис. 6) площадь ячеек со значением “A”.

Площадь SA = 1*(Count leaf (00,02,03,32)) + 4*(Count leaf (2)) = 8.

Карта 1 A A A B Карта 2 1 2 2 A A B B 1 1 2 A A B B 1 2 2 A B B B 2 2 2 Квадродерево “AB” Квадродерево “12” Квадродерево “AB” + “12” Рис. 6. Оверлейная операция на квадродереве.

Оверлейная задача на квадродеревьях заключается в совмещении квад родеревьев двух карт и получении нового квадродерева, индексирующего обе карты. Для этого требуется одновременно обойти оба дерева, следуя вет вям, существующим в обоих деревьях. В тех узлах, где у одного из деревьев будет отсутствовать ветвление, значение атрибута переносится на все после дующие подуровни. В результате образуется “широкое” дерево, содержащее оба атрибута (рис. 6).

Многие операции ГИС, работающие с иерархическими структурами данных, требуют наличия способов определения смежности ячеек. Для этого будем использовать представление координат ячеек в системе счисления с основанием 4 и разделять биты так, как это делалось в сканировании растра по Мортону. При этом используется tesseral – арифметика, в которой перенос между разрядами осуществляется через две позиции. Например, разность 1000 – 1 = 0010, а сумма 1 + 0001 = 0100.

Будем различать два случая: когда проверяемые ячейки имеют коды одинаковой и разной длины. Блоки одинакового размера являются смежны ми, если их представления в tesseral – арифметике отличаются на 1 или 10.

Например, блоки 01 и 03 являются смежными, так как 0011 – 0001 = 10. Бло ки 033 и 211 также смежные, так как 001111 + 10 = 100101. Блоки 01 и 30 не являются смежными, так как 1100 – 0001 = 1001.

Для определения смежности блоков разного размера коды приводятся к основанию 2. Далее с кодом большей длины суммируются ±01 и ±10. Из по лучившихся четырех кодов отбрасываются как невозможные все отрицатель ные (по переносу). Оставшиеся коды сдвигом вправо приводятся к длине ко да меньшей длины. Два блока являются смежными, если трансформирован ный и обрезанный код большей длины равен короткому коду.

Например, требуется определить, являются ячейки 02 и 2 смежными.

Приведем коды в двоичной системе: 024 = 00102, 24 = 102. Прибавим к длин ному коду ±01 и ±10. Получим 0010+1=0001, 0010+10=1000, 0010–10=0000.

Разность 0010–1 отрицательная, эта комбинация отбрасывается как невоз можная. Выбрав два старших разряда из оставшихся результатов, получим коды 002=04 и 102=24. Один из получившихся кодов равен короткому коду, поэтому ячейки 2 и 02 являются смежными.

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

Процесс построения индекса для цифровой карты включает следующие шаги. Сначала для каждого объекта базы данных находится наименьший лист квадродерева, полностью включающий объект. Некоторые крупные объекты могут лежать более чем в одном квадранте первого уровня квадро дерева. В этом случае объекты помечаются значением “NULL”. Остальные объекты помечаются кодом включающего листа квадродерева. Затем объек ты сортируются по возрастанию получившегося ключа, а сам индексный файл в свою очередь индексируется обычным способом (рис. 7).

Построенные таким способом индексы используются для поиска объек тов, пересекающих заданный полигон или линию. Для этого определяется минимальный лист квадродерева, включающий заданный объект. Подняв шись из полученного узла до вершины дерева и выполнив обход поддерева, корнем которого является этот узел, получим список листьев дерева, внутри которых объекты могут пересекаться с заданным объектом.

Очевидно, пространственные индексы, построенные на квадродеревьях, более эффективны по сравнению с независимым упорядочиванием объектов по x и y, так как в этом случае учитывается пространственный характер дан ных. Индексирование квадродеревьями наиболее целесообразно для мелких объектов (особенно для точек). Большим объектам обычно соответствуют крупные блоки. Для них часто требуется определять пересечения с другими объектами.

Объект Индекс Река ПАРК Деревня Парк NULL Деревня Колодец Колодец Рис. 7. Индексирование цифровой карты квадродеревом.

Проблему индексации крупных объектов можно решить с использова нием R-деревьев (R – rectangle, прямоугольник), в которых также использу ется концепция минимального вмещающего прямоугольника. Здесь требует ся найти два прямоугольника таких, что внутри них расположено максималь но возможное число объектов. При этом нужно стремиться, чтобы количест во объектов в прямоугольника должно быть приблизительно одинаковым.

Прямоугольники могут пересекаться, но площадь пересечения должна быть настолько малой, насколько возможно. Далее эта процедура рекурсивно по вторяется (рис. 8).

12 9 Рис. 8. Индексирование цифровой карты R–деревом.

2. Алгоритмы вычислительной геометрии В геоинформационных системах сложные алгоритмы анализа часто строятся из простых алгоритмов. Рассмотрим сначала некоторые простые ал горитмы, а далее покажем, как из этих простых алгоритмов строятся слож ные аналитические процедуры.

2.1. Пересечение линий Операция нахождение пересечения линий является одной из базовых в ГИС–анализе. Она используется в оверлейных операциях с полигонами, при соединении и разъединении (merge и dissolve) линий и полигонов. Эта опера ция является базисной при определении нахождения точки в полигоне, при удалении расщепленных полигонов. Поэтому эффективные алгоритмы опре деления пересечения линий важны в любой векторной ГИС.

Рассмотрим простейший пример: требуется определить, пересекается ли отрезок AB (4, 2) – (2, 0) с отрезком CD (0, 4) – (4, 0) и если да, то в какой точке? Для этого нужно найти уравнения прямых AB и CD и решить их со вместно (рис. 9-а). Уравнение прямой y=a+bx может быть найдено по двум точкам, через которые она проходит. Коэффициент наклона прямой b=(y2–y1) / (x2–x1). Используя любую из точек, через которые проходит пря мая, найдем a=yi–bxi. Уравнение первой линии y=x–2, а второй линии y=4–x.

Сложив два уравнения, получим точку пересечения (3, 1).

4 2 0 2 4 02 а) б) Рис. 9. Точка пересечения прямых: а) внутри отрезков;

б) снаружи.

В общем виде две линии,заданные уравнениями y=a1+b1x и y=a2+b2x, пересекаются в точке x = –(a1–a2) / (b1–b2);

y = a1+b1x. Однако таким спосо бом можно найти только точку пересечения непараллельных линий беско нечной длины. Возможно отрезки не пересекаются, а пересекаются продол женные по этим отрезкам прямые (рис. 9-б). Отрезки пересекаются, если для точки пересечения (x, y) и точек A, B, C, D выполнены условия:

(xA– x)(x– xB) >= 0;

(xC– x)(x– xD) >= 0;

(yA– y)(y– yB) >= 0;

(yC– y)(y– yD) >= 0.

Необходимо учитывать специальные случаи. Для вертикальных линий угол наклона b стремится к бесконечности, поэтому точку пересечения ищут особым способом. Если обе линии вертикальные, они не пересекаются. Если вертикальная только одна из линий, то подстановкой решается система урав нений y=const и y=a2+b2x. Невертикальные параллельные линии также вызы вают сбой в работе алгоритма, поэтому перед решением системы уравнений следует проверять b1–b2 на равенство нулю.

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

Сложность алгоритма вычисления пересечения полилиний может быть снижена, если предварительно проверять на пересечение минимальные огра ничивающие прямоугольники полилиний. Эти прямоугольники определяют ся минимальными и максимальными координатами x и y. Две полилинии не пересекаются, если не пересекаются их ограничивающие прямоугольники.

Можно применить этот подход и для определения пересечения отдельных сегментов полилиний. Два отрезка AB и CD не пересекаются, если не пере секаются интервалы (xA, xB) и (xC, xD) или не пересекаются интервалы (yA, yB) и (yC, yD).

Следующий метод, впервые использованный в ГИС ArcInfo, основан на разбиении полилинии на секции, в которых линия монотонно возрастает или убывает по x и по y (рис. 10-а). Разбиение происходит в точках локального минимума или максимума по x или по y. Горизонтальная или вертикальная линия пересекает такую секцию только в одной точке. Это дает возможность уменьшить трудоемкость алгоритма поиска пересечения полилиний. Если для двух секций найдена точка пересечения, не нужно проверять оставшиеся пары точек, т.к. это пересечение единственное при условии, что вторые про изводные в секциях не меняют знак. Это ограничение может быть разрешено либо разбиением секции в критических точках, либо полным перебором пар сегментов для таких секций. Модифицированный таким способом алгоритм в некоторых случаях позволяет получить вычислительную сложность порядка O(n1 + n2).

На рис. 10-б представлены два различных случая пересечения секций. В одном случае секции пересекаются только в одной точке, в другом – в не скольких точках. Определим условия, при которых точка пересечения един ственна. Если две секции одновременно возрастают или убывают по одному направлению, одна из них возрастает, а другая убывает по другому, то поли линии на этих секциях пересекаются не более чем в одной точке. На рис. 10-в серым цветом выделены условия, при которых можно применять вышеопи санный метод оптимизации алгоритма поиска пересечений полилиний.

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

Y Y max(y) X возр.

Y возр.

X возр.

Y убыв.

X max(x) X X убыв.

б) Y убыв.

min(x) X убыв.

Y возр.

x возр. x возр. X убыв. x убыв.

y возр. y убыв. y возр. y убыв.

max(x) x возр.

min(x) X возр.

X возр.

y возр.

Y убыв.

Y возр.

X возр. x возр.

Y возр. y убыв.

min(y) x убыв.

y возр.

x убыв.

y убыв.

а) в) Рис. 10. Оптимизация алгоритма определения пересечения полилиний, основанная на разбиении на монотонные секции: а) разбиение на секции;

б) различные варианты пересечения секций;

в) схема определения единственности точки пересечения секций.

2.2. Операции с полигонами Перейдем теперь к операциям с полигонами, заданными последователь ностью вершин. Рассмотрим задачу определения площади полигона. Чаще всего применяется алгоритм, основанный на разбиении многоугольника на трапеции, ограниченные линией сегмента полигона, перпендикулярами, опущенными из вершин сегмента на ось x, и осью x (рис. 11-а). Для сегмента, соединяющего вершины (xA,yA) и (xB,yB), площадь такой трапеции равна S=(xB– xA)* (yB– yA) / 2.

Вычислим площади трапеций для всех сегментов полигона и просумми руем их. Для сегментов, у которых xi > xi+1, площадь берется отрицательной (рис. 11-б). Следует заметить, что полигон – замкнутая фигура, поэтому нуж но учитывать сегмент, соединяющий последнюю вершину с первой.

а) б) Рис. 11. Вычисление площади полигона: а) исходная фигура;

б) разбиение на трапеции.

Таким способом можно вычислить площади не только для выпуклых многоугольников, но и для вогнутых, а также для полигонов, имеющих ды ры. Алгоритм непригоден для вычисления площадей полигонов, имеющих самопересечения границ. Для полигонов, оцифрованных против часовой стрелки, площадь получается отрицательной. Проблемы возникают также при отрицательных значениях координат y вершин полигона. В таком случае можно либо прибавить к координатам y достаточно большое число, либо опускать перпендикуляры из вершин на прямую y=const, где const меньше самой малой y-координаты в полигоне.

Если используется система координат с большими значениями x и y (на пример, в системе координат Гаусса-Крюгера в районе Красноярска дейст вуют координаты x=6200000;

y=16500000), то при многократном суммирова нии площадей трапеций будет накапливаться вычислительная ошибка. Отно сительная погрешность получается особенно высокой для малых полигонов.

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

В модели “дуга-узел” полигоны формируются из дуг. При этом кодиру ется расположение полигона относительно направления цифрования дуги.

Очевидно, определяемую границей двух смежных полигонов площадь доста точно вычислить один раз. Затем для правых полигонов эта площадь сумми руется со знаком “плюс”, а для левых – со знаком “минус”.

Рассмотрим следующую задачу, часто встречающуюся в процедурах ГИС-анализа. Для заданной точки A(u, v) и полигона P=(xi, yi) i=1..n требуется определить, находится точка внутри полигона или снаружи. Задача может быть решена с использованием топологических свойств полигонов. Из точки A проведем вертикальную линию (xi, yi) – (xi, ) и вычислим количество пе ресечений этой линии с сегментами границы полигона. Если это число не четное, точка лежит внутри полигона. Если число четное – точка лежит вне полигона (рис. 12).

Точка Полигон Число пересечений Рис. 12. Схема определения принадлежности точки полигону.

Вертикальная линия x=u пересекает сегмент с концевыми точками (xi, yi) и (xi+1, yi+1) когда эти точки расположены по разную сторону от вертикальной линии. Уравнение прямой сегмента y=a+bx определяется по его концевым точкам, а пересечение сегмента с вертикальной прямой x=u находится в точ ке (u, a+b*u). Вертикальные сегменты границы являются для этого способа определения принадлежности точки полигону специальным случаем. Когда xi=xi+1 и xi<>u, линия и сегмент границы не пересекаются. Когда xi=xi+1 и xi=u, прибавим к количеству пересечений 0 или 2. Для сохранения детерми нированности алгоритма нужно также проверять, расположена ли точка на границе полигона (a+b*u=v). Иначе алгоритм в разных случаях будет выда вать разные результаты (рис. 13).

2 пересечения 3 пересечения Точка снаружи Точка внутри Рис. 13. Неопределенность при расположении точки на границе объекта.

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

Перейдем к задаче определения центральной, репрезентативной точки полигона. В ГИС с этой целью часто используется понятие центроида – точ ки, являющейся центром тяжести полигона. Как видно из рис. 14, центроид расположен не всегда внутри полигона. Однако эта операция может приме няться, например, в картографической генерализации при замене площадных объектов точечными.

Рис. 14. Расположение центроидов полигонов.

Координаты центроида региона с площадью A вычисляются по следую щим формулам: xц= ((yi – yi+1) (xi2 + xi xi+1 + xi+12) / 6A);

yц= ((xi – xi+1) (yi + yi yi+1 + yi+12) / 6A). В некоторых программах координаты центроида вычис ляются как среднее значение x и y.

Рассмотрим следующую полигональную операцию – вычисление скеле та полигона. Скелет полигона определяется как сеть линий, построенная та ким образом, что каждая точка сети расположена на равном расстоянии от ближайших двух сегментов границы полигона. При помощи этой операции можно определить оптимальные места для подписи полигона. Скелет полу чается путем “сжатия” полигона (рис. 15). Сегменты границы полигона сдви гаются внутрь полигона на равные расстояния, поэтому можно считать эту операцию обратной построению буфера.

Рис. 15. Схема получения скелета полигона.

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

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

Оверлеи полигонов изоморфны операциям теории множеств. Когда на кладываются два полигона A и B, получается графическая интерпретация объединения или пересечения множеств A и B. На рис. 16 показаны шестна дцать возможных оверлейных операций с двумя полигонами, выраженные через объединение, пересечение и отрицание помеченных цифрами 1 … множеств.

На основе оверлейных операций строятся некоторые другие функции ГИС. При визуализации данных интерес представляют только объекты, по падающие в “окно” пользователя, а остальные объекты для ускорения ото бражения должны быть пропущены. Для этого на слои карты накладывается прямоугольник – экстент карты, вне которого объекты не отображаются. При построении буферов вокруг точек, сегментов полилиний и полигонов созда ются круги и прямоугольники, которые впоследствии оверлейной операцией сливаются в один объект.

Оверлейные операции применяются при площадной интерполяции.

Здесь требуется распределить некоторую величину, связанную с полигоном A, между пересечением A B и разностью A – B пропорционально их пло щади. При этом считается, что плотность этой величины по всему полигону постоянна.

1 ¬(AB) 2 A ¬B 3 A B 4 B¬A Рис. 16. Связь оверлеев полигонов с операциями теории множеств.

Перейдем теперь к способам реализации оверлейных операций. Будем рассматривать наиболее распространенный в ГИС-анализе случай, когда на кладываются два слоя с непересекающимися полигонами. Представим, что в одном из слоев содержатся “красные” полигоны, а в другом – “синие”. Тогда задача заключается в поиске полигонов на комбинированном “фиолетовом” слое. Атрибуты этого слоя содержат конкатенацию характеристик “синих” и “красных” полигонов. Количество полигонов, получившихся в результате наложения слоев, заранее предсказать нельзя. Из пересекающихся “синего” и “красного” полигонов может получиться сколь угодно много “фиолетовых”.

Чтобы получить оверлей двух полигонов, вначале необходимо вычис лить все пересечения между их границами. На рис. 17 изображен “красный” полигон с атрибутом “A” (тонкие линии) и синий полигон с атрибутом “1” (толстые линии). Внешняя часть на обеих картах имеет атрибут “0”. Каждый полигон представлен одной дугой, для каждой из них известно, с какой сто роны расположен полигон. После вычисления пересечений этих дуг образу ются шесть новых дуг и четыре новых полигона, наследовавших атрибуты 00, A0, A1, 01. Для новых дуг также известно, какие полигоны лежат справа и слева.

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

Дуга Справа Слева в) a) б) Рис. 17. Оверлей полигонов в модели “дуга–узел”: а) исходные объекты;

б) вычисление пересечений дуг;

в) метки смежности дуг и полигонов.

Рассмотрим более сложный пример (рис. 18). Здесь накладываются слой с объектом “1” и слой с тремя объектами “A”, “B” и “C”. Вычислим пересе чения дуг объектов и получим метки правых и левых полигонов получив шихся новых дуг. Как видно из рисунка, дуги 3, 6 и 8 не имеют пересечений с объектом “1”. Определим метки для третьей дуги. Легко видеть, что метки внутри полигона могут передаваться от дуги к дуге. Правый полигон третьей дуги тот же, что правые полигоны второй и четвертой дуги – “А1”. Левый полигон шестой дуги тот же, что левые полигоны второй и четвертой дуги – “А1”. Для третьей дуги, как для части сети “красного” слоя, известна “крас ная” часть метки левого полигона – это “B”. “Синяя” часть левой метки бе рется из метки правого полигона. В результате получается метка “B1”. Ана логично, метка правого полигона для шестой дуги будет “B1”.

Рис. 18. Оверлей полигонов с непересекающимися границами.

Восьмая дуга не пересекается с границей полигона “1” и не является смежной с другими дугами сети “красного” слоя. Для изолированной дуги при помощи алгоритма “точка в полигоне” следует определить вмещающий полигон “синего” слоя. Получим для восьмой дуги правый полигон “С1” и левый полигон “01”. Последний шаг оверлейного алгоритма заключается в формировании полигонов из новых дуг путем обхода полигона от дуги к дуге до тех пор, пока полигон не замкнется.

Точность представления координат сегментов в машинной форме более высока, чем погрешности оцифровки и векторизации. Поэтому при поиске пересечений сегментов полигона могут возникать ошибки, связанные с от сутствием сведений о топологической структуре объектов. На рис. 19-а и на рис 19-б показаны два различных случая пересечения линий.

а) Расщепленные полигоны б) в) Рис. 19. Проблемы поиска пересечения полигонов: а) пересекающиеся сегменты;

б) ложные пересечения сегментов смежных полигонов;

в) расщепленные полигоны.

Из рисунков видно, что в первом случае линии пересекаются на самом деле, а во втором случае пересекающиеся участки линий представляют одну и ту же границу. Необходимо составлять оверлейные алгоритмы таким обра зом, чтобы различать эти ситуации. Полигоны, образующиеся при оверлее двух полигонов с ошибочно векторизованными общими границами, называ ются расщепленными. Для двух полигонов, состоящих из n1 и n2 точек, мо жет образоваться до (2n1n2/(n1+n2) – 3) расщепленных полигонов.

Расщепленные полигоны могут быть устранены либо в процессе овер лейной операции, либо после ее выполнения. В большинстве коммерческих ГИС используется первый подход, заключающийся в “нечетком” представ лении линий. При этом для каждой линии задается уровень толерантности, связанный с возникающей из-за ошибок векторизации неопределенностью геометрии линии. Поиск пересечений ведется для “полос”, заданных самой линией и уровнем толерантности (рис. 20-а). Следует заметить, что опреде ление пересечений для нечетких линий не является транзитивным.

A B а) б) Рис. 20. Удаление расщепленных полигонов: а) до оверлейной операции;

б) после оверлейной операции.

Для устранения расщепленных полигонов после оверлейной операции необходимо определить критерии, по которым расщепленный полигон мож но отличить от настоящего. Расщепленные полигоны обычно имеют неболь шую площадь и вытянутую форму. Они чаще всего состоят из двух дуг. Рас щепленные полигоны характеризуются также “перемежающимися” атрибу тами. Если синяя дуга с атрибутами “1” и “2” накладывается на “красную” дугу с атрибутами “A” и “B”, расщепленные полигоны будут иметь атрибуты “B1” и “A2” (рис. 20-б). Дуги расщепленных полигонов заканчиваются в че тырехвалентных узлах, а в реальных полигонах валентность узлов обычно равна трем. Состоящий из двух дуг расщепленный полигон можно заменить одной дугой, проходящей через центр полигона.

3. Моделирование поверхностей В отличие от цифровых представлений точечных, линейных и площад ных объектов, трехмерные объекты требуют особых форм представления, т.к.

их местоположение описывается не только двумерными, но и высотными ко ординатами. К наиболее распространенному типу трехмерных объектов от носится топографический рельеф земной поверхности. При помощи трех мерных объектов могут быть также смоделированы карты плотности населе ния, атмосферного давления, влажности и т.п. Однако трехмерные модели традиционно связывают с цифровыми моделями рельефа (digital elevation model - DEM). Далее мы будем рассматривать модели рельефа, подразумевая возможность моделирования и других непрерывных феноменов и явлений.

Цифровые модели рельефа позволяют по конечному набору выборочных точек определять возвышение, крутизну склона, направление ската в произ вольной точке на местности. Возможно выявление особенностей местности – бассейнов рек, дренажных сетей, пиков, впадин и т.п. Такие модели широко применяются во многих процедурах ГИС-анализа: при выборе места строи тельства зданий и коммуникаций, в анализе дренажных сетей, в анализе ви димости, при выборе маршрута движении по пересеченной местности. Осо бенно широко цифровые модели рельефа применяются в гидрологии.

Поверхности являются непрерывными феноменами в противополож ность дискретным объектам, выражаемым точками, линиями и полигонами.

Но существуют способы представления поверхностей, в которых использует ся конечное количество точек. Разные подходы к выбору узловых точек, в которых известно значение возвышения поверхности, определяют две наибо лее распространенные модели данных. В геоинформационных системах по верхности обычно описываются при помощи растровых моделей и триангу ляционных сетей. В растровых моделях выборочные точки расположены в узлах регулярной растровой решетки, а в триангуляционных сетях – распола гаются нерегулярно так, чтобы наилучшим образом “обогнуть” поверхность (отсюда название – triangulated irregular networks – TIN).

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

В обеих моделях могут быть вычислены производные к поверхностям, из которых наиболее часто используются угол и экспозиция склона. Угол на клона поверхности в некоторой точке обычно измеряют в градусах или про центах. Плоские регионы имеют нулевую крутизну склонов. Чем круче горы, тем больше угол наклона. Экспозиция склона характеризует направление наибольшего угла наклона в некоторой точке поверхности. Следует помнить, что подобные измерения имеют смысл только для прямоугольных систем ко ординат, а для системы координат широта/долгота результаты будут неточны (особенно в северных широтах).

Исходными данными для создания цифровых моделей рельефа могут быть данные полевых измерений, изолинии рельефа топографических карт, профили фотограмметрических стерео–изображений, гидрографические и гипсометрические карты и т.п. Обратная операция также возможна. На осно ве растровых DEM или триангуляционных сетей можно построить карты рельефа в изолиниях, рассчитать профили поверхности между заданными точками. Эффектным средством визуализации топографической поверхности являются трехмерные изображения.

3.1. Растровые цифровые модели местности В случае, когда выборочные точки располагаются в узлах регулярной решетки, цифровая модель рельефа может быть построена при помощи рас тровой модели. Как известно, эта модель имеет преимущества перед объект ными моделями в простоте алгоритмов обработки и анализа данных, обу словленной простотой организации данных. Растровые DEM являются са мым простым способом представления топографических данных и широко распространены.

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

Z z2 z z3 Y z Z ( x, y) Z (x, y) z z z1 z Y z Z X1 X а) б) в) Рис. 21. Оценка возвышения в произвольной точке:

а) – по ближайшей узловой точке;

б) – аппроксимация МНК;

в) – билинейная интерполяция Более гладкую поверхность можно получить, если аппроксимировать значения высоты в области, ограниченной четырьмя точками сети (рис. 21-б).

При этом необходимо учитывать, что полученная методом наименьших квадратов поверхность не обязательно проходит через узлы решетки, следо вательно, в полученной поверхности вдоль соединяющих узлы линий будут разрывы. Будем искать приближение в окрестности узлов в виде плоскости Z (x, y) = c0 + c1x + c2y. Коэффициенты c0 + c1x + c2y находятся из СЛАУ.

1 x y z c i i i i x x xy c = zx.

i i i i c y xy y zy i i i i Поверхность без разрывов (с разрывом первой производной) можно по лучить, используя билинейную интерполяцию (рис. 21-в). Выберем такую систему координат, что x1 = x2, y2 = y3, x3 – x2 = 1, y2 – y1 = 1. Найдем возвы шения Z5 = Z2 + (Z3 – Z2) x и Z6 = Z1 + (Z4 – Z1) x. Тогда Z = Z6 + (Z5 – Z6) y.

Наклон и направление в произвольной точке растровой DEM вычисля ются с использованием соседних точек. Обычно используется окно 3 x 3. На клон поверхности определяется как отношение изменения возвышения к из менению горизонтального местоположения и выражается в процентах или градусах (рис. 22-а). Наклон измеряется в направлении наиболее крутого из менения возвышения. Чаще всего это направление не совпадает с направле нием строк и столбцов растра (рис. 22-б). Для вычисления угла наклона бу дем использовать формулу Sград = arctan [(Z/ x)2 + (Z/ y)2]1/2. Направле ние склона (aspect) определяется как A = arctan [– (Z/ y) / (Z/ x)].

h S% = h/ l 100%;

Sград = arctan(h / l).

l б) а) Рис. 22. Вычисление наклона поверхности:

а) – вычисление угла наклона поверхности;

б) – направление наиболее крутого изменения возвышения.

Рассмотрим следующий способ определения угла наклона поверхности в растровой ячейке DEM. Вычислим отношение Z / x по значениям ячеек Z и Z5, а Z / y – по ячейкам Z4 и Z5 (рис. 3-а). В этом примере расстояния между центроидами ячеек равно десяти, поэтому Z/x = (49 – 40) / 20 = 0.45;

Z/y = (45 – 48) / 20 = –0.15. Отсюда угол наклона поверхности Sград = arctan [(0.45)2 + (–0.15)2]1/2 = 25.38°. Вычислим направление склона:

A = arctan [ – (–0.15) / (0.45) ] = 18.43°.

Ядро Z / x Ядро Z / y б) а) Рис. 23. Вычисление угла наклона поверхности в растровой ячейке по четырем соседним: а) – схема вычисления;

б) –ядра преобразования.

Вычислить угол наклона поверхности во всех ячейках растрового слоя можно при помощи двух преобразований, ядра которых приведены на рис. 23-б. Эти преобразования позволяют получить для ячейки значения Z/x и Z/y, по которым вычисляется угол наклона. Исходный растро вый слой m x n точек обрабатывается скользящим окном размера 3 x 3. В ре зультате этой фокальной операции получается растровый слой, в ячейках ко торого содержатся значения угла наклона, и имеющий размер (m – 2) (n – 2).

Как видно из рис. 23-а, при определении угла наклона поверхности в растровой ячейке DEM по четырем соседним не используются угловые ячей ки окна – z1, z3, z6, z8. Рассмотрим способ вычисления угла наклона поверх ности по конечным разностям третьего порядка (рис. 24-а).

Y Ядро Z / x Ядро Z / y X б) а) Рис. 24. Вычисление угла наклона поверхности в ячейке по конечным разностям 3-го порядка: а) – схема вычисления;

б) –ядра преобразования.

Здесь Z/x = [1*(47–42) + 2*(49–40) + 1*(52–44)] / 80 = 0.39;

Z/y = [1*(47–52) + 2*(45–48) + 1*(42–44)] / 80 = –0.16. Отсюда угол наклона по верхности Sград = arctan [(0.39)2 + (–0.16)2]1/2 = 22.9°, а направление склона А = arctan [ – (–0.16) / (0.39) ] = 22.36°.

По растровому слою, содержащему углы наклона поверхности, может быть построена тематическая карта. Для этого задается легенда – шкала уг лов наклона с соответствующими цветами отображения. Пример такой тема тической карты приведен на рис. 25-а. На физических картах с целью улуч шения восприятия изображения рельеф показывается с отмывкой – пластиче ским полутоновым изображением рельефа путем наложения теней. На циф ровых картах такой эффект можно получить, раскрашивая ячейки растровой DEM в соответствии с экспозицией склонов. Если источник освещения рас положен на северо-западе, светлыми цветами раскрашиваются ячейки с экс позицией 270°–360°, темными цветами – 90°–180°, средними цветами – 0°–90° и 180°–270° (рис. 25-б).

Угол а) б) Рис. 25. Карты рельефа: а) – тематическая карта углов наклона топографической поверхности;

б) – изображение рельефа с отмывкой по экспозиции.

Растровые DEM могут использоваться в гидрологии для построения мо делей стока, определения водоразделов рек (рис. 26). Направление стока из ячейки растра определяется ее высотой и высотой соседних ячеек. При этом могут рассматриваться все восемь соседей (движение королевы) или четыре соседа по вертикали и горизонтали (движение ладьи). Будем считать, что из данной ячейки имеется сток на соседнюю, если высота соседней ячейки меньше высоты других соседних ячеек и высоты данной ячейки.

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

Направление стока Водораздел бассейна реки Дренажная сеть б) а) Рис. 26. Использование растровых DEM в гидрологии: а) - направление стока, водоразделы и дренажная сеть, рассчитанные на растровой сетке;

б) – тематическая карта дренажной сети и водоразделов.

Пусть наблюдатель находится в ячейке с высотой 100 м. Справа от на блюдателя расположены ячейки с высотами { 110;

125;

115;

160;

… }. Легко видеть, что ячейки с высотами 110, 125, 160 будут наблюдателю видны, а ячейка с высотой 115 – нет. На рис. 27-б показана трехмерная модель с рас считанными для заданного положения наблюдателя зонами видимости.

Видимая часть Невидимая часть Видимая часть Наблюдатель Взгляд Наблюдатель а) б) Рис. 27.Анализ видимости: а) – схема определения видимости;

б) – пример зон видимости для заданного положения наблюдателя (белый цвет).

3.2. Нерегулярные триангуляционные сети (TIN) Нерегулярные триангуляционные сети (Triangulation Irregular Network – TIN) являются альтернативой растровым DEM и используются во многих геоинформационных системах, системах автоматизированного картографи рования, пакетах построения контуров. Модели TIN разработаны в 1970-х годах как простой способ построения поверхностей по нерегулярно располо женным точкам.

Модель TIN обладает некоторыми преимуществами перед растровыми DEM. В первую очередь, расположение точек адаптировано к местности: в равнинных участках точки расположены реже, а гористых – чаще. Выбороч ные точки соединяются прямыми отрезками, образующими треугольники, внутри которых поверхность задается плоскостью. Поверхность непрерывна, треугольники соединены между собой. Структуры данных в TIN-моделях бо лее компактны и экономичны: TIN-модели из сотен точек может соответст вовать растровая DEM из десятков тысяч точек. Несмотря не простоту моде ли, создание TIN требует решения ряда сложных задач: как размещать выбо рочные точки, как соединять точки в треугольники, как моделировать по верхность внутри треугольника.

Рассмотрим задачу выбора размещения точек триангуляции на следую щем примере: имеется растровая DEM или оцифрованные изолинии рельефа, необходимо выбрать точки таким образом, чтобы наиболее точно предста вить поверхность в TIN-модели. Рассмотрим алгоритмы выбора точек DEM.

Алгоритм Фаулера – Литтла основан на поиске особых точек поверхно сти – пиков, хребтов, впадин и т.п. Поверхность проверяется скользящим ок ном размера 3 x 3. При этом соседи центральной ячейки помечаются плюсом, если их высота больше, и минусом если меньше. Очевидно, точка является пиком, если все восемь ее соседних ячеек помечены минусом. Аналогично, точка является впадиной, если все восемь ее соседних ячеек помечены плю сом. Точка является ущельем или перевалом, если плюсы и минусы вокруг точки образуют хотя бы два полных цикла:

{ + + – – – – + + } = 2 цикла;

{ + – + – + – + – } = 4 цикла.

Далее слой обрабатывается окном 2 x 2 таким образом, что каждая ячейка по очереди становится во все позиции окна. Точка является гребнем горы, если в каждом из четырех окон она не самая низкая. Аналогично, точка принадлежит протоку, если во всех четырех окнах она не самая высокая. За тем, начиная от ущелий или впадин, ищутся связные протоки, пока не будут достигнуты пики (поиск можно вести и в обратном направлении). В резуль тате получаются соединенные линиями пики, протоки, впадины, ущелья и перевалы. По выбранным точкам создаются треугольники.

Алгоритм VIP (очень важная точка) в отличие от предыдущего алгорит ма, в котором идентифицируются основные особенности местности, прове ряет поверхность локально, используя окно 3 x 3. Это упрощение впервые использовано в ГИС ESRI Arc/Info. Ячейка в растровой DEM имеет восемь соседей, образующих четыре тройки (рис. 28-а).

Тройка Тройка 1: { 8 – 6 – 4 } Тройка 3: { 10 – 6 – 8 } 10 9 10 6 6 Тройка 6 V = 0 V > 4 5 2 Тройка Тройка б) а) Рис. 28. Алгоритм VIP: а) – четыре тройки ячеек;

б) –расчет вариаций.

Для каждой тройки ячеек по соответствующей вариограмме рассчитыва ется коэффициент вариации (рис. 8-б). Первая тройка имеет нулевой коэф фициент вариации, вторая и четвертая – низкий, а третья – высокий коэффи циент вариации. Далее оценивается средняя вариация значения узла растро вой DEM. Узлы с высокими показателями вариации включаются в результи рующее разбиение, остальные отбрасываются.

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

По узлам регулярной сети легко построить исходное разбиение на тре угольники. В начале работы алгоритма разбиение полностью соответствует исходной растровой DEM. Далее все точки разбиения поочередно проверя ются следующим способом. Точка временно удаляется, соответственно из меняется и разбиение. Затем определяются треугольники, содержащие уда ляемую точку, оценивается разность возвышений этой точки, полученных из DEM и из трех верши треугольника. Эта разность записывается в базе дан ных и удаленная точка восстанавливается. Когда таким образом будут обра ботаны все точки, нужно удалить точки с наименьшими значениями запом ненных в базе данных разностей.

После того, как выбрано необходимое количество узлов TIN, нужно вы брать способ разбиения поверхности на треугольники. При этом желательно получить близкие к равносторонним треугольники, чтобы произвольная точ ка поверхности была как можно ближе к узлам TIN, где значения возвыше ния известны точно. Рассмотрим два способа разбиения на треугольники.

Триангуляция может быть получена путем упорядочивания точек по расстоянию между ними. Для этого вычисляются и сортируются по возраста нию расстояния между всеми парами точек. Ближайшие пары точек соеди няются линией, если эта линия не пересекает полученные на предыдущих шагах линии. Процесс завершается, когда невозможно создать ни одной но вой линии. В результате получится TIN-разбиение, в котором будет много остроугольных треугольников.

Этого недостатка лишена триангуляция Делоне. По определению три точки формируют треугольник в триангуляции Делоне тогда и только тогда, когда в окружности, описанной вокруг этого треугольника нет других точек разбиения. Разобьем поверхность на области, в которых каждая точка распо ложена ближе всего к некоторому узлу сети – генерирующей точке. Полу ченные границы называют полигонами Тиссена или полигонами Вороного.

Две точки соединяются линией в триангуляции Делоне, если их полигоны Тиссена имеют общую границу. Этот метод позволяет получить требуемые “жирные” треугольники.

Это треугольник Добавляемая Выпуклый нужно удалить точка многоугольник а) б) в) г) Рис. 29. Алгоритм Уотсона: а) – исходное разбиение;

б) – выбор треугольников, описанная вокруг которых окружность содержит новую точку;

в) – удаление треугольников;

г) – разбиение выпуклого многоугольника.

Триангуляция Делоне может быть получена при помощи алгоритма Уотсона. В начале работы этого алгоритма создается супертреугольник, со держащий все точки разбиения. Точки последовательно добавляются в суще ствующее разбиение. Опишем процедуру образования нового разбиения при добавлении новой точки (рис. 29). Сначала выбираются треугольники, опи санная вокруг которых окружность содержит добавляемую точку. По опре делению эти треугольники не могут входить в триангуляцию Делоне, поэто му их следует удалить из разбиения. Выбранные треугольники разбиваются на отрезки, дублирующиеся отрезки удаляются. Оставшиеся отрезки форми руют границу выпуклого многоугольника, который разбивается на новые треугольники простым соединением вершин с добавляемой точкой. По окон чании работы алгоритма супертреугольник удаляется.

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

30-а). Во втором способе (рис. 30-б) для каждой точки разбиения сохраняется ее уникальный номер, координаты и список точек, с которыми она соединена прямыми (по часовой стрелке).

9 5 4 ID треугольника 10 ID Координаты (X 10, 1 ;

Y 10, 1 ;

Z 10, 1 ) Координаты (X 8 ;

Y 8 ;

Z 8 ) точек (X 10, 2 ;

Y 10, 2 ;

Z 10, 2 ) точки (X 10, 3 ;

Y 10, 3 ;

Z 10, 3 ) Соседи 6, 7, 9, 4, Сосед №1 Сосед №2 б) Сосед №3 а) Рис. 30. Модели хранения TIN: а) – по треугольникам;

б) – по точкам.

Рассмотрим способ вычисления наклона поверхности и экспозиции склона в TIN-модели. Для этого вычислим нормали к каждому треугольнику разбиения. Треугольник задается тремя точками (xa,ya), (xb,yb), (xc,yc). Тогда нормаль P к плоскости треугольника может быть выражена через векторное произведение двух его сторон как ( yb - ya ) (zc - za ) - (zb - za ) ( yc - ya ) pa P = (zb - za ) (xc - xa ) - (xb - xa ) (zc - za ) = pb.

- xa ) ( yc - ya ) - ( yb - ya ) (xc - xa ) pc (xb Теперь вычислим направляющие косинусы и проекцию нормали P к го ризонтальной плоскости:

pa pb pc pa pb g = ;

h = ;

i = ;

e = ;

f =.

2 2 2 2 2 2 2 2 2 2 2 2 pa + pb + pc pa + pb + pc pa + pb + pc pa + pb pa + pb После этого можно вычислить угол наклона поверхности треугольника и экспозицию склона к северному направлению (рис. 31):

pc ;

= arccos( f ) = arccos pb = arccos(i) = arccos.

2 2 2 2 pa + pb + pc pa + pb Z (0, 0, 1) C P (g, h, i) Y (0, 1, 0) B (e, f, 0) A X (1, 0, 0) Рис. 31. Графическое представление положения угла наклона и экспозиции склона треугольника.

Возвышение произвольной точки внутри треугольника определяется по уравнению плоскости, заданной вершинами треугольника. Плоскость с нор мальным вектором (мы его уже вычисляли) P = {pa, pb, pc}, проходящая через точку M0(x0, y0, z0), описывается уравнением pa ( x – x0) + pb ( y – y0) + pc ( z – z0) = 0. Отсюда по известным значениям x и y находятся возвышения произвольных точек.

Анализ стока в триангуляционных сетях осуществляется на основе рас считанных для каждого треугольника угла наклона и экспозиции склона. При этом можно выделить два подхода к моделированию стока. В первом из них с треугольниками сети обращаются как с дискретными элементами: вода сте кает с одного треугольника на другой в направлении максимального наклона.

Второй способ предполагает работу с поверхностью как с мозаикой плоско стей треугольников. Здесь вода течет в треугольнике, собирается и движется вдоль дуг.

По растровым DEM и по триангуляционным сетям можно построить карты распределения характеристики по территории. Для этого нужно задать список значений z, по которым будут проведены изолинии. Например, можно взять максимальное и минимальное значения характеристики, разделить по лученный интервал на несколько частей и с полученным шагом провести изолинии.

Изолинии заданной высоты строятся следующим способом. Рассмотрим все ячейки растровой DEM (или все треугольники триангуляционной сети).

Будем искать пересечения ячеек (треугольников) с горизонтальной плоско стью Z=const. Из полученных отрезков нужно собрать полилинии. Для этого следует объединять смежные отрезки. Полилинии могут быть сглажены, на пример, при помощи сплайн–функций (рис. 32).

а) б) Рис. 32. Построение изолиний в растровой модели рельефа (а) и в триангуляционной сети (б).

Цифровые модели местности могут использоваться для построения про филей местности между выбранными точками. Для решения этой задачи нужно найти точки пересечения линии профиля с ячейками DEM или тре угольниками TIN и оценить возвышение поверхности в полученных точках.

Полученные точки соединяются и образуют полилинию: вдоль горизонталь ной оси показывается протяженность профиля, а вдоль вертикальной – воз вышения (рис. 33).

Расстояние, мили Рис. 33. Профиль топографической поверхности.

Возвышение, футы Литература 1. Варфоломеев И.В., Савельев А.С. Представление и обработка про странственных данных в ГИС: Методические указания / Красноярск:

КГТУ, 2001.

2. Коновалова Н.В., Капралов Е.Г. Введение в ГИС. Петрозаводск:

Изд-во Петрозаводского Госуниверситета, 1995.

3. Кошкарев А.В. Толковый мини-словарь основных терминов по гео информатике (с английскими эквивалентами) // ГИС-обозрение, 1994.

№1, с59-62.

4. Кошкарев А.В., Тикунов В.С. Геоинформатика / Под ред.

Д.В. Лисицкого. М.: Картгеоцентр - Геоиздат, 1993.

5. Тикунов В.С. Моделирование в картографии: Учеб. М.: Изд-во МГУ, 1997.

6. Цикритзис Д., Лоховски Ф. Модели данных. М. Финансы и статистика, 1985.

7. Якубайлик О.Э Методы и приемы пространственного анализа в гео информационных системах: Учебное пособие / Красноярск: Изда тельство КрасГУ, 2001.

Содержание 1. Структуры пространственных данных ГИС........................................... 1.1. Хранение растровых данных.............................................................. 1.2. Иерархические структуры данных..................................................... 1.3. Алгоритмы на квадродеревьях......................................................... 1.4. Пространственные индексы.............................................................. 2. Алгоритмы вычислительной геометрии................................................. 2.1. Пересечение линий............................................................................ 2.2. Операции с полигонами.................................................................... 2.3. Оверлей полигонов............................................................................ 3. Моделирование поверхностей................................................................. 3.1. Растровые цифровые модели местности......................................... 3.2. Нерегулярные триангуляционные сети (TIN)................................. Содержание....................................................................................................




© 2011 www.dissers.ru - «Бесплатная электронная библиотека»

Материалы этого сайта размещены для ознакомления, все права принадлежат их авторам.
Если Вы не согласны с тем, что Ваш материал размещён на этом сайте, пожалуйста, напишите нам, мы в течении 1-2 рабочих дней удалим его.