兰勃特投影是等面积投影。
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 = ;
y = ; //基准经度
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 / ( - e2);
nn = a / Math.Sqrt( - 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 + ( - tt + cc) * aa3 / 6.0 +
(5.0 - 18.0 * tt + tt * tt + * 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 = ( - Math.Sqrt( - e2)) / ( + Math.Sqrt( - e2));
e12 = e1 * e1;
e13 = e1 * e12;
e14 = e12 * e12; mm0 = a * (( - 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 + * 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;
}