Конформное преобразование: различия между версиями
ErnieBoyd (обсуждение | вклад) Нет описания правки |
ArmRus (обсуждение | вклад) м (опечатка названия столбца расчетов) |
||
(не показано 36 промежуточных версий 2 участников) | |||
Строка 1: | Строка 1: | ||
{{Статья| | {{Статья|Опубликована|helmert2d}} | ||
Конформное преобразование на плоскости широко используется в геодезии при создании местных координатных систем на небольшие территории, ограниченные размерами населённого пункта. | {{Аннотация|Конформное преобразование на плоскости широко используется в геодезии при создании местных координатных систем на небольшие территории, ограниченные размерами населённого пункта.}} | ||
== Введение == | == Введение == | ||
Строка 7: | Строка 7: | ||
Следующие формулы связывают координаты точек ''x'', ''y'', заданные в местной системе координат (МСК), и координаты ''X'', ''Y'', заданные в государственной системе координат (ГСК): | Следующие формулы связывают координаты точек ''x'', ''y'', заданные в местной системе координат (МСК), и координаты ''X'', ''Y'', заданные в государственной системе координат (ГСК): | ||
<math>\begin{ | <math>\begin{align} | ||
X & = | X & = X_0 + m \big[ (x - x_0) \cos \theta - (y - y_0) \sin \theta \big] \\ | ||
Y & = | Y & = Y_0 + m \big[ (x - x_0) \sin \theta + (y - y_0) \cos \theta \big] | ||
\end{ | \end{align}</math> | ||
где ''m'' — масштабный множитель, ''θ'' — угол разворота, ''X''<sub>0</sub>, ''Y''<sub>0</sub>, ''x''<sub>0</sub>, ''y''<sub>0</sub> — координаты одного из геодезических пунктов в ГСК и МСК, как правило. Этот набор параметров называется «ключ». | где ''m'' — масштабный множитель, ''θ'' — угол разворота, ''X''<sub>0</sub>, ''Y''<sub>0</sub>, ''x''<sub>0</sub>, ''y''<sub>0</sub> — координаты одного из геодезических пунктов в ГСК и МСК, как правило. Этот набор параметров называется «ключ». | ||
Строка 20: | Строка 20: | ||
Конформное преобразование представляется следующей математической моделью: | Конформное преобразование представляется следующей математической моделью: | ||
<math>\begin{ | <math>\begin{align} | ||
X' & = | X' & = a_0 + a_1 x - b_1 y \\ | ||
Y' & = | Y' & = b_0 + b_1 x + a_1 y | ||
\end{ | \end{align}</math> | ||
Определению подлежат четыре параметра: ''a''<sub>0</sub>, ''b''<sub>0</sub>, ''a''<sub>1</sub>, ''b''<sub>1</sub>,. | Определению подлежат четыре параметра: ''a''<sub>0</sub>, ''b''<sub>0</sub>, ''a''<sub>1</sub>, ''b''<sub>1</sub>,. | ||
Строка 65: | Строка 65: | ||
=== Шаг 4: вычисление ''a''<sub>0</sub> и ''b''<sub>0</sub> === | === Шаг 4: вычисление ''a''<sub>0</sub> и ''b''<sub>0</sub> === | ||
<math>\begin{ | <math>\begin{align} | ||
a_0 & = | a_0 & = X_c - a_1 x_c + b_1 y_c \\ | ||
b_0 & = | b_0 & = Y_c - b_1 x_c - a_1 y_c | ||
\end{ | \end{align}</math> | ||
=== Шаг 5: вычисление невязок === | === Шаг 5: вычисление невязок === | ||
Строка 94: | Строка 94: | ||
== Пример вычисления параметров == | == Пример вычисления параметров == | ||
Даны координаты четырёх пунктов: | Даны координаты четырёх пунктов ''x'', ''y'' в исходной системе, ''X'', ''Y'' в конечной системе с весами ''p'': | ||
{| class="wikitable" | {| class="wikitable" | ||
Строка 102: | Строка 102: | ||
| 1 || 1334.71 || 285.94 || 83477.64 || 87377.60 || 1.0 | | 1 || 1334.71 || 285.94 || 83477.64 || 87377.60 || 1.0 | ||
|- align="right" | |- align="right" | ||
| 2 || 563.67 || | | 2 || 563.67 || −5197.34 || 82557.14 || 81916.51 || 1.0 | ||
|- align="right" | |- align="right" | ||
| 3 || 4444.27 || 1153.79 || 86610.19 || 88160.39 || 1.0 | | 3 || 4444.27 || 1153.79 || 86610.19 || 88160.39 || 1.0 | ||
|- align="right" | |- align="right" | ||
| 4 || | | 4 || −252.07 || 2881.90 || 81962.05 || 90016.34 || 1.0 | ||
|- align="right" | |- align="right" | ||
| colspan="5" | ∑''p'' = || 4.0 | | colspan="5" | ∑''p'' = || 4.0 | ||
Строка 121: | Строка 121: | ||
| 1 || 1334.71 || 285.94 || 83477.64 || 87377.60 | | 1 || 1334.71 || 285.94 || 83477.64 || 87377.60 | ||
|- align="right" | |- align="right" | ||
| 2 || 563.67 || | | 2 || 563.67 || −5197.34 || 82557.14 || 81916.51 | ||
|- align="right" | |- align="right" | ||
| 3 || 4444.27 || 1153.79 || 86610.19 || 88160.39 | | 3 || 4444.27 || 1153.79 || 86610.19 || 88160.39 | ||
|- align="right" | |- align="right" | ||
| 4 || | | 4 || −252.07 || 2881.90 || 81962.05 || 90016.34 | ||
|- align="right" | |- align="right" | ||
| ∑ = || 6090.58 || | | ∑ = || 6090.58 || −875.71 || 334607.02 || 347470.84 | ||
|} | |} | ||
Строка 136: | Строка 136: | ||
! ''x<sub>c</sub>'' !! ''y<sub>c</sub>'' !! ''X<sub>c</sub>'' !! ''Y<sub>c</sub>'' | ! ''x<sub>c</sub>'' !! ''y<sub>c</sub>'' !! ''X<sub>c</sub>'' !! ''Y<sub>c</sub>'' | ||
|- align="right" | |- align="right" | ||
| 1522.645 || | | 1522.645 || −218.9275 || 83651.755 || 86867.71 | ||
|} | |} | ||
Строка 145: | Строка 145: | ||
! ''i'' !! ∆''x<sub>i</sub>'' !! ∆''y<sub>i</sub>'' !! ∆''X<sub>i</sub>'' !! ∆''Y<sub>i</sub>'' | ! ''i'' !! ∆''x<sub>i</sub>'' !! ∆''y<sub>i</sub>'' !! ∆''X<sub>i</sub>'' !! ∆''Y<sub>i</sub>'' | ||
|- align="right" | |- align="right" | ||
| 1 || | | 1 || −187.935 || 504.8675 || −174.115 || 509.89 | ||
|- align="right" | |- align="right" | ||
| 2 || | | 2 || −958.975 || −4978.4125 || −1094.615 || −4951.20 | ||
|- align="right" | |- align="right" | ||
| 3 || 2921.625 || 1372.7175 || 2958.435 || 1292.68 | | 3 || 2921.625 || 1372.7175 || 2958.435 || 1292.68 | ||
|- align="right" | |- align="right" | ||
| 4 || | | 4 || −1774.715 || 3100.8275 || −1689.705 || 3148.63 | ||
|} | |} | ||
Строка 161: | Строка 161: | ||
! ''p<sub>i</sub>'' ∆''x<sub>i</sub>'' ∆''Y<sub>i</sub>'' | ! ''p<sub>i</sub>'' ∆''x<sub>i</sub>'' ∆''Y<sub>i</sub>'' | ||
! ''p<sub>i</sub>'' ∆''y<sub>i</sub>'' ∆''X<sub>i</sub>'' | ! ''p<sub>i</sub>'' ∆''y<sub>i</sub>'' ∆''X<sub>i</sub>'' | ||
! ''p<sub>i</sub>'' (∆''x<sub>i</sub>''<sup>2</sup> + ∆'' | ! ''p<sub>i</sub>'' (∆''x<sub>i</sub>''<sup>2</sup> + ∆''y<sub>i</sub>''<sup>2</sup>) | ||
|- align="right" | |- align="right" | ||
| 1 || 32722.3 || 257426.9 || | | 1 || 32722.3 || 257426.9 || −95826.2 || −87905.0 || 290210.8 | ||
|- align="right" | |- align="right" | ||
| 2 || 1049708.4 || 24649116.0 || 4748077.0 || 5449445.0 || 25704224.1 | | 2 || 1049708.4 || 24649116.0 || 4748077.0 || 5449445.0 || 25704224.1 | ||
Строка 169: | Строка 169: | ||
| 3 || 8643437.7 || 1774484.5 || 3776726.2 || 4061095.5 || 10420246.0 | | 3 || 8643437.7 || 1774484.5 || 3776726.2 || 4061095.5 || 10420246.0 | ||
|- align="right" | |- align="right" | ||
| 4 || 2998744.8 || 9763358.5 || | | 4 || 2998744.8 || 9763358.5 || −5587920.9 || −5239483.7 || 12764744.5 | ||
|- align="right" | |- align="right" | ||
| ∑ = || 12724613.2 || 36444385.8 || 2841056.2 || 4183151.8 || 49179425.3 | | ∑ = || 12724613.2 || 36444385.8 || 2841056.2 || 4183151.8 || 49179425.3 | ||
Строка 185: | Строка 185: | ||
! ''a''<sub>1</sub> !! ''b''<sub>1</sub> | ! ''a''<sub>1</sub> !! ''b''<sub>1</sub> | ||
|- align="right" | |- align="right" | ||
| 0.99978799 || | | 0.99978799 || −0.02728978 | ||
|} | |} | ||
{| class="wikitable" style="float:right" | {| class="wikitable" style="float:right" | ||
Строка 203: | Строка 203: | ||
| 1 || 0.002 || 0.001 | | 1 || 0.002 || 0.001 | ||
|- align="right" | |- align="right" | ||
| 2 || 0.016 || | | 2 || 0.016 || −0.013 | ||
|- align="right" | |- align="right" | ||
| 3 || | | 3 || −0.032 || −0.016 | ||
|- align="right" | |- align="right" | ||
| 4 || 0.013 || 0.028 | | 4 || 0.013 || 0.028 | ||
Строка 216: | Строка 216: | ||
! ''m'' !! ''θ'' | ! ''m'' !! ''θ'' | ||
|- align="right" | |- align="right" | ||
| 1.00016037 || | | 1.00016037 || −1°33′48.72″ | ||
|} | |} | ||
Строка 225: | Строка 225: | ||
! ''X''<sub>0</sub> !! ''Y''<sub>0</sub> !! ''x''<sub>0</sub> !! ''y''<sub>0</sub> !! ''m'' !! ''θ'' | ! ''X''<sub>0</sub> !! ''Y''<sub>0</sub> !! ''x''<sub>0</sub> !! ''y''<sub>0</sub> !! ''m'' !! ''θ'' | ||
|- align="right" | |- align="right" | ||
| 83477.64 || 87377.60 || 1334.71 || 285.94 || 1.00016037 || | | 83477.64 || 87377.60 || 1334.71 || 285.94 || 1.00016037 || −1°33′48.72″ | ||
|} | |} | ||
== Пример программной реализации == | |||
Вот листинг простой утилиты на языке C: | |||
<syntaxhighlight lang="c"> | |||
#include <stdio.h> | |||
#include <stdlib.h> | |||
#include <math.h> | |||
int main(int argc, char *argv[]) | |||
{ | |||
char buf[1024], name[32]; | |||
double x[2], y[2], wgt; | |||
double xc[2], yc[2]; | |||
double dx[2], dy[2], dz[2]; | |||
double a[2][2], scale, rotation; | |||
double s[5]; | |||
int i; | |||
FILE *fp0, *fp1; | |||
if (argc < 4) { | |||
printf("usage: findkey input-file key-file var-file\n"); | |||
exit(EXIT_FAILURE); | |||
} | |||
if ((fp0 = fopen(argv[1], "r")) == NULL) { | |||
printf("can't open %s\n", argv[1]); | |||
exit(EXIT_FAILURE); | |||
} | |||
/* подсчитать сумму координат */ | |||
for (i = 0; i < 5; i++) | |||
s[i] = 0.; | |||
while (fgets(buf, 1024, fp0) != NULL) { | |||
sscanf(buf, "%s %lf %lf %lf %lf %lf", | |||
name, &x[0], &x[1], &y[0], &y[1], &wgt); | |||
s[0] += x[0] * wgt; | |||
s[1] += x[1] * wgt; | |||
s[2] += y[0] * wgt; | |||
s[3] += y[1] * wgt; | |||
s[4] += wgt; | |||
} | |||
rewind(fp0); | |||
/* найти центр масс */ | |||
for (i = 0; i < 2; i++) { | |||
xc[i] = s[i] / s[4]; | |||
yc[i] = s[2 + i] / s[4]; | |||
} | |||
/* подсчитать сумму произведений */ | |||
for (i = 0; i < 5; i++) | |||
s[i] = 0.; | |||
while (fgets(buf, 1024, fp0) != NULL) { | |||
sscanf(buf, "%s %lf %lf %lf %lf %lf", | |||
name, &x[0], &x[1], &y[0], &y[1], &wgt); | |||
/* вычислить разности */ | |||
dx[0] = x[0] - xc[0]; | |||
dx[1] = x[1] - xc[1]; | |||
dy[0] = y[0] - yc[0]; | |||
dy[1] = y[1] - yc[1]; | |||
/* суммировать */ | |||
s[0] += dx[0] * dy[0] * wgt; | |||
s[1] += dx[1] * dy[1] * wgt; | |||
s[2] += dx[0] * dy[1] * wgt; | |||
s[3] += dx[1] * dy[0] * wgt; | |||
s[4] += (dx[0] * dx[0] + dx[1] * dx[1]) * wgt; | |||
} | |||
rewind(fp0); | |||
/* найти первичные параметры */ | |||
a[1][0] = (s[0] + s[1]) / s[4]; | |||
a[1][1] = (s[2] - s[3]) / s[4]; | |||
a[0][0] = yc[0] - a[1][0] * xc[0] + a[1][1] * xc[1]; | |||
a[0][1] = yc[1] - a[1][1] * xc[0] - a[1][0] * xc[1]; | |||
/* найти вторичные параметры */ | |||
scale = hypot(a[1][0], a[1][1]); | |||
rotation = atan2(a[1][1], a[1][0]); | |||
/* вывести параметры в файл ключа */ | |||
if ((fp1 = fopen(argv[2], "w")) == NULL) { | |||
printf("can't create %s\n", argv[2]); | |||
exit(EXIT_FAILURE); | |||
} | |||
fprintf(fp1, "%.3f\n", a[0][0]); | |||
fprintf(fp1, "%.3f\n", a[0][1]); | |||
fprintf(fp1, "%.12f\n", a[1][0]); | |||
fprintf(fp1, "%.12f\n", a[1][1]); | |||
fprintf(fp1, "%.12f\n", scale); | |||
fprintf(fp1, "%+.10f\n", rotation * 180. / M_PI); | |||
fclose(fp1); | |||
/* вывести данные вместе с невязками */ | |||
if ((fp1 = fopen(argv[3], "w")) == NULL) { | |||
printf("can't create %s\n", argv[3]); | |||
exit(EXIT_FAILURE); | |||
} | |||
while (fgets(buf, 1024, fp0) != NULL) { | |||
sscanf(buf, "%s %lf %lf %lf %lf %lf", | |||
name, &x[0], &x[1], &y[0], &y[1], &wgt); | |||
/* "наблюдённые" dx[], dy[] */ | |||
dx[0] = x[0] - xc[0]; | |||
dx[1] = x[1] - xc[1]; | |||
dy[0] = y[0] - yc[0]; | |||
dy[1] = y[1] - yc[1]; | |||
/* "вычисленные" dy[] */ | |||
dz[0] = a[1][0] * dx[0] - a[1][1] * dx[1]; | |||
dz[1] = a[1][1] * dx[0] + a[1][0] * dx[1]; | |||
fprintf(fp1, "%s %.3f %.3f %.3f %.3f %g %.3f %.3f\n", | |||
name, x[0], x[1], y[0], y[1], wgt, dy[0] - dz[0], dy[1] - dz[1]); | |||
} | |||
fclose(fp1); | |||
fclose(fp0); | |||
return 0; | |||
} | |||
</syntaxhighlight> | |||
Сохраним код в файл '''findkey.c'''. Исполняемый модуль можно создать, например, компилятором '''gcc''': | |||
<syntaxhighlight lang="bash"> | |||
$ gcc -o findkey findkey.c -lm | |||
</syntaxhighlight> | |||
Для MS Windows можно загрузить уже скомпилированную [http://gis-lab.info/other/findkey.zip программу]. | |||
Утилита '''findkey''' запускается в командной строке с именами трёх файлов в качестве аргументов: входной файл данных, выходные файл параметров и файл невязок. | |||
Подготовим файл данных '''data.dat''', структура которого соответствует таблице исходных данных, т. е. содержит в колонках названия пунктов, входные координаты ''x'', ''y'', выходные координаты ''X'', ''Y'', веса: | |||
<pre> | |||
1 1334.71 285.94 83477.64 87377.60 1.0 | |||
2 563.67 -5197.34 82557.14 81916.51 1.0 | |||
3 4444.27 1153.79 86610.19 88160.39 1.0 | |||
4 -252.07 2881.90 81962.05 90016.34 1.0 | |||
</pre> | |||
После запуска программы | |||
<syntaxhighlight lang="bash"> | |||
$ findkey data.dat key.dat var.dat | |||
</syntaxhighlight> | |||
получим параметры '''key.dat''' | |||
<pre> | |||
82135.407 | |||
87128.144 | |||
0.9997879942 | |||
-0.0272897781 | |||
1.0001603698 | |||
-1.56353244 | |||
</pre> | |||
и невязки '''var.dat''' | |||
<pre> | |||
1 1334.710 285.940 83477.640 87377.600 1 0.002 0.001 | |||
2 563.670 -5197.340 82557.140 81916.510 1 0.016 -0.013 | |||
3 4444.270 1153.790 86610.190 88160.390 1 -0.032 -0.016 | |||
4 -252.070 2881.900 81962.050 90016.340 1 0.013 0.028 | |||
</pre> | |||
== Заключение == | == Заключение == | ||
Положенные в основу статьи формулы называются в учебниках геодезии «неполными» в противоположность «полным». Дело в том, что при большом удалении объекта от осевого меридиана исходной проекции Гаусса-Крюгера возникает значительный градиент масштаба изображения в направлении запад-восток. Чтобы компенсировать возникающие при этом специфические искажения конформного отображения, в «полные» выражения добавляют необходимые члены разложений формул проекций. Разумеется, такой подход несовершенен, как любые ограниченные разложения. В статье [ | Положенные в основу статьи формулы называются в учебниках геодезии «неполными» в противоположность «полным». Дело в том, что при большом удалении объекта от осевого меридиана исходной проекции Гаусса-Крюгера возникает значительный градиент масштаба изображения в направлении запад-восток. Чтобы компенсировать возникающие при этом специфические искажения конформного отображения, в «полные» выражения добавляют необходимые члены разложений формул проекций. Разумеется, такой подход несовершенен, как любые ограниченные разложения. В статье [http://gis-lab.info/qa/local-cs.html Добавление местной координатной системы в GIS] в качестве альтернативы предлагается переход от ГСК к проекции с нулевым градиентом масштаба в центре объекта или вблизи него, что делает «полные» формулы ненужными. | ||
== Примечания == | == Примечания == | ||
Строка 236: | Строка 401: | ||
<references /> | <references /> | ||
== Ссылки | == Ссылки == | ||
* [http://www.google.ru/url?q=http://communities.bentley.com/products/geospatial/desktop/m/geospatial_desktop_gallery/206218/download.aspx&sa=U&ei=YG47UdKNPOeG4gT92oD4Ag&ved=0CDgQFjAJ&usg=AFQjCNFbjYGDGtp3cTXaKerd9LuE8l9_zQ Convert Local Coordinate Systems to Standard Coordinate Systems, McCoy J., 2012] | |||
* [http://user.gs.rmit.edu.au/rod/files/publications/COTRAN_1.pdf Coordinate transformations in surveying and mapping, Deakin R.E., 2004] | |||
* [http://kartoweb.itc.nl/geometrics/coordinate%20transformations/coordtrans.html Coordinate transformations, Knippers R., 2009] | |||
* [http://gis-lab.info/qa/affine-math.html Аффинные преобразования - математика] | * [http://gis-lab.info/qa/affine-math.html Аффинные преобразования - математика] | ||
* [http://gis-lab.info/qa/local-cs.html Добавление местной координатной системы в GIS] | |||
* [http://gis-lab.info/qa/polynom.html Полиномиальные преобразования] | * [http://gis-lab.info/qa/polynom.html Полиномиальные преобразования] | ||
* [http://gis-lab.info/qa/polynom-calc-examples.html Полиномиальные преобразования - примеры реализации] | * [http://gis-lab.info/qa/polynom-calc-examples.html Полиномиальные преобразования - примеры реализации] | ||
* [http://gis-lab.info/qa/rmse.html Среднеквадратичная ошибка (RMSE)] | * [http://gis-lab.info/qa/rmse.html Среднеквадратичная ошибка (RMSE)] |
Текущая версия от 14:57, 22 марта 2020
по адресу http://gis-lab.info/qa/helmert2d.html
Конформное преобразование на плоскости широко используется в геодезии при создании местных координатных систем на небольшие территории, ограниченные размерами населённого пункта.
Введение
Следующие формулы связывают координаты точек x, y, заданные в местной системе координат (МСК), и координаты X, Y, заданные в государственной системе координат (ГСК):
где m — масштабный множитель, θ — угол разворота, X0, Y0, x0, y0 — координаты одного из геодезических пунктов в ГСК и МСК, как правило. Этот набор параметров называется «ключ».
Исходный материал для определения ключа — пары координат пунктов геодезической сети, полученные из независимого уравнивания одних и тех же измерений в МСК и в ГСК[1]. В зависимости от класса пунктам (вернее, соответствующим парам уравнений) назначаются веса p.
Алгоритм нахождения параметров
Конформное преобразование представляется следующей математической моделью:
Определению подлежат четыре параметра: a0, b0, a1, b1,.
Очевидно, конформное преобразование является частным случаем аффинного.
Шаг 1: вычисление взвешенных средних
Шаг 2: перенос осей в центр масс
Шаг 3: вычисление a1 и b1
Шаг 4: вычисление a0 и b0
Шаг 5: вычисление невязок
Невязки позволяют выявить точки, плохо укладывающиеся в полученную модель. Классическая процедура удаляет такие «отлетающие» точки, после чего вычисление параметров повторяется без них. Робастные процедуры переназначают веса уравнениям в соответствии с невязками, и повторное вычисление повторяется с полным набором точек при том, что «отлетающие» влияют на результат незначительно.
Кроме того, невязки необходимы для оценки точности измерений и результатов.
Шаг 6: вычисление ключа
Вычислим масштабный множитель и угол разворота:
Выберем j-й пункт с малыми невязками по возможности в середине массива точек. Его координаты в обеих системах Xj, Yj, xj, yj становятся параметрами X0, Y0, x0, y0.
Пример вычисления параметров
Даны координаты четырёх пунктов x, y в исходной системе, X, Y в конечной системе с весами p:
i | xi | yi | Xi | Yi | pi |
---|---|---|---|---|---|
1 | 1334.71 | 285.94 | 83477.64 | 87377.60 | 1.0 |
2 | 563.67 | −5197.34 | 82557.14 | 81916.51 | 1.0 |
3 | 4444.27 | 1153.79 | 86610.19 | 88160.39 | 1.0 |
4 | −252.07 | 2881.90 | 81962.05 | 90016.34 | 1.0 |
∑p = | 4.0 |
i | pixi | piyi | piXi | piYi |
---|---|---|---|---|
1 | 1334.71 | 285.94 | 83477.64 | 87377.60 |
2 | 563.67 | −5197.34 | 82557.14 | 81916.51 |
3 | 4444.27 | 1153.79 | 86610.19 | 88160.39 |
4 | −252.07 | 2881.90 | 81962.05 | 90016.34 |
∑ = | 6090.58 | −875.71 | 334607.02 | 347470.84 |
Взвешенные средние:
xc | yc | Xc | Yc |
---|---|---|---|
1522.645 | −218.9275 | 83651.755 | 86867.71 |
Перенос осей в центр масс:
i | ∆xi | ∆yi | ∆Xi | ∆Yi |
---|---|---|---|---|
1 | −187.935 | 504.8675 | −174.115 | 509.89 |
2 | −958.975 | −4978.4125 | −1094.615 | −4951.20 |
3 | 2921.625 | 1372.7175 | 2958.435 | 1292.68 |
4 | −1774.715 | 3100.8275 | −1689.705 | 3148.63 |
i | pi ∆xi ∆Xi | pi ∆yi ∆Yi | pi ∆xi ∆Yi | pi ∆yi ∆Xi | pi (∆xi2 + ∆yi2) |
---|---|---|---|---|---|
1 | 32722.3 | 257426.9 | −95826.2 | −87905.0 | 290210.8 |
2 | 1049708.4 | 24649116.0 | 4748077.0 | 5449445.0 | 25704224.1 |
3 | 8643437.7 | 1774484.5 | 3776726.2 | 4061095.5 | 10420246.0 |
4 | 2998744.8 | 9763358.5 | −5587920.9 | −5239483.7 | 12764744.5 |
∑ = | 12724613.2 | 36444385.8 | 2841056.2 | 4183151.8 | 49179425.3 |
S1 | S2 | S3 | S4 | S5 |
Параметры конформного преобразования:
|
Невязки:
i | vxi | vyi |
---|---|---|
1 | 0.002 | 0.001 |
2 | 0.016 | −0.013 |
3 | −0.032 | −0.016 |
4 | 0.013 | 0.028 |
Масштаб и разворот:
m | θ |
---|---|
1.00016037 | −1°33′48.72″ |
Сконструируем ключ на основе первого геодезического пункта:
X0 | Y0 | x0 | y0 | m | θ |
---|---|---|---|---|---|
83477.64 | 87377.60 | 1334.71 | 285.94 | 1.00016037 | −1°33′48.72″ |
Пример программной реализации
Вот листинг простой утилиты на языке C:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
int main(int argc, char *argv[])
{
char buf[1024], name[32];
double x[2], y[2], wgt;
double xc[2], yc[2];
double dx[2], dy[2], dz[2];
double a[2][2], scale, rotation;
double s[5];
int i;
FILE *fp0, *fp1;
if (argc < 4) {
printf("usage: findkey input-file key-file var-file\n");
exit(EXIT_FAILURE);
}
if ((fp0 = fopen(argv[1], "r")) == NULL) {
printf("can't open %s\n", argv[1]);
exit(EXIT_FAILURE);
}
/* подсчитать сумму координат */
for (i = 0; i < 5; i++)
s[i] = 0.;
while (fgets(buf, 1024, fp0) != NULL) {
sscanf(buf, "%s %lf %lf %lf %lf %lf",
name, &x[0], &x[1], &y[0], &y[1], &wgt);
s[0] += x[0] * wgt;
s[1] += x[1] * wgt;
s[2] += y[0] * wgt;
s[3] += y[1] * wgt;
s[4] += wgt;
}
rewind(fp0);
/* найти центр масс */
for (i = 0; i < 2; i++) {
xc[i] = s[i] / s[4];
yc[i] = s[2 + i] / s[4];
}
/* подсчитать сумму произведений */
for (i = 0; i < 5; i++)
s[i] = 0.;
while (fgets(buf, 1024, fp0) != NULL) {
sscanf(buf, "%s %lf %lf %lf %lf %lf",
name, &x[0], &x[1], &y[0], &y[1], &wgt);
/* вычислить разности */
dx[0] = x[0] - xc[0];
dx[1] = x[1] - xc[1];
dy[0] = y[0] - yc[0];
dy[1] = y[1] - yc[1];
/* суммировать */
s[0] += dx[0] * dy[0] * wgt;
s[1] += dx[1] * dy[1] * wgt;
s[2] += dx[0] * dy[1] * wgt;
s[3] += dx[1] * dy[0] * wgt;
s[4] += (dx[0] * dx[0] + dx[1] * dx[1]) * wgt;
}
rewind(fp0);
/* найти первичные параметры */
a[1][0] = (s[0] + s[1]) / s[4];
a[1][1] = (s[2] - s[3]) / s[4];
a[0][0] = yc[0] - a[1][0] * xc[0] + a[1][1] * xc[1];
a[0][1] = yc[1] - a[1][1] * xc[0] - a[1][0] * xc[1];
/* найти вторичные параметры */
scale = hypot(a[1][0], a[1][1]);
rotation = atan2(a[1][1], a[1][0]);
/* вывести параметры в файл ключа */
if ((fp1 = fopen(argv[2], "w")) == NULL) {
printf("can't create %s\n", argv[2]);
exit(EXIT_FAILURE);
}
fprintf(fp1, "%.3f\n", a[0][0]);
fprintf(fp1, "%.3f\n", a[0][1]);
fprintf(fp1, "%.12f\n", a[1][0]);
fprintf(fp1, "%.12f\n", a[1][1]);
fprintf(fp1, "%.12f\n", scale);
fprintf(fp1, "%+.10f\n", rotation * 180. / M_PI);
fclose(fp1);
/* вывести данные вместе с невязками */
if ((fp1 = fopen(argv[3], "w")) == NULL) {
printf("can't create %s\n", argv[3]);
exit(EXIT_FAILURE);
}
while (fgets(buf, 1024, fp0) != NULL) {
sscanf(buf, "%s %lf %lf %lf %lf %lf",
name, &x[0], &x[1], &y[0], &y[1], &wgt);
/* "наблюдённые" dx[], dy[] */
dx[0] = x[0] - xc[0];
dx[1] = x[1] - xc[1];
dy[0] = y[0] - yc[0];
dy[1] = y[1] - yc[1];
/* "вычисленные" dy[] */
dz[0] = a[1][0] * dx[0] - a[1][1] * dx[1];
dz[1] = a[1][1] * dx[0] + a[1][0] * dx[1];
fprintf(fp1, "%s %.3f %.3f %.3f %.3f %g %.3f %.3f\n",
name, x[0], x[1], y[0], y[1], wgt, dy[0] - dz[0], dy[1] - dz[1]);
}
fclose(fp1);
fclose(fp0);
return 0;
}
Сохраним код в файл findkey.c. Исполняемый модуль можно создать, например, компилятором gcc:
$ gcc -o findkey findkey.c -lm
Для MS Windows можно загрузить уже скомпилированную программу.
Утилита findkey запускается в командной строке с именами трёх файлов в качестве аргументов: входной файл данных, выходные файл параметров и файл невязок.
Подготовим файл данных data.dat, структура которого соответствует таблице исходных данных, т. е. содержит в колонках названия пунктов, входные координаты x, y, выходные координаты X, Y, веса:
1 1334.71 285.94 83477.64 87377.60 1.0 2 563.67 -5197.34 82557.14 81916.51 1.0 3 4444.27 1153.79 86610.19 88160.39 1.0 4 -252.07 2881.90 81962.05 90016.34 1.0
После запуска программы
$ findkey data.dat key.dat var.dat
получим параметры key.dat
82135.407 87128.144 0.9997879942 -0.0272897781 1.0001603698 -1.56353244
и невязки var.dat
1 1334.710 285.940 83477.640 87377.600 1 0.002 0.001 2 563.670 -5197.340 82557.140 81916.510 1 0.016 -0.013 3 4444.270 1153.790 86610.190 88160.390 1 -0.032 -0.016 4 -252.070 2881.900 81962.050 90016.340 1 0.013 0.028
Заключение
Положенные в основу статьи формулы называются в учебниках геодезии «неполными» в противоположность «полным». Дело в том, что при большом удалении объекта от осевого меридиана исходной проекции Гаусса-Крюгера возникает значительный градиент масштаба изображения в направлении запад-восток. Чтобы компенсировать возникающие при этом специфические искажения конформного отображения, в «полные» выражения добавляют необходимые члены разложений формул проекций. Разумеется, такой подход несовершенен, как любые ограниченные разложения. В статье Добавление местной координатной системы в GIS в качестве альтернативы предлагается переход от ГСК к проекции с нулевым градиентом масштаба в центре объекта или вблизи него, что делает «полные» формулы ненужными.
Примечания
- ↑ ГКИНП (ОНТА)-01-271-03 Руководство по созданию и реконструкции городских геодезических сетей с использованием спутниковых систем ГЛОНАСС/GPS.
Ссылки
- Convert Local Coordinate Systems to Standard Coordinate Systems, McCoy J., 2012
- Coordinate transformations in surveying and mapping, Deakin R.E., 2004
- Coordinate transformations, Knippers R., 2009
- Аффинные преобразования - математика
- Добавление местной координатной системы в GIS
- Полиномиальные преобразования
- Полиномиальные преобразования - примеры реализации
- Среднеквадратичная ошибка (RMSE)