UTM投影正转算法实现

大大源码 • • 其他
    在Gissky上发现了戴勤奋老师关于坐标转换的程序以及公式，程序的界面和转换精度都值得称赞，只是可惜没有提供实现的源码，最近由于工作需要正好需要实现这部分功能，通过对戴勤奋老师提供的公式以及OGP（欧洲石油勘探组织）最新的《Coordinate Conversions and Transformations including Formulas》文档进行了研究，利用程序将UTM投影正转实现，现发布出来供大家共同学习，也希望大家指出程序中的不足。
 



/****************************************************************************
 ** 函数名：LLtoUTM
 ** 输入：   a                 椭球体长半轴
 **          f                 椭球体扁率  f=(a-b)/a 其中b代表椭球体的短半轴
 **          Lat               经过UTM投影之前的纬度
 **          Long              经过UTM投影之前的经度
 **          LongOrigin        中央经度线
 **          FN                纬度起始点，北半球为0，南半球为10000000.0m
 **          UTMNorthing 经过UTM投影后的纬度方向的坐标
 **          UTMEasting  经过UTM投影后的经度方向的坐标      
 ** 功能描述：UTM投影正转
 ** 作者：   CiCi
**  单位：   中国地质大学（北京）地球科学与资源学院
 ** 创建日期：2007年3月28日
 ** 版本：1.0
 ***************************************************************************/

LLtoUTM(
double
 a,
double
 f,
const
 
double
 Lat, 
const
 
double
 Long,
          
double
 LongOrigin,
double
 FN,
double
 
&
UTMNorthing,
double
 
&
UTMEasting)


{
      //e表示WGS84第一偏心率,eSquare表示e的平方,
     double eSquare =(2*f-f*f) ;
     double k0 = 0.9996;
     double e2Square;
     double V, T, C, A, M;
 
      //确保longtitude位于-180.00----179.9之间
     double LongTemp = (Long+180)-int((Long+180)/360)*360-180; 
     double LatRad = Lat*PI/180;
     double LongRad = LongTemp*PI/180;
     double LongOriginRad;
     LongOriginRad = LongOrigin * PI/180;
     e2Square = (eSquare)/(1-eSquare);


     V = a/sqrt(1-eSquare*sin(LatRad)*sin(LatRad));
     T = tan(LatRad)*tan(LatRad);
     C = e2Square*cos(LatRad)*cos(LatRad);
     A = cos(LatRad)*(LongRad-LongOriginRad);
     M = a*((1-eSquare/4-3*eSquare*eSquare/64-5*eSquare*eSquare*eSquare/256)*LatRad
         -(3*eSquare/8+3*eSquare*eSquare/32+45*eSquare*eSquare*eSquare/1024)*sin(2*LatRad)
         +(15*eSquare*eSquare/256+45*eSquare*eSquare*eSquare/1024)*sin(4*LatRad)
         -(35*eSquare*eSquare*eSquare/3072)*sin(6*LatRad));
 
     UTMEasting = (double)(k0*V*(A+(1-T+C)*A*A*A/6
                        + (5-18*T+T*T+72*C-58*e2Square)*A*A*A*A*A/120)+ 500000.0);
     UTMNorthing = (double)(k0*(M+V*tan(LatRad)*(A*A/2+(5-T+9*C+4*C*C)*A*A*A*A/24
                        +(61-58*T+T*T+600*C-330*e2Square)*A*A*A*A*A*A/720)));
     //南半球纬度起点为10000000.0m
     UTMNorthing=UTMNorthing+FN;
}

本程序实现的公式请参考《Coordinate Conversions and Transformations including Formulas》文档第35页，如果你需要此文档，请和我联系！也欢迎大家提出程序的不足。
[转载请注明出处，谢谢！]
版权声明：本文为oglmatrix原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
THE END
< <上一篇
下一篇>>