兰勃特投影C#实现

兰勃特投影是等面积投影。

static double PI = 3.1415926;

//---------------------------------------------------------------------------
//根据经纬度得到xy值
//lambda0,phi0,lambda,phi:基准点,目标点纬度、经度，单位度。
//xi,yi:长度，米     结果单位：公里
static public PointXY LBToXY(PointLB pmtLonLat0, PointLB pmtLonLat1)
{
    PointXY result;
    double mm, mm0;
    double aa, aa2, aa3, aa4, aa5, aa6;
    double ep2, nn, tt, cc;
    double lambda0, phi0, lambda, phi;

    const double k0 = 0.9996;
    const double a = 6378137.0;                   //地球椭球长半轴,WGS-84基准面参数(单位公里)
    //b   :double = 6356.75231425;                //地球椭球短半轴,WGS-84基准面参数(单位公里)
    //f   :double = 0.00335281066399172;          //扁平率,由a,b算出  f  := 1 - (b/a);
    //e   :double = 0.0818191908334155;           //偏心率,由f算出    sqrt(e2)
    const double e2 = 0.00669437998863487;        //e*e               e2 := 2*f - f*f;
    const double e4 = 4.4814723432235E-5;         //e2*e2
    const double e6 = 3.0000678774096E-7;         //e4*e2
    //e8  :double = 2.00835943630771E-9;          //e4*e4
    //   const double PI=3.1415926535897;


    lambda0 = pmtLonLat0.lon;
    phi0 = pmtLonLat0.lat;
    lambda = pmtLonLat1.lon;
    phi = pmtLonLat1.lat;
    double x, y;
    x = 99999999;
    y = 99999999;

    //基准经度
    lambda0 = lambda0 * PI / 180.0;
    phi0 = phi0 * PI / 180.0;

    //转为弧度
    lambda = lambda * PI / 180.0;    //经度
    phi = phi * PI / 180.0;            //纬度

    mm = a * ((1.0 - e2 / 4.0 - 3.0 * e4 / 64.0 - 5.0 * e6 / 256.0) * phi -
            (3.0 * e2 / 8.0 + 3.0 * e4 / 32.0 + 45.0 * e6 / 1024.0) * Math.Sin(2.0 * phi) +
            (15.0 * e4 / 256.0 + 45.0 * e6 / 1024.0) * Math.Sin(4.0 * phi) -
            (35.0 * e6 / 3072.0) * Math.Sin(6.0 * phi));

    mm0 = a * ((1.0 - e2 / 4.0 - 3.0 * e4 / 64.0 - 5.0 * e6 / 256.0) * phi0 -
            (3.0 * e2 / 8.0 + 3.0 * e4 / 32.0 + 45.0 * e6 / 1024.0) * Math.Sin(2.0 * phi0) +
            (15.0 * e4 / 256.0 + 45.0 * e6 / 1024.0) * Math.Sin(4.0 * phi0) -
            (35.0 * e6 / 3072.0) * Math.Sin(6.0 * phi0));

    aa = (lambda - lambda0) * Math.Cos(phi);
    aa2 = aa * aa;
    aa3 = aa2 * aa;
    aa4 = aa2 * aa2;
    aa5 = aa4 * aa;
    aa6 = aa3 * aa3;

    ep2 = e2 / (1 - e2);
    nn = a / Math.Sqrt(1 - e2 * Math.Sin(phi) * Math.Sin(phi));
    tt = Math.Tan(phi) * Math.Tan(phi);
    cc = ep2 * Math.Cos(phi) * Math.Cos(phi);

    x = k0 * nn * (aa + (1 - tt + cc) * aa3 / 6.0 +
            (5.0 - 18.0 * tt + tt * tt + 72 * cc - 58.0 * ep2) * aa5 / 120.0);
    y = k0 * (mm - mm0 + nn * Math.Tan(phi) * (aa2 / 2.0 + (5.0 - tt + 9.0 * cc + 4.0 * cc * cc) * aa4 / 24.0 +
            (61.0 - 58.0 * tt + tt * tt + 600.0 * cc - 330.0 * ep2) * aa6 / 720.0));
    y = -y;

    result.x = x * 1.0 / 1000.0;
    result.y = y * 1.0 / 1000.0;

    return result;
}

