【控制测量学】-高斯投影正算公式以及java代码技术頭條網

正算公式(将经纬度转化为坐标):

java代码(附有源代码和修改后的代码):

源代码:

<code>/**  
   * 由经纬度反算成高斯投影坐标

   *    

   * @param longitude  
   * @param latitude  
   * @return  
   */ 
public static double[] GaussToBLToGauss(
  \t\t\t\t\t\t\t\tdouble longitude, double latitude) 
{  
  \tint ProjNo = 0;  int ZoneWide; // //带宽
   \tdouble[] output = new double[2];  
  \tdouble longitude1, latitude1, longitude0, X0, Y0, xval, yval;  
  \tdouble a, f, e2, ee, NN, T, C, A, M, iPI;  
  \tiPI = 0.0174532925199433; // //3.1415926535898/180.0;  
  \tZoneWide = 6; // //6度带宽
   \ta = 6378245.0;  
  \tf = 1.0 / 298.3; // 54年北京坐标系参数
   \t// //a=6378140.0; f=1/298.257; //80年西安坐标系参数
   ProjNo = (int) (longitude / ZoneWide);  
   longitude0 = ProjNo * ZoneWide + ZoneWide / 2;  
   longitude0 = longitude0 * iPI;  
   longitude1 = longitude * iPI; // 经度转换为弧度
   latitude1 = latitude * iPI; // 纬度转换为弧度
   e2 = 2 * f - f * f;  
   ee = e2 * (1.0 - e2);  
   NN = a    / Math.sqrt(1.0 - e2 * Math.sin(latitude1)      * Math.sin(latitude1));  
   T = Math.tan(latitude1) * Math.tan(latitude1);  
   C = ee * Math.cos(latitude1) * Math.cos(latitude1);  
   A = (longitude1 - longitude0) * Math.cos(latitude1);  
   M = a    * ((1 - e2 / 4 - 3 * e2 * e2 / 64 - 5 * e2 * e2 * e2 / 256)      * latitude1      
               - (3 * e2 / 8 + 3 * e2 * e2 / 32 + 45 * e2 * e2 * e2        / 1024) * Math.sin(2 * latitude1)      
               + (15 * e2 * e2 / 256 + 45 * e2 * e2 * e2 / 1024)      * Math.sin(4 * latitude1) 
               - (35 * e2 * e2 * e2 / 3072)      * Math.sin(6 * latitude1));  
   // 因为是以赤道为Y轴的，与我们南北为Y轴是相反的，所以xy与高斯投影的标准xy正好相反;  
   xval = NN    * (A + (1 - T + C) * A * A * A / 6 + (5 - 18 * T + T * T + 72      * C - 58 * ee)      
                   * A * A * A * A * A / 120);  
   yval = M    + NN    * Math.tan(latitude1)    * (A * A / 2 + (5 - T + 9 * C + 4 * C * C) 
                   * A * A * A * A / 24 + (61      - 58 * T + T * T + 600 * C - 330 * ee)      
                   * A * A * A * A * A * A / 720);  
   X0 = 1000000L * (ProjNo + 1) + 500000L;  
   Y0 = 0;  xval = xval + X0;  
   yval = yval + Y0;  
   output[0] = xval;  
   output[1] = yval;  
   return output;  

}/<code>

java代码是网上找的, 原作者对该代码很有自信, 我根据与正算公式的比较, 发现了几个不同点, 对代码做了修改.

不同点1:ee不同

代码中的"ee = e2 * (1.0 - e2)",这对应了正算公式中的e'的平方, 代码中的f就是正算公式中的扁率α. 计算ee是否与e'的平方一致, e的平方=(a²-b²)/a², e'的平方=(a²-b²)/b², 过程如下:

一. e2=2*f-f*f=2*(a-b)/a-(a-b)/a * (a-b)/a=(a²-b²)/a²;// e2=公式中e的平方,正确

二. ee=e2*(1.0-e2)=(a²-b²)b²/(a² * a²);//ee!=e'的平方

不同点2: 代码与公式某常量值不同

公式中的X=.....(...+270C-330TC)...; Y=....(...+14C-58TC)....;

对应代码分别为...(...+600*C-330*ee)..; ...(...+72*C-58*ee);

修改后的代码:

<code>/**  
\t * 由经纬度反算成高斯投影坐标

   *   
   * @param longitude  
   * @param latitude  
   * @return  
   */ 
public static double[] GaussToBLToGauss(double longitude, double latitude)  

{  
   int ProjNo = 0;  int ZoneWide; // //带宽
   double[] output = new double[2]; 
   double longitude1, latitude1, longitude0, X0, Y0, xval, yval;  
   double a, f, e2, ee, NN, T, C, A, M, iPI;  
   iPI = 0.0174532925199433; // //3.1415926535898/180.0;  
   ZoneWide = 6; // //6度带宽
   a = 6378245.0;  
   f = 1.0 / 298.3; // 54年北京坐标系参数
   // //a=6378140.0; f=1/298.257; //80年西安坐标系参数
   ProjNo = (int) (longitude / ZoneWide);  
   longitude0 = ProjNo * ZoneWide + ZoneWide / 2;  
   longitude0 = longitude0 * iPI;  
   longitude1 = longitude * iPI; // 经度转换为弧度
   latitude1 = latitude * iPI; // 纬度转换为弧度
   e2 = 2 * f - f * f;  
   ee = e2 / (1.0 - e2);  
   NN = a    / Math.sqrt(1.0 - e2 * Math.sin(latitude1)      * Math.sin(latitude1));  
   T = Math.tan(latitude1) * Math.tan(latitude1);  
   C = ee * Math.cos(latitude1) * Math.cos(latitude1);  
   A = (longitude1 - longitude0) * Math.cos(latitude1);  
   M = a    * ((1 - e2 / 4 - 3 * e2 * e2 / 64 - 5 * e2 * e2 * e2 / 256)      * latitude1      
               - (3 * e2 / 8 + 3 * e2 * e2 / 32 + 45 * e2 * e2 * e2        / 1024) * Math.sin(2 * latitude1)      
               + (15 * e2 * e2 / 256 + 45 * e2 * e2 * e2 / 1024)      * Math.sin(4 * latitude1) 
               - (35 * e2 * e2 * e2 / 3072)      * Math.sin(6 * latitude1));  
  // 因为是以赤道为Y轴的，与我们南北为Y轴是相反的，所以xy与高斯投影的标准xy正好相反;  
  xval = NN    * (A + (1 - T + C) * A * A * A / 6 + (5 - 18 * T + T * T + 14      * C - 58 * ee)      
                  * A * A * A * A * A / 120);  yval = M    + NN    * Math.tan(latitude1)    
                  * (A * A / 2 + (5 - T + 9 * C + 4 * C * C) * A * A * A * A / 24 
                 + (61      - 58 * T + T * T + 270 * C - 330 * ee)      * A * A * A * A * A * A / 720);  
  X0 = 1000000L * (ProjNo + 1) + 500000L;  
  Y0 = 0;  
  xval = xval + X0;  
  yval = yval + Y0;  
  output[0] = xval;  
  output[1] = yval;  
  return output; 
}/<code>