GIS坐标转换及其编程实现

长平狐 发布于 2013/12/25 17:25
阅读 481
收藏 3

最近,在做这个坐标转换的东西,涉及到大地测量学等很深奥的东西,在这里我就不讲解那些难懂的理论了,在此,我将会把代码贴出来和大家分享,其实要编写出这个代码,还真得把大地测量相关的知识弄熟,否则是无法理解代码的。好了废话少说。

源代码

/**
  * 空间大地直角坐标->大地坐标
  */
 public GeoPoint XYZ_BLH(int ellipse, Point3D point)
 {
  GeoPoint geoPoint = new GeoPoint();
  if(geoPoint == null || point == null)
  {
   System.out.println("对象为空!");
  }
  double a = 0,b = 0,e1 = 0,e2 = 0; //a为长半轴,b为短半轴,e1为第一偏心率,e2为第一偏心率的平方
  double deta = 0.0000001;
  switch(ellipse) //选择椭球体
  {
   case 0:  //WGS84椭球体
   {
    a = 6378137.0; //长半轴
    b = 6356752.3142;//短半轴
    e1 = Math.sqrt(Math.pow(a, 2)-Math.pow(b, 2))/a;
    e2 = 0.006694384999588;
    break;
   }
   case 1:  //北京54椭球
   {
    a = 6378245.0; //长半轴
    b = 6356863.0187730473;//短半轴
    e1 = Math.sqrt(Math.pow(a, 2)-Math.pow(b, 2))/a;
    e2 = 0.006693421622966;
    break;
   }
   case 2:  //西安80椭球
   {
    a = 6378140.0; //长半轴
    b = 6356755.2881575287;//短半轴
    e1 = Math.sqrt(Math.pow(a, 2)-Math.pow(b, 2))/a;
    e2 = 0.006694384999588;
    break;
   }
   default:break;
  }
  //double W = Math.sqrt(1.0-e1*e1*Math.pow(Math.sin(point.getX()), 2));
  //double N = a/W;
  double X = point.getX();
  double Y = point.getY();
  double Z = point.getZ();
  double m = Math.sqrt(Math.pow(X, 2)+Math.pow(Y, 2));
  double L = geoPoint.getL(); //大地经度
  double B = geoPoint.getB(); //大地纬度
  double H = geoPoint.getH(); //大地高
  L = Math.atan(Y/X);
  if (L < 0)
  {
   L += Math.PI;
  }
  double e2_ = e2/(1-e2);
  double c = a*Math.sqrt(1+e2_);
  double ce2 = c*e2;
  double k = 1 + e2_;
  double front = Z/m;
  double temp = front;
  int count = 0; //迭代次数
  do
  {
   front = temp;
   m = Math.sqrt(Math.pow(X, 2)+Math.pow(Y, 2));
   temp = Z/m + ce2*front/(m*Math.sqrt(k+ Math.pow(front, 2)));
  }
  while(Math.abs(front-temp)<deta&&count<100000); //是否在误差范围之内
  B = Math.atan(temp);//求纬度
  if (B<0)
  {
   B += Math.PI;
  }
  double W = Math.sqrt(1-e1*e1*Math.sin(B)*Math.sin(B));
  double N = a/W;
  //N = (a*m - c*c)/(2*b*Z);
  System.out.println("N = " + N);
  H = m/Math.cos(B)-N;//求高
  //H = Z/Math.sin(B)-N*(1.0-e1*e1);
  //H = X/(Math.cos(B)*Math.cos(L))-N;
  //转换为角度值
  L = Math.toDegrees(L);
  B = Math.toDegrees(B);
  //H = Math.toDegrees(H);
  geoPoint.setL(L);
  geoPoint.setB(B);
  geoPoint.setH(H);
  return geoPoint;
 }
 
 /**
  * 大地坐标->空间大地直角坐标
  * @param ellipse 椭球体
  * @param geoPoint 转换前的坐标
  * @return 返回空间直角坐标
  */
 public Point3D BLH_XYZ(int ellipse, GeoPoint geoPoint)
 {
  Point3D point1 = new Point3D(1,1,1);
  if(geoPoint == null || point1 == null)
  {
   System.out.println("对象为空!");
  }
  double a = 0,b,e1 = 0,e2 = 0; //a为长半轴,b为短半轴,e1为第一偏心率,e2为第一偏心率的平方
  switch(ellipse) //选择椭球体
  {
   case 0:  //WGS84椭球体
   {
    a = 6378137.0; //长半轴
    b = 6356752.3142;//短半轴
    e1 = Math.sqrt(Math.pow(a, 2)-Math.pow(b, 2))/a;
    e2 = 0.006694384999588;
    break;
   }
   case 1:  //北京54椭球
   {
    a = 6378245.0; //长半轴
    b = 6356863.0187730473;//短半轴
    e1 = Math.sqrt(Math.pow(a, 2)-Math.pow(b, 2))/a;
    e2 = 0.006693421622966;
    break;
   }
   case 2:  //西安80椭球
   {
    a = 6378140.0; //长半轴
    b = 6356755.2881575287;//短半轴
    e1 = Math.sqrt(Math.pow(a, 2)-Math.pow(b, 2))/a;
    e2 = 0.006694384999588;
    break;
   }
   default:break;
  }
  double B = geoPoint.getB()*Math.PI/180; //转换为弧度制
  double L = geoPoint.getL()*Math.PI/180; 
  double H = geoPoint.getH();
  //计算卯酉圈曲率半径
  double N = a/Math.sqrt(1.0-e1*e1*Math.pow(Math.sin(B), 2));
  //计算空间直角坐标
  double X = (N + H)*Math.cos(B)*Math.cos(L);
  double Y = (N + H)*Math.cos(B)*Math.sin(L);
  double Z = (N*(1-e1*e1)+ H)*Math.sin(B);
  point1.X = X;
  point1.Y = Y;
  point1.Z = Z;
  return point1;
 }
 
 /**
  * 高斯-克吕格坐标正算(从大地坐标到平面直角坐标)
  * @param ellipse 椭球体
  * @param zoneWide 带宽(3度或6度)
  * @param geoPoint 大地坐标
  * @param point  直角坐标
  */
 public void BL_xy(int ellipse,int zoneWide,GeoPoint geoPoint,Point point)
 {
  if(geoPoint == null || point == null)
  {
   System.out.println("对象为空!");
  }
  double a = 0,b,e1 = 0,e2 = 0; //a为长半轴,b为短半轴,e1为第一偏心率,e2为第一偏心率的平方
  switch(ellipse) //选择椭球体
  {
   case 0:  //WGS84椭球体
   {
    a = 6378137.0; //长半轴
    b = 6356752.3142;//短半轴
    e1 = Math.sqrt(Math.pow(a, 2)-Math.pow(b, 2))/a;
    e2 = Math.sqrt(Math.pow(a, 2)-Math.pow(b, 2))/b;
    //e2 = 0.006694384999588;
    break;
   }
   case 1:  //北京54椭球
   {
    a = 6378245.0; //长半轴
    b = 6356863.0187730473;//短半轴
    e1 = Math.sqrt(Math.pow(a, 2)-Math.pow(b, 2))/a;
    e2 = Math.sqrt(Math.pow(a, 2)-Math.pow(b, 2))/b;
    //e2 = 0.006693421622966;
    break;
   }
   case 2:  //西安80椭球
   {
    a = 6378140.0; //长半轴
    b = 6356755.2881575287;//短半轴
    e1 = Math.sqrt(Math.pow(a, 2)-Math.pow(b, 2))/a;
    e2 = Math.sqrt(Math.pow(a, 2)-Math.pow(b, 2))/b;
    //e2 = 0.006694384999588;
    break;
   }
   default:break;
  }
  
  double B = geoPoint.getB();
  double L = geoPoint.getL();
  double centerL = 0.0;//投影带中央经线的度数
  int n = 0;   //投影带的带号(6°带是1-60,3°带是1-120)
  double l = 0;  //该地所在经度与中央经度的差值
  if (6 == zoneWide) //如果是6度分带
  {
   n = (int)L/6;
   if (L%6 > 0)
   {
    n = n + 1;
   }
   centerL = 6*n-3;
   l = L - centerL;
  }
  if (3 == zoneWide)//如果是3度分带
  {
   n = (int)(L-1.5)/3;
   if ((L-1.5)%3>0)
   {
    n = n +1;
   }
   centerL = 3*n;
   l = L - centerL;
  }
  //转化为弧度
  B = B*Math.PI/180;
  L = L*Math.PI/180;
  l = l*Math.PI/180;
  double X;//从赤道起算的子午线弧长
  //计算子午线弧长的系数
  double A0 = 1.0+3.0/4.0*Math.pow(e1, 2)+45.0/64.0*Math.pow(e1, 4)
  +350.0/512.0*Math.pow(e1, 6)+11025.0/16384.0*Math.pow(e1, 8);
  double A2 = -3.0/4*Math.pow(e1, 2)+60.0/64*Math.pow(e1, 4)+
  525.0/512*Math.pow(e1, 6)+17640.0/16384.0*Math.pow(e1, 8)/2.0;
  double A4 = 15.0/64.0*Math.pow(e1, 4)+210.0/512.0*Math.pow(e1, 6)+
  8820.0/16384.0*Math.pow(e1, 8)/4.0;
  double A6 = -35.0/512.0*Math.pow(e1, 6)+2520.0/16384.0*Math.pow(e1, 8)/6.0;
  double A8 = 315.0/16384.0*Math.pow(e1, 8)/8.0;
  //计算子午线弧长X
  X = a*(1.0-Math.pow(e1, 2))*A0*B + A2*Math.sin(2*B) + A4*Math.sin(4*B)
   + A6*Math.sin(6*B) + A8*Math.sin(8*B);
  double t = Math.tan(B);
  double anke = e2*Math.cos(B);
  double N = a/Math.sqrt(1.0-Math.pow(e1, 2)*Math.pow(Math.sin(B), 2));
  //坐标计算
  double x = point.getX();
  double y = point.getY();
  x = X + 1.0/2.0*N*t*Math.pow(Math.cos(B), 2)*Math.pow(l, 2)+
  1.0/24.0*N*t*(5.0-Math.pow(t, 2)+9*Math.pow(anke, 2)+4*Math.pow(anke, 4))*Math.pow(Math.cos(B), 4)*Math.pow(l, 4)
  + 1.0/720*N*t*(61.0-58*Math.pow(t, 2)+Math.pow(t, 4)+270.0*Math.pow(anke, 2)-330.0*Math.pow(anke, 2)*Math.pow(t, 2))*Math.pow(Math.cos(B), 6);
  
  y = N*Math.cos(B)*l + 1.0/6*N*(1.0-Math.pow(t, 2)+Math.pow(anke, 2))*Math.pow(Math.cos(B), 3)*Math.pow(l, 3)+
  1.0/120*N*(5.0-18*t*t+Math.pow(t, 4)+14*Math.pow(anke, 2)-58*anke*anke*t*t)*
  Math.pow(Math.cos(B), 3)*Math.pow(l, 5);
  y += 500000.0;//加上500KM
  point.setX(x);
  point.setY(y);
 }
 
 /**
  * 高斯平面直角坐标转换为大地坐标
  * @param ellipse 椭球体
  * @param zoneWide 带宽
  * @param point  平面坐标
  * @param geoPoint 大地坐标
  */
 public void xy_BL(int ellipse,double centerL,Point point,GeoPoint geoPoint)
 {
  if(geoPoint == null || point == null)
  {
   System.out.println("对象为空!");
  }
  double a = 0,b,e1 = 0,e2 = 0; //a为长半轴,b为短半轴,e1为第一偏心率,e2为第二偏心率
  switch(ellipse) //选择椭球体
  {
   case 0:  //WGS84椭球体
   {
    a = 6378137.0; //长半轴
    b = 6356752.3142;//短半轴
    e1 = Math.sqrt(Math.pow(a, 2)-Math.pow(b, 2))/a;
    e2 = Math.sqrt(Math.pow(a, 2)-Math.pow(b, 2))/b;
    //e2 = 0.006694384999588;
    break;
   }
   case 1:  //北京54椭球
   {
    a = 6378245.0; //长半轴
    b = 6356863.0187730473;//短半轴
    e1 = Math.sqrt(Math.pow(a, 2)-Math.pow(b, 2))/a;
    e2 = Math.sqrt(Math.pow(a, 2)-Math.pow(b, 2))/b;
    //e2 = 0.006693421622966;
    break;
   }
   case 2:  //西安80椭球
   {
    a = 6378140.0; //长半轴
    b = 6356755.2881575287;//短半轴
    e1 = Math.sqrt(Math.pow(a, 2)-Math.pow(b, 2))/a;
    e2 = Math.sqrt(Math.pow(a, 2)-Math.pow(b, 2))/b;
    //e2 = 0.006694384999588;
    break;
   }
   default:break;
  }
  double x = point.getX();
  double y = point.getY()-500000;
  double Bf = 0.0;//底点纬度
  double B0 = 1.0+3.0/4.0*e1*e1 + 45.0/64.0*Math.pow(e1, 4)+350.0/312.0*Math.pow(e1, 6)+
     11025.0/16384.0*Math.pow(e1, 8) + 43659.0/65536.0*Math.pow(e1, 10);
  Bf = x/(a*(1.0-Math.pow(e1, 2))*B0);
  double tf = Math.tan(Bf);
  double Mf = a*(1.0-e1*e1)/Math.sqrt(Math.pow(1.0-Math.pow(e1, 2), 3));
  double Nf = (a/Math.sqrt(1.0-Math.pow(e1, 2)))/Math.sqrt(1.0+Math.pow(e2*Math.cos(Bf),2));
  double anke = e2*Math.cos(Bf);
  double B = geoPoint.getB();
  double L = geoPoint.getL();
  //开始坐标计算
  B = Bf - tf/(2*Mf*Nf*Math.cos(Bf))*Math.pow(y, 2) + tf/(24*Mf*Math.pow(Nf, 3))*
  (5.0+3.0*tf*tf+anke*anke-9.0*anke*anke*tf*tf)*Math.pow(y, 4)-
  1.0/(720.0*Math.pow(Nf, 5)*Math.cos(Bf))*(61.0+90.0*tf*tf+45*Math.pow(tf,4))*Math.pow(y, 6);
  L = y/(Nf*Math.cos(Bf)) - (1.0+2*tf*tf+anke*anke)*Math.pow(y, 3)/(6.0*Math.pow(Nf, 3)*Math.cos(Bf))
  + (5.0+28*tf*tf+24*Math.pow(tf, 4)+6*anke*anke+8*anke*anke*tf*tf)*Math.pow(y, 5)/(120.0*Math.pow(Nf, 5)*Math.cos(Bf));
  
  B = B*180/Math.PI;
  L = L*180/Math.PI;
  L = L + centerL;//中央经线+经度差
  geoPoint.setB(B);
  geoPoint.setL(L);
 }

这些代码可能有什么不合理的地方,还有就是由空间直角坐标转换为大地坐标的时候,大地高计算有很大偏差,一致找不出来是什么原因


原文链接:http://blog.csdn.net/zhouxuguang236/article/details/6700706
加载中
返回顶部
顶部