GIS-Lab: Полиномиальные преобразования

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

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

Математические выкладки описывающие аналитическое решение задачи приводятся в отдельной статье.

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

Сделаем нашей тестовой точкой точку с координатами: x = 500 и y = -300. Используя возможности ERDAS IMAGINE получим его вариант предсказания: 596703.345492103, 5437370.8467262

Так же, используя тот же инструмент ERDAS IMAGINE, посмотрим вычисленные им коэффициенты преобразования. Примечание: по какой-то причине, ERDAS IMAGINE показыват коэффициенты не прямой, а обратной задачи, поэтому при вычислении например в Excel, нужно будет поменять местами x,y и x',y'. Итак, для нашего тестового набора данных коэффициенты расчитанные ERDAS IMAGINE следующие:

Для быстрой проверки можно воспользоваться данной таблицей в MS Excel, которая использует вышеприведенные расчеты для вычисления трансформации и предсказания новых координат. С помощью этой таблицы, используем тестовый набор данных и сравним коэффициенты преобразования с полученными в ERDAS IMAGINE:

Результат - не отличается от результата ERDAS IMAGINE.

И решим прямую задачу (x,y -> x2,y2) для нахождения коэффициентов a и b:

Или обратную задачу (x,y <- x2,y2), для этого нам также понадобится переопределить матрицу:

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

Создадим тестовую точку и проверим результат:

Что отличается от результатов ERDAS лишь на доли метра. Что и требовалось доказать.

Файлы представляют из себя лист рассчётов для среды MathCad (версия 11 и выше), где на примере показано как можно трансформировать координаты с импольльзованием полиномиальных преборазований 1-го и 2-го порядков (скачать). Прислал Александр Г.

Для расчетов также понадобится дополнительная библиотека для осуществления операций с матрицами.

map:record
     pnts:array of record                    // набор точек для привязки карты
     xr,yr,xg,yg:extended;                 // xr,yr - растровые координаты; xg,yg - географические координаты
     end;
     coeff_x:array of extended;
     coeff_y:array of extended;
     end;

function getx(x,y:extended):extended;var m:word;
begin
    if length(main.map.coeff_x)>0 then getx:=main.map.coeff_x[1]+main.map.coeff_x[2]*x+main.map.coeff_x[3]*y+main.map.coeff_x[4]*x*x+main.map.coeff_x[5]*x*y+main.map.coeff_x[6]*y*y else getx:=-1;
end;
function gety(x,y:extended):extended;
begin
    if length(main.map.coeff_y)>0 then gety:=main.map.coeff_y[1]+main.map.coeff_y[2]*x+main.map.coeff_y[3]*y+main.map.coeff_y[4]*x*x+main.map.coeff_y[5]*x*y+main.map.coeff_y[6]*y*y else gety:=-1;
end;

procedure map_calculate_ceeff;var aa:TReal2DArray;var i,j:byte;
begin
// вычисление коэффициентов полинома
//--------------------------------------------------------
// заполнение обратной матрицы
 setlength(aa,7,7);
 for i:=1 to 6 do begin
   aa[1,i]:=1;
   aa[2,i]:=main.map.pnts[i-1].xg;
   aa[3,i]:=main.map.pnts[i-1].yg;
   aa[4,i]:=main.map.pnts[i-1].xg*main.map.pnts[i-1].xg;
   aa[5,i]:=main.map.pnts[i-1].xg*main.map.pnts[i-1].yg;
   aa[6,i]:=main.map.pnts[i-1].yg*main.map.pnts[i-1].yg;
 end;
 // вычисление обратной матрицы
 if Inverse(aa,6)=false then begin main.show_err('Ошибка. Необходим другой набор контрольных точек.');exit;end;

 // умножение с обратной матрицей X
 setlength(main.map.coeff_x,7);for j:=1 to 6 do main.map.coeff_x[j]:=0;
 for i:=1 to 6 do for j:=1 to 6 do main.map.coeff_x[i]:=main.map.coeff_x[i]+aa[j,i]*main.map.pnts[j-1].xr;
 // умножение с обратной матрицей Y
 setlength(main.map.coeff_y,7);for j:=1 to 6 do main.map.coeff_y[j]:=0;
 for i:=1 to 6 do for j:=1 to 6 do main.map.coeff_y[i]:=main.map.coeff_y[i]+aa[j,i]*main.map.pnts[j-1].yr;
end;

Полиномиальные преобразования - примеры реализации