Страница 1 из 1

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

Добавлено: 22 фев 2019, 13:18
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]]
Откуда может быть такое приличное смещение от истинного расположения?

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

Добавлено: 22 фев 2019, 13:40
Игорь Белов
  1. Создайте набор точек в проекции без аффинного преобразования

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

    "Выборг-0", 8, 1001, 7, 27.95, 0, 1, 250000, -11057.626
    Нет, даже не так. Просто создайте произвольный набор точек на плоскости.
  2. Получите второй набор точек, использовав аффинное преобразование.
  3. Убедитесь, что обратное аффинное преобразование второго набора воспроизводит первый.

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

Добавлено: 22 фев 2019, 15:35
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

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

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

Добавлено: 22 фев 2019, 16:37
gamm
Zorgis писал(а):
22 фев 2019, 15:35
А где можно посмотреть исходный алгоритм обратного аффинного преобразования?
а чего там смотреть, там система двух линейных уравнений с двумя неизвестными из школьной программы. Единственно, для численной устойчивости, данные нужно центровать (вычесть среднее, и поделить на размах). А-матрица 2х2, x,y,b - вектора. Обратная матрица inv(A) есть в вики.
прямое - y=Ax+b
обратное - x=Inv(A)(y-b)

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

Добавлено: 22 фев 2019, 19:37
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;
        }