//根据xy值获得经纬度,采用UTM投影,WGS-84(GPS standard)地球数据
//lambda0,phi0:基准点纬度、经度，单位度，
//lambda:经度  phi:纬度单位度
//p1为离p0点的距离（公里）
static public PointLB XYtoLB(PointLB p0, PointXY p1)
{
    PointLB result;
    double phi1;
    double x, y, e1, e12, e13, e14;
    double mm, mm0, mu, ep2, cc1, tt1, nn1, rr1;
    double dd, dd2, dd3, dd4, dd5, dd6;

    double phi, lambda, lambda0, phi0;

    const double k0 = 0.9996;
    const double a = 6378137.0;                //地球椭球长半轴,WGS-84基准面参数(单位公里)
    const double e2 = 0.00669437998863487;    //e*e               e2 := 2*f - f*f;
    const double e4 = 4.4814723432235E-5;    //e2*e2
    const double e6 = 3.0000678774096E-7;    //e4*e2
    const double pi = 3.1415926535897;

    phi = -1.0;
    lambda = -1.0;
    lambda0 = p0.lon;
    phi0 = p0.lat;

    lambda0 = lambda0 * pi / 180.0;    //基准点经度,由度转为弧度
    phi0 = phi0 * pi / 180.0;        //基准点纬度，当phi0为零时即为UTM投影，以赤道做为基准点

    x = p1.x * 1000.0;    //得到离基准点真实x值
    y = -p1.y * 1000.0;    //得到离基准点真实y值

    e1 = (1 - Math.Sqrt(1 - e2)) / (1 + Math.Sqrt(1 - e2));
    e12 = e1 * e1;
    e13 = e1 * e12;
    e14 = e12 * e12;

    mm0 = a * ((1 - e2 / 4.0 - 3.0 * e4 / 64.0 - 5.0 * e6 / 256.0) * phi0 -
            (3.0 * e2 / 8.0 + 3.0 * e4 / 32.0 + 45.0 * e6 / 1024.0) * Math.Sin(2.0 * phi0) +
            (15.0 * e4 / 256.0 + 45.0 * e6 / 1024.0) * Math.Sin(4.0 * phi0) -
            (35.0 * e6 / 3072.0) * Math.Sin(6.0 * phi0));
    mm = mm0 + y / k0;
    mu = mm / (a * (1.0 - e2 / 4.0 - 3.0 * e4 / 64.0 - 5.0 * e6 / 256.0));

    phi1 = mu + (3.0 * e1 / 2.0 - 27.0 * e13 / 32.0) * Math.Sin(2.0 * mu) +
            (21.0 * e12 / 16.0 - 55.0 * e14 / 32.0) * Math.Sin(4.0 * mu) +
            (151.0 * e13 / 96.0) * Math.Sin(6.0 * mu) +
            (1097.0 * e14 / 512.0) * Math.Sin(8.0 * mu);

    // 计算lambda 和 phi
    ep2 = e2 / (1.0 - e2);
    cc1 = ep2 * Math.Cos(phi1) * Math.Cos(phi1);
    tt1 = Math.Tan(phi1) * Math.Tan(phi1);
    nn1 = a / Math.Sqrt(1.0 - e2 * Math.Sin(phi1) * Math.Sin(phi1));
    rr1 = a * (1.0 - e2) / Math.Pow(1.0 - e2 * Math.Sin(phi1) * Math.Sin(phi1), 1.5);
    dd = x / (nn1 * k0);

    dd2 = dd * dd;
    dd3 = dd * dd2;
    dd4 = dd2 * dd2;
    dd5 = dd3 * dd2;
    dd6 = dd4 * dd2;

    lambda = lambda0 + (dd - (1.0 + 2.0 * tt1 + cc1) * dd3 / 6.0 +
            (5.0 - 2.0 * cc1 + 28.0 * tt1 - 3.0 * cc1 * cc1 + 8.0 * ep2 + 24 * tt1 * tt1) * dd5 / 120.0) / Math.Cos(phi1);
    phi = phi1 - (nn1 * Math.Tan(phi1) / rr1) *
            (dd2 / 2.0 - (5.0 + 3.0 * tt1 + 10.0 * cc1 - 4.0 * cc1 * cc1 - 9.0 * ep2) * dd4 / 24.0 +
            (61.0 + 90.0 * tt1 + 298.0 * cc1 + 45.0 * tt1 * tt1 - 252.0 * ep2 - 3.0 * cc1 * cc1) * dd6 / 720.0);

    //转为度
    result.lon = lambda * 180.0 / pi;    //经度
    result.lat = phi * 180.0 / pi;        //纬度

    return result;
}