Перевод СК Выборга формата Mapinfow.prj в WKT плюс аффинное преобразование

[Zorgis]

Добрый день,
Есть МСК "Выборг" в формате Mapinfow.prj.
"Выборг", 1008, 1001, 7, 27.95, 0, 1, 250000, -11057.626, 7, 0.999477237, 0.031204196, -485909.751, -0.031204196, 0.999477237, -6699985.630
Утилитой TransCoor геометрии трансформируются корректно.
Нужно программно преобразовать с помощью обратного аффинного преобразования и трансформации СК "Выборг" в формате WKT -> WGS84.
Но программно получается смещение в районе 3-4 км на север. В восточном все нормально.
Функция реафинного преобразования:

Код: Выделить всё

Параметры:  0.999477237, 0.031204196, -485909.751, -0.031204196,  0.999477237,   -6699985.630
public Double[] reAffine(Double[] points, Double[] koeff) {
    double k, kx, ky;
    k  = koeff[0] * koeff[4] - koeff[1] * koeff[3];
    kx = koeff[1] * koeff[5] - koeff[2] * koeff[4];
    ky = koeff[0] * koeff[5] - koeff[2] * koeff[3];

    for (int i = 0; i < points.length; i += 2)
    {
        points[i] = (koeff[4] * points[i] - koeff[1] * points[i + 1] + kx) / k;
        points[i + 1] = -1 * (koeff[3] * points[i] - koeff[0] * points[i + 1] + ky) / k;
    }
    return points;
}

Раскладывая параметры из MapInfo в WKT:
"Выборг", 1008, 1001, 7, 27.95, 0, 1, 250000, -11057.626, 7, 0.999477237, 0.031204196, -485909.751, -0.031204196, 0.999477237, -6699985.630
1001 - датум Pulkovo 1942
Параметры датума:
Pulkovo 1942 : имя датума; KRASSOVSKY: эллипсоид Крассовского; TOWGS84: 24, -123, -94, -0.02, 0.25, 0.13, 1.1
Ед.измерения: метры
5 параметров:
27.95, 0, 1, 250000, -11057.626
Аффинные параметры в WKT не используются, т.к. считаются отдельно.

Код: Выделить всё

WKT:PROJCS["Выборг",
  GEOGCS["MI Custom",
    DATUM["Pulkovo_1942",
      SPHEROID["Krassowsky 1940", 6378245.0, 298.3, AUTHORITY["EPSG","7024"]],
      TOWGS84[24.0, -123.0, -94.0, -0.02, 0.25, 0.13, 1.1]],
    PRIMEM["Greenwich", 0.0],
    UNIT["degree", 0.017453292519943295],
    AXIS["Longitude", EAST],
    AXIS["Latitude", NORTH]],
  PROJECTION["Transverse_Mercator"],
  PARAMETER["central_meridian", 27.95],
  PARAMETER["latitude_of_origin", 0.0],
  PARAMETER["scale_factor", 1.0],
  PARAMETER["false_easting", 250000.0],
  PARAMETER["false_northing", -11057.626],
  UNIT["m", 1.0],
  AXIS["x", EAST],
  AXIS["y", NORTH]]

Откуда может быть такое приличное смещение от истинного расположения?

[Игорь Белов]

1. Создайте набор точек в проекции без аффинного преобразования

   Код: Выделить всё

   "Выборг-0", 8, 1001, 7, 27.95, 0, 1, 250000, -11057.626

   Нет, даже не так. Просто создайте произвольный набор точек на плоскости.
2. Получите второй набор точек, использовав аффинное преобразование.
3. Убедитесь, что обратное аффинное преобразование второго набора воспроизводит первый.

[Zorgis]

да, интересно
на входе:
1000.15, 1000.15,
1010.15, 1000.15,
1010.15, 1010.15,
1000.15, 1010.15;
на выходе:
1000.1499999999775, 16162.700469148733,
1010.1500000000362,16162.700632283417,
1010.149999999978, 16172.690894568761,
1000.1499999999775,16172.690731434077

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

[gamm]

А где можно посмотреть исходный алгоритм обратного аффинного преобразования?

а чего там смотреть, там система двух линейных уравнений с двумя неизвестными из школьной программы. Единственно, для численной устойчивости, данные нужно центровать (вычесть среднее, и поделить на размах). А-матрица 2х2, x,y,b - вектора. Обратная матрица inv(A) есть в вики.
прямое - y=Ax+b
обратное - x=Inv(A)(y-b)

[trir]

Код: Выделить всё

        public double[] ReAffine(double[] points, double[] koeff)
        {
            double k, kx, ky, nx, ny;
            k  = koeff[0] * koeff[4] - koeff[1] * koeff[3];
            kx = koeff[1] * koeff[5] - koeff[2] * koeff[4];
            ky = koeff[0] * koeff[5] - koeff[2] * koeff[3];

            for (int i = 0; i < points.Length; i += 2)
            {
                nx = (koeff[4] * points[i] - koeff[1] * points[i + 1] + kx) / k;
                ny = -1 * (koeff[3] * points[i] - koeff[0] * points[i + 1] + ky) / k;
                points[i] = nx;
                points[i + 1] = ny;
            }
            return points;
        }