Конформное преобразование: различия между версиями

Материал из GIS-Lab
Перейти к навигации Перейти к поиску
Нет описания правки
м (опечатка названия столбца расчетов)
 
(не показано 36 промежуточных версий 2 участников)
Строка 1: Строка 1:
{{Статья|Черновик}}
{{Статья|Опубликована|helmert2d}}


Конформное преобразование на плоскости широко используется в геодезии при создании местных координатных систем на небольшие территории, ограниченные размерами населённого пункта.
{{Аннотация|Конформное преобразование на плоскости широко используется в геодезии при создании местных координатных систем на небольшие территории, ограниченные размерами населённого пункта.}}


== Введение ==
== Введение ==
Строка 7: Строка 7:
Следующие формулы связывают координаты точек ''x'', ''y'', заданные в местной системе координат (МСК), и координаты ''X'', ''Y'', заданные в государственной системе координат (ГСК):
Следующие формулы связывают координаты точек ''x'', ''y'', заданные в местной системе координат (МСК), и координаты ''X'', ''Y'', заданные в государственной системе координат (ГСК):


<math>\begin{array}{lcl}
<math>\begin{align}
X & = & X_0 + m \big[ (x - x_0) \cos \theta - (y - y_0) \sin \theta \big] \\
X & = X_0 + m \big[ (x - x_0) \cos \theta - (y - y_0) \sin \theta \big] \\
Y & = & Y_0 + m \big[ (x - x_0) \sin \theta + (y - y_0) \cos \theta \big]
Y & = Y_0 + m \big[ (x - x_0) \sin \theta + (y - y_0) \cos \theta \big]
\end{array}</math>
\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{array}{lcl}
<math>\begin{align}
X' & = & a_0 + a_1 x - b_1 y \\
X' & = a_0 + a_1 x - b_1 y \\
Y' & = & b_0 + b_1 x + a_1 y
Y' & = b_0 + b_1 x + a_1 y
\end{array}</math>
\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{array}{lcl}
<math>\begin{align}
a_0 & = & X_c - a_1 x_c + b_1 y_c \\
a_0 & = X_c - a_1 x_c + b_1 y_c \\
b_0 & = & Y_c - b_1 x_c - a_1 y_c
b_0 & = Y_c - b_1 x_c - a_1 y_c
\end{array}</math>
\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 || -5197.34 || 82557.14 || 81916.51 || 1.0
| 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 || -252.07 || 2881.90 || 81962.05 || 90016.34 || 1.0
| 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 || -5197.34 || 82557.14 || 81916.51
| 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 || -252.07 || 2881.90 || 81962.05 || 90016.34
| 4 || −252.07 || 2881.90 || 81962.05 || 90016.34
|- align="right"
|- align="right"
| ∑ = || 6090.58 || -875.71 || 334607.02 || 347470.84
| ∑ = || 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 || -218.9275 || 83651.755 || 86867.71
| 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 || -187.935 || 504.8675 || -174.115 || 509.89
| 1 || −187.935 || 504.8675 || −174.115 || 509.89
|- align="right"
|- align="right"
| 2 || -958.975 || -4978.4125 || -1094.615 || -4951.20
| 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 || -1774.715 || 3100.8275 || -1689.705 || 3148.63
| 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> + ∆''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 || -95826.2 || -87905.0 || 290210.8
| 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 || -5587920.9 || -5239483.7 || 12764744.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.02728978
| 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 || -0.013
| 2 || 0.016 || −0.013
|- align="right"
|- align="right"
| 3 || -0.032 || -0.016
| 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°33′48.72″
| 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 || -1°33′48.72″
| 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''', структура которого соответствует таблице исходных данных, т.&nbsp;е. содержит в колонках названия пунктов, входные координаты ''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>


== Заключение ==
== Заключение ==


Положенные в основу статьи формулы называются в учебниках геодезии «неполными» в противоположность «полным». Дело в том, что при большом удалении объекта от осевого меридиана исходной проекции Гаусса-Крюгера возникает значительный градиент масштаба изображения в направлении запад-восток. Чтобы компенсировать возникающие при этом специфические искажения конформного отображения, в «полные» выражения добавляют необходимые члены разложений формул проекций. Разумеется, такой подход несовершенен, как любые ограниченные разложения. В статье [[Добавление_местной_координатной_системы_в_GIS|Добавление местной координатной системы в GIS]] в качестве альтернативы предлагается переход от ГСК к проекции с нулевым градиентом масштаба в центре объекта или вблизи него, что делает «полные» формулы ненужными.
Положенные в основу статьи формулы называются в учебниках геодезии «неполными» в противоположность «полным». Дело в том, что при большом удалении объекта от осевого меридиана исходной проекции Гаусса-Крюгера возникает значительный градиент масштаба изображения в направлении запад-восток. Чтобы компенсировать возникающие при этом специфические искажения конформного отображения, в «полные» выражения добавляют необходимые члены разложений формул проекций. Разумеется, такой подход несовершенен, как любые ограниченные разложения. В статье [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 pixiXi piyiYi pixiYi piyiXi 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

Параметры конформного преобразования:

a1 b1
0.99978799 −0.02728978
a0 b0
82135.407 87128.144

Невязки:

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 в качестве альтернативы предлагается переход от ГСК к проекции с нулевым градиентом масштаба в центре объекта или вблизи него, что делает «полные» формулы ненужными.

Примечания

  1. ГКИНП (ОНТА)-01-271-03 Руководство по созданию и реконструкции городских геодезических сетей с использованием спутниковых систем ГЛОНАСС/GPS.

Ссылки