用c语言坐标正反算,高斯正反算程序代码——C语言

以下是用C语言实现的坐标正反算和高斯正反算程序代码：

坐标正反算：

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

#define pi 3.14159265358979324
#define a 6378245.0
#define ee 0.00669342162296594323

static double transformLat(double x, double y)
{
    double ret = -100.0 + 2.0 * x + 3.0 * y + 0.2 * y * y + 0.1 * x * y + 0.2 * sqrt(fabs(x));
    ret += (20.0 * sin(6.0 * x * pi) + 20.0 * sin(2.0 * x * pi)) * 2.0 / 3.0;
    ret += (20.0 * sin(y * pi) + 40.0 * sin(y / 3.0 * pi)) * 2.0 / 3.0;
    ret += (160.0 * sin(y / 12.0 * pi) + 320 * sin(y * pi / 30.0)) * 2.0 / 3.0;
    return ret;
}

static double transformLon(double x, double y)
{
    double ret = 300.0 + x + 2.0 * y + 0.1 * x * x + 0.1 * x * y + 0.1 * sqrt(fabs(x));
    ret += (20.0 * sin(6.0 * x * pi) + 20.0 * sin(2.0 * x * pi)) * 2.0 / 3.0;
    ret += (20.0 * sin(x * pi) + 40.0 * sin(x / 3.0 * pi)) * 2.0 / 3.0;
    ret += (150.0 * sin(x / 12.0 * pi) + 300.0 * sin(x / 30.0 * pi)) * 2.0 / 3.0;
    return ret;
}

void wgs2gcj(double wgLat, double wgLon, double *gcjLat, double *gcjLon)
{
    if (outOfChina(wgLat, wgLon))
    {
        *gcjLat = wgLat;
        *gcjLon = wgLon;
        return;
    }
    double dLat = transformLat(wgLon - 105.0, wgLat - 35.0);
    double dLon = transformLon(wgLon - 105.0, wgLat - 35.0);
    double radLat = wgLat / 180.0 * pi;
    double magic = sin(radLat);
    magic = 1 - ee * magic * magic;
    double sqrtMagic = sqrt(magic);
    dLat = (dLat * 180.0) / ((a * (1 - ee)) / (magic * sqrtMagic) * pi);
    dLon = (dLon * 180.0) / (a / sqrtMagic * cos(radLat) * pi);
    *gcjLat = wgLat + dLat;
    *gcjLon = wgLon + dLon;
}

void gcj2wgs(double gcjLat, double gcjLon, double *wgLat, double *wgLon)
{
    if (outOfChina(gcjLat, gcjLon))
    {
        *wgLat = gcjLat;
        *wgLon = gcjLon;
        return;
    }
    double dLat = transformLat(gcjLon - 105.0, gcjLat - 35.0);
    double dLon = transformLon(gcjLon - 105.0, gcjLat - 35.0);
    double radLat = gcjLat / 180.0 * pi;
    double magic = sin(radLat);
    magic = 1 - ee * magic * magic;
    double sqrtMagic = sqrt(magic);
    dLat = (dLat * 180.0) / ((a * (1 - ee)) / (magic * sqrtMagic) * pi);
    dLon = (dLon * 180.0) / (a / sqrtMagic * cos(radLat) * pi);
    *wgLat = gcjLat - dLat;
    *wgLon = gcjLon - dLon;
}

int outOfChina(double lat, double lon)
{
    if (lon < 72.004 || lon > 137.8347)
        return 1;
    if (lat < 0.8293 || lat > 55.8271)
        return 1;
    return 0;
}

高斯正反算：

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

#define pi 3.14159265358979324
#define a 6378245.0
#define ee 0.00669342162296594323
#define m0 0.9996
#define x_pi (3.14159265358979324 * 3000.0 / 180.0)

