Преобразование координат

October 28, 2009

Разбираясь в своих архивах, нашел много всего почти позабытого. В том числе утилитки для преобразования координат из WGS84 в Пулково 1942 и "морских" координат. Разумеется, эти преобразования можно выполнить с помощью универсальных библиотек, но они огромны и не всегда эффективны, не говоря уж про зависимости. Так что любителям минимализма и встраиваемых решений посвящается.

#include <stdio.h>
#include <math.h>
#include <stdlib.h>

double pifact = 0.017453292519943295; // PI/180 - перевод градусов в радианы

// параметры референц-эллипсоида Красовского
double aval = 6378245.0;    // полуось эллипсоида
double f = 298.3;           // величина сжатия
// параметры системы координат 42 года (СК Пулково 42)
double cmlon = 45;          // центральный меридиан 8-й зоны
double orglat = 0;          // центральная параллель (экватор)
double scale = 1;           // масштабный коэффициент
double fe = 500000;         // восточное смещение (без номера зоны)
double fn = 0;              // северное смещение (в нашем полушарии 0)
int latNS = 1;              // полушарие (1 - северное, -1 - южное)
int lonEW = 1;              // смещение от Гринвича (1 - E, -1 - W)

int main(int argc, char *argv[])
{
  // пересчет координат GPS из геодезических в картографические
  double g = 1/f;
  double esq = 2*g-g*g;
  double epsq = esq/(1-esq);
  // координаты точки в десятичных градусах
  double inlat = 56.33847;    // широта
  double inlon = 43.98538;    // долгота

  printf("latitude = %g\n",inlat);
  printf("longitude = %g\n",inlon);

  inlat = inlat * pifact * latNS;
  orglat = orglat * pifact * latNS;
  inlon = inlon * pifact * lonEW;
  cmlon = cmlon * pifact * lonEW;
  double t = tan(inlat) * tan(inlat);//T
  //printf("T = %g\n",t);
  double c = epsq * cos(inlat) * cos(inlat);//C
  //printf("C = %g\n",c);
  double a = (inlon - cmlon) * cos(inlat);//A
  //printf("A = %g\n",a);
  double mint;
  mint = inlat * (1.-esq/4.-3*esq*esq/64.-5*esq*esq*esq/256);
  mint = mint - sin(2.*inlat)*(3*esq/8+3*esq*esq/32+45*esq*esq*esq/1024);
  mint = mint + sin(4.*inlat)*esq*esq*(15/256 + 45*esq/1024);
  double m = 1.*aval*(mint - sin(6.*inlat)*35*esq*esq*esq/3072);//M
  //printf("M = %g\n",m);
  mint = orglat * (1.-esq/4.-3*esq*esq/64.-5*esq*esq*esq/256);
  mint = mint - sin(2.*orglat)*(3*esq/8+3*esq*esq/32+45*esq*esq*esq/1024);
  mint = mint + sin(4.*orglat)*esq*esq*(15/256 + 45*esq/1024);
  double m0 = 1.*aval*(mint - sin(6.*orglat)*35*esq*esq*esq/3072);//M0
  //printf("M0 = %g\n",m0);
  double nu = 1.*aval/sqrt(1.-esq*sin(inlat)*sin(inlat));//NU
  //printf("NU = %g\n",nu);
  double easting = (5-18*t+t*t+72*c-58*epsq)*a*a/120 + (1-t+c)/6;
  easting = (1 +easting*a*a)*a*scale*nu + fe;
  double northing = a*a*(61-58*t+t*t+600*c-330*epsq)/720 + (5-t+9*c+4*c*c)/24;
  northing = (.5 + a*a*northing)*a*a*nu*tan(inlat)+m-m0;
  northing = fn + scale * northing;

  // координаты точки в СК Пулково 42, метры
  // в соответствии с правилами геодезии, приняты следующие обозначения
  // X - отсчитывается вверх от экватора (ось ординат)
  // Y - отсчитывается вправо от начального меридиана зоны (ось абсцисс)
  printf("X = %g\n",northing);
  printf("Y = %g\n",easting);

  return 0;        // завершение программы
}

#include <stdio.h>
#include <math.h>
#include <stdlib.h>

int main(int argc, char *argv[])
{
 double xpos = 4784270;
 double ypos = 7803506;

 double RadToDeg = 57.2957795132;
 double DegToRad = 0.0174532925199;
 double b = 6356752.3142;
 double PI = 3.141592654;
 double HALF_PI = 1.570796327;

 double MerToGeoLong = xpos * RadToDeg / b;
 printf("GeoLong = %f\n",MerToGeoLong);

 double MerToGeoLat = RadToDeg * (2 * atan(exp(ypos / b)) - HALF_PI);
 printf("GeoLat = %f\n",MerToGeoLat);
}