Геопривязка карт
Топографические карты на территорию СССР составлялись в проекции Гаусса-Крюгера. Задание проекции позволяет ввести систему плоских прямоугольных координат — соответствующая сетка наносится на лист карты.
Когда вы работаете с листом карты в бумажном виде, можно определить координаты любой точки, проведя несложные измерения относительно линий сетки, измерить расстояния, углы и площади. Использование ГИС значительно упрощает рутинные операции, поскольку все объекты в ГИС-проекте помещены в некоторую систему координат — но для этого растровую карту сначала нужно туда поместить, то есть привязать.
Если забыли, как скачивать и открывать данные в QGIS — см. Практику 1.
Определение целевой системы координат
Привязка растра начинается с определения целевой системы координат. Мы работаем с топографической картой, изданной в 1990 г. в рамках системы топографического картографирования СССР.
Топографические карты СССР составлялись в проекции Гаусса-Крюгера на основе датума СК-42, использующего эллипсоид Красовского, с разграфкой по 6-градусным зонам. Описание целевой системы координат должно включать все эти компоненты.
В правом верхнем углу в окне поиска введите:
Pulkovo 1942 / Gauss-Krugerи нажмите Search.
Примечание
Датум СК-42 в англоязычной литературе обозначается как Pulkovo 1942. Также в описаниях систем координат орфография немецкого языка традиционно не соблюдается: для поиска систем координат проекции Гаусса-Крюгера следует вводить Gauss-Kruger, а не Gauss-Krüger или Gauss-Krueger.
В открывшемся окне используйте кнопку Next Page, чтобы переходить на следующие страницы и найти системы координат, соответствующие нужной вам зоне.
Обратите внимание на три важных момента:
- Для вашего листа карты подходит несколько разных систем координат.
- Название соответствующей системы координат не обязательно содержит номер зоны — вместо него может быть указан средний меридиан зоны (Central Meridian, CM). Центральный меридиан для нужной зоны можно определить самостоятельно, используя знания из курса топографии.
- Среди описаний встречаются системы координат 3-градусных зон. Они иногда используются, но ваша карта составлена на основе 6-градусной разграфки, поэтому 3-градусные зоны вам не подходят.
Найдя нужные системы координат, откройте страницу каждой из них в новой вкладке браузера.
Ответьте на вопрос: чем отличаются системы координат? Чтобы это определить, посмотрите описания систем координат в форматах Proj4 или ESRI WKT — эти форматы дают подробную информацию о параметрах системы координат: проекции (например, Gauss-Kruger), датуме (например, Pulkovo 1942), системе координат (плоские прямоугольные координаты), смещении координат (false easting, или x_0 в Proj4).
Прямоугольные координаты, подписываемые вдоль северной и южной рамок топографических карт СССР, формируются так: берётся фактическое значение прямоугольной координаты, отсчитанное от среднего меридиана зоны, смещённого к востоку на 500 км (например, 348); слева к этому номеру приписывается номер зоны (например, 8); полученное число (8348) обозначает координату линии сетки в километрах. Полностью числа подписываются только в углах карты, в остальных случаях — только последние две цифры (например, 48).
При описании таких координат в ГИС-среде может использоваться полная форма записи, включающая номер зоны, либо сокращённая (англ. truncated), не включающая номер зоны. Отличие между вариантами — в единственном параметре: смещении координаты X (x_0 в Proj4, или false easting в WKT). Удобнее пользоваться полной формой записи — для точности и однозначности работы с системой координат.
Типичная ошибка
В проекции Гаусса-Крюгера ось X традиционно направлена на север, ось Y — на восток. В QGIS и большинстве других геоинформационных программ, напротив, ось X направлена на восток, а ось Y — на север. При ручном вводе координат опорных точек эту пару легко перепутать местами — тогда все точки лягут не туда, а привязка не будет складываться в разумную картинку.
О методах передискретизации
При геопривязке (как и при любом пересчёте растра в новую сетку пикселей — например, при изменении разрешения снимка) значения новых пикселей нужно как-то получить из значений исходных. Простейший способ — взять значение ближайшего пикселя (метод «ближайшего соседа»): быстро, но даёт «ступенчатую» картинку. Более гладкий результат даёт билинейная интерполяция (используются 4 соседних пикселя), ещё более гладкий — кубическая и бикубическая интерполяция (используются 16 соседних пикселей вместо 4).
Для данного упражнения этот выбор не критичен — используйте метод передискретизации «Билинейная (ядро 2×2)», указанный далее в настройках трансформации.
Аффинные преобразования и их связь с геопривязкой
Определив систему координат и способ пересчёта пикселей, разберёмся, что именно происходит, когда QGIS «натягивает» скан карты на местность. Аффинное преобразование — отображение плоскости (или пространства) в себя, при котором параллельные прямые переходят в параллельные, пересекающиеся — в пересекающиеся, скрещивающиеся — в скрещивающиеся. Именно этим преобразованием мы и занимаемся при геопривязке — вот почему её называют ещё геометрической коррекцией.
Геометрическими искажениями изображения называют различия между координатами и размерами реального изображения и идеального, которое теоретически можно было бы получить идеальным сенсором в идеальных условиях съёмки. Когда природа искажений неизвестна или сложна, их компенсируют полиномиальным преобразованием n-го порядка. Простейший случай — полином 1-й степени, то есть аффинное (линейное) преобразование:
\[ u = a_0 + a_1 x + a_2 y, \qquad v = b_0 + b_1 x + b_2 y \]Такое преобразование описывает параллельный перенос, поворот, масштабирование, зеркальное отражение и их комбинации — а также перевод изображения в картографическую проекцию или пересчёт между проекциями (собственно то, чем мы и занимаемся в этой практике).
Рис. П4.а — типы искажений, компенсируемых аффинным (линейным) преобразованием
Не все искажения снимков линейны — оптика, рельеф местности и особенности сенсора могут давать искажения, которые линейное преобразование компенсирует лишь приближённо. В этом случае повышают порядок полинома (квадратичный, кубический), но за это приходится платить бо́льшим числом опорных точек.
Сколько нужно опорных точек?
Коэффициенты преобразования вычисляются по наземным контрольным точкам (НКТ) — парам точек на снимке (исходные координаты) и на карте местности (эталонные координаты). Минимальное число пар НКТ зависит от степени полинома n: N = (n+1)(n+2)/2.
| Преобразование | Минимум пар НКТ | Рекомендуемое число пар |
|---|---|---|
| Линейное (аффинное) | 3 | 6 |
| Квадратичное | 6 | 12 |
| Кубическое | 10 | 20 |
В практике геопривязки листа карты (см. ниже) мы используем 5 точек для линейного преобразования — это больше минимума (3) и близко к рекомендованному числу (6), что и обеспечивает разумный запас точности. Точки для привязки стоит выбирать в местах с хорошо выраженными особенностями (пересечения линий сетки, чёткие ориентиры) и равномерно по всей площади листа, включая его края — скопление точек в одном углу плохо ограничивает решение в остальной части листа.
Как это считается
Коэффициенты a₀…a₂ и b₀…b₂ находятся как решение системы линейных уравнений методом наименьших квадратов: минимизируется суммарная невязка между эталонными координатами НКТ и координатами, которые давало бы преобразование. Сама система линейных уравнений — в матричной форме, где неизвестными выступают коэффициенты преобразования, а известными — координаты НКТ; матричная форма метода наименьших квадратов разбирается отдельно в лекционной части курса.
Привязка листа топографической карты
- Запустите QGIS.
- Сохраните проект: «Файл» → «Сохранить проект», выберите рабочую директорию.
- Запустите инструмент привязки растров: в меню «Слой» выберите «Геопривязка…».
- Включите подписи идентификаторов точек: в окне привязки — «Параметры» → «Параметры геопривязки…» → отметьте «Включить идентификаторы».
- Добавьте растровое изображение в инструмент привязки — перетащите файл из проводника или используйте кнопку
на панели инструментов, затем укажите систему координат исходного растра.
В окне привязки по умолчанию активен инструмент добавления опорных точек «Добавить точку» — каждая точка, которую вы поставите этим инструментом, и есть одна из НКТ (наземных контрольных точек) из раздела «Сколько нужно опорных точек?» выше: пара «пиксель на скане ↔ координата на местности», по которым QGIS вычислит коэффициенты a₀…a₂ и b₀…b₂ преобразования. С помощью мыши увеличьте изображение северо-западного угла карты, найдите пересечение линий сетки и определите его координаты, затем щёлкните по пересечению и введите координаты во всплывающем окне.
Кроме пересечений линий сетки, на топографической карте есть метки географических координат — они подписаны у каждого угла рамки карты и у перекрестия в центре листа. Эти точки тоже можно использовать как опорные; отличие в том, что при вводе географических координат нужно выбрать географическую систему координат — для отечественных топографических карт это Pulkovo 1942 (EPSG:4284).
Аналогичным образом добавьте ещё точки — всего должно получиться пять опорных точек: четыре в углах листа карты и одна в центре. Это больше минимума для линейного преобразования (3) и близко к рекомендованному числу (6) — см. таблицу в разделе «Сколько нужно опорных точек?».
В таблице внизу окна «X источника» и «Y источника» — координаты точек привязки на исходном растре в пикселях (от левого верхнего угла); «X назначения» и «Y назначения» — координаты на местности. В столбцах «dX», «dY» и «Невязка» отображаются несхождения; пока там нули, поскольку параметры трансформации ещё не заданы.
Откройте настройки трансформации растра. Параметр «Тип трансформации» — это и есть выбор степени полинома преобразования (см. раздел «Аффинные преобразования и их связь с геопривязкой» выше): для пяти точек и относительно простых искажений скана подходит линейное преобразование. В списке QGIS есть и другие варианты — ниже показано, как каждый из них по-разному деформирует одну и ту же правильную сетку:
Рис. П4.б — та же исходная сетка (Original) после разных типов трансформации, доступных в диалоге привязки QGIS
Linear и Helmert сохраняют форму прямоугольника (Helmert допускает ещё и поворот); Polynomial 1 — это и есть аффинное преобразование из раздела выше: прямые остаются прямыми и параллельными, но сетка может стать параллелограммом. Polynomial 2/3 и сплайн допускают нелинейный, «резиновый» изгиб сетки — полезно для сильно покоробленных сканов, но требует намного больше опорных точек. Projective сохраняет прямые линии, но не параллельность — им моделируют перспективные искажения (например, снимок, сделанный не строго в надир).
Укажите, куда сохранить привязанный файл (параметр «Целевой растр») — например, в папку referenced_images рабочей директории под именем o38-085_modified.tif. Отметьте опции «Сохранить контрольные точки» и «Открыть результат в QGIS». Остальные параметры оставьте по умолчанию.
Примечание
Окно параметров трансформации при добавлении растра позволяет выбрать набор параметров перехода с одного датума на другой (в нашем случае — со старого отечественного датума СК-42 на стандартный в ГИС датум WGS 84). В большинстве случаев можно использовать вариант по умолчанию; подробное обсуждение преобразования датумов выходит за рамки курса.
Закройте окно настроек привязки — привязка не запустится автоматически, для запуска нужно нажать отдельную кнопку. Но прежде изучите величины ошибок в таблице внизу (в пикселях исходного растра). В практике привязки обычно стараются добиться ошибки не более 0,5 пиксела (с оговорками). Наш картографический материал отсканирован с невысоким разрешением, и такой точности достичь не получится — для учебного упражнения достаточно точности порядка 1,5–2 пикселов. Если точность какой-либо точки значительно хуже, удалите её и установите заново.
Типичная ошибка
Не путайте невязку в пикселях исходного растра с невязкой на местности. Один и тот же пиксель на низком разрешении скана может соответствовать десяткам метров на местности — «хорошая» невязка 1,5–2 пиксела не означает точность привязки 1,5–2 метра.
После запуска привязки окно QGIS примет вид, аналогичный представленному ниже — привязанная карта ложится на своё место поверх спутниковой подложки.
Повышение точности трансформации
Если итоговая невязка вас не устраивает, есть несколько стандартных способов её уменьшить:
- исключить контрольную точку с наибольшей среднеквадратичной невязкой и пересчитать коэффициенты преобразования заново;
- пересмотреть допустимый порог ошибки;
- повысить порядок трансформирования (например, перейти от линейного к квадратичному) — но помните, что это требует больше пар точек (см. таблицу выше);
- оставить только те точки, в которых вы уверены больше всего, — даже если их останется минимально допустимое количество.
Оцифровка контура и сравнение с современными данными
Привязанная карта — уже полноценный слой ГИС-проекта, из которого можно оцифровать (векторизовать вручную) любой линейный объект — например, береговую линию на дату издания карты (1990 год). Сравнив этот исторический контур с контуром, полученным по спутниковому снимку в Практике 3 (2023 год), можно количественно оценить смещение берега за десятилетия.
Чтобы измерить это расхождение по всей длине береговой линии, а не «на глаз», удобно разместить вдоль одной из линий равномерный набор точек и для каждой точки вычислить расстояние до ближайшей точки на другой линии.
1. Точки вдоль линии. Инструмент «Вектор → Геометрия → Точки вдоль геометрии» создаёт слой точек, равномерно распределённых вдоль линий исходного векторного слоя, с шагом, заданным параметром «Расстояние» (выбирайте шаг исходя из длины и детальности вашей линии — например, от 10 до 50 метров).
2. Расстояние до ближайшей точки другой линии. Откройте таблицу атрибутов слоя точек и создайте новое поле через калькулятор полей:
distance(closest_point(geometry(get_feature('map','id','1')), $geometry), $geometry)Разберём выражение по частям: get_feature('map', 'id', '1') находит нужный объект (линию) в слое map; geometry(...) получает его геометрию; closest_point(geometry, $geometry) находит ближайшую к текущей точке ($geometry) точку на этой линии; внешняя функция distance(...) вычисляет расстояние между текущей точкой и найденной ближайшей точкой. В результате каждая точка вдоль одной линии получает атрибут — расстояние до другой линии в единицах проекции (метрах).
Настройте отображение точек так, чтобы размер или цвет маркера зависел от значения поля с расстоянием (data-defined override в свойствах стиля слоя) — это сразу покажет на карте, где расхождение между историческим и современным контуром максимально, то есть где берег отступил или нарастил больше всего.
Типичная ошибка
Функция closest_point находит ближайшую точку геометрически, вдоль всей линии — она не обязана находиться «напротив» исходной точки в интуитивном смысле. На сильно изогнутых участках берега (бухты, косы) такая пара точек может давать длину, которую нельзя напрямую трактовать как перпендикулярное смещение берега. Для таких участков результат стоит проверять визуально, а не только по числу.