void gauss2wgs(double gaussX, double gaussY, double *wgLat, double *wgLon)
{
    double x = gaussX / m0;
    double y = gaussY / m0;
    double z = sqrt(1 - ee * sin(y * pi / 180.0) * sin(y * pi / 180.0));
    double z1 = (1 - ee) / z;
    double d = x / (a * z1);
    double r = pi / 180.0;
    double lat = d;
    double v = 0;
    double sign = y < 0 ? -1.0 : 1.0;
    for (int i = 0; i < 5; i++)
    {
        double t = lat;
        double c = sqrt(1 - ee * sin(t * r) * sin(t * r)) * tan(t * r);
        double f = 1 + c * c;
        double m = a * (1 - ee) / pow(f, 1.5);
        double n = a / sqrt(f);
        double tmp = n / m * c;
        v = sign * d * d / 2.0;
        lat = d + (y - v + tmp * tmp * d / 3.0) / m;
    }
    *wgLat = lat;
    *wgLon = x_pi * (x / a / z1 + (3.0 * tmp - 0.0024 * tmp * tmp * tmp) * cos(y * r) + 0.002 * cos(3.0 * y * r) - 0.0003 * cos(5.0 * y * r)) / pi;
}

void wgs2gauss(double wgLat, double wgLon, double *gaussX, double *gaussY)
{
    double lat = wgLat;
    double lon = wgLon;
    double L = lon * pi / 180.0;
    double B = lat * pi / 180.0;
    double n = (a - 6378140) / (a + 6378140);
    double e2 = 1 - (a - 6378140) * (a - 6378140) / (a * a);
    double ee2 = (a * a - 6378140 * 6378140) / (6378140 * 6378140);
    double cosB = cos(B);
    double sinB = sin(B);
    double tanB = sinB / cosB;
    double N = a / sqrt(1 - e2 * sinB * sinB);
    double t = tanB * tanB;
    double t2 = t * t;
    double c = ee2 * cosB * cosB;
    double c2 = c * c;
    double a0 = 1 + 3 * t / 4 + 45 * t2 / 64 + 175 * t * t2 / 256 + 11025 * t2 * t2 / 16384;
    double a2 = -3 * t / 2 + 15 * t2 / 16 + 525 * t * t2 / 512 - 2205 * t2 * t2 / 2048;
    double a4 = 15 * t2 / 64 - 175 * t * t2 / 256 + 945 * t2 * t2 / 4096;
    double a6 = -35 * t2 / 512 + 315 * t2 * t2 / 2048;
    double a8 = 315 * t2 * t2 / 131072;
    double s = a0 * B - sinB * cosB * (a2 + c * a4 + c2 * a6 + c2 * c * a8);
    double A = L - 117 * pi / 180;
    double A2 = A * A;
    double A3 = A2 * A;
    double A4 = A2 * A2;
    double A5 = A2 * A3;
    double A6 = A2 * A4;
    double A7 = A2 * A5;
    double A8 = A2 * A6;
    double x = 6367558.4969 * B - (662.789 + 0.58 * cos(2 * B * pi / 180) + 1.25 * cos(4 * B * pi / 180) + 0.21 * cos(6 * B * pi / 180)) * 100000 * (A2 / 2 - (5 + 3 * tanB * tanB + 10 * c - 4 * c * c - 9 * ee2) * A4 / 24 + (61 + 90 * tanB * tanB + 298 * c + 45 * tanB * tanB * tanB * tanB - 252 * ee2 - 3 * c * c) * A6 / 720);
    double y = (1 + (1 - ee2) * cosB * cosB - 3 * ee2 * sinB * sinB) * N * tanB * A - N * (1 + (1 - ee2) * cosB * cosB - ee2 * sinB * sinB) * A3 / 6 + N * (5 - 18 * tanB * tanB + t2 + 14 * c - 58 * ee2) * tanB * A5 / 120;
    *gaussX = x + 500000;
    *gaussY = y + 2700000;
}

注意：以上代码仅供参考，实际应用中需要进行一定的测试和修改以保证正确性和稳定性。