Преобразование координат
Разбираясь в своих архивах, нашел много всего почти позабытого. В том числе утилитки для преобразования координат из 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); }
Comments
gdal - для растров
ogr - для векторов
Оба пакета понимают уйму самых различных форматов ГИС-данных.