经纬度坐标与高斯坐标的转换代码
经纬度坐标与高斯坐标的转换代码
/功能说明：将绝对高斯坐标(y,x)转换成绝对的地理坐标(wd,jd)。/ // double y; 输入参数: 高斯坐标的横坐标，以米为单位
// double x; 输入参数: 高斯坐标的纵坐标，以米为单位
// short DH; 输入参数: 带号，表示上述高斯坐标是哪个带的
// double *L; 输出参数: 指向经度坐标的指针，其中经度坐标以秒为单位// double *B; 输出参数: 指向纬度坐标的指针，其中纬度坐标以秒为单位voidGaussToGeo(double y, double x, short DH, double *L, double *B, double LP) {
double l0; // 经差
double tf; // tf = tg(Bf0),注意要将Bf转换成以弧度为单位
double nf ; // n = y * sqrt( 1 + etf ** 2) / c, 其中etf = e'**2 * cos(Bf0) ** 2 double t_l0; // l0，经差，以度为单位
double t_B0; // B0，纬度，以度为单位
double Bf0; // Bf0
double etf; // etf,其中etf = e'**2 * cos(Bf0) ** 2
double X_3 ;
double PI=3.14159265358979;
double b_e2=0.0067385254147;
doubleb_c=6399698.90178271;
X_3 = x / 1000000.00 - 3 ; // 以兆米（1000000）为单位
// 对于克拉索夫斯基椭球，计算Bf0
Bf0 = 27.11115372595 + 9.024******** * X_3 - 0.00579740442 * pow(X_3,2)
- 0.00043532572 * pow(X_3,3) + 0.00004857285 * pow(X_3,4)
+ 0.00000215727 * pow(X_3,5) - 0.00000019399 * pow(X_3,6) ;
tf = tan(Bf0*PI/180); // tf = tg(Bf),注意这里将Bf转换成以弧度为单位
etf = b_e2 * pow(cos(Bf0*PI/180),2); // etf = e'**2 * cos(Bf) ** 2
nf = y * sqrt( 1 + etf ) / b_c; // n = y * sqrt( 1 + etf ** 2) / c
// 计算纬度，注意这里计算出来的结果是以度为单位的
t_B0 = Bf0 - (1.0+etf) * tf / PI * (90.0 * pow(nf,2)
- 7.5 * (5.0 + 3 * pow(tf,2) + etf - 9 * etf * pow(tf,2)) * pow(nf,4)
+ 0.25 * (61 + 90 * pow(tf,2) + 45 * pow(tf,4)) * pow(nf,6)) ;
// 计算经差，注意这里计算出来的结果是以度为单位的
t_l0 = (180 * nf - 30 * ( 1 + 2 * pow(tf,2) + etf ) * pow(nf,3)
+ 1.5 * (5 + 28 * pow(tf,2) + 24 * pow(tf,4)) * pow(nf,5))
/ ( PI * cos(Bf0*PI/180) ) ;
l0 = (t_l0 * 3600.0); // 将经差转成秒
if (LP == -1000)
{
*L = (double)((DH * 6 - 3) * 3600.0 + l0); // 根据带号计算出以秒为单位的绝对经度，返回指针
else
{
*L = LP * 3600.0 + l0; // 根据带号计算出以秒为单位的绝对经度，返回指针
}
//----------------------------------
*B = (double)(t_B0 * 3600.0) ; // 将纬差转成秒，并返回指针
}
/ 功能说明：（1）将地理坐标(wd,jd)转换成绝对的高斯坐标(y,x)
（2）本函数支持基于六度带（或三度带）、克拉索夫斯基椭球进行转换/
/ 适用范围：本函数适用于将地球东半球中北半球（即东经0度到东经180度，北纬0度至90度）范围
内所有地理坐标到高斯坐标的转换/
/ 使用说明：调用本函数后返回的结果应在满足精度的条件下进行四舍五入/
// double jd; 输入参数: 地理坐标的经度，以秒为单位
// double wd; 输入参数: 地理坐标的纬度，以秒为单位
// short DH; 输入参数: 三度带或六度带的带号
/ 六度带(三度带)的带号是这样得到的：从东经0度到东经180度自西向东按每6度(3度)顺序编号
(编号从1开始)，这个顺序编号就称为六度带(三度带)的带号。因此，六度带的带号的范围
是1-30，
三度带的带号的范围是1-60。
如果一个点在图号为TH的图幅中，那麽该点所处的六度带的带号就可以这样得到：将该图号的
第3、4位组成的字符串先转换成数字，再减去30。例如某点在图幅06490701中，该点所在的带号就
是49-30，即19。
如果调用本函数去进行一般的从地理坐标到基于六度带高斯坐标的变换（非邻带转换），则参
数DH的选取按前一段的方法去确定。
如果调用本函数去进行基于六度带邻带转换，则参数DH的选取先按上述方法去确定，然后看是
往前一个带还是后一个带进行邻带转换再确定是加1还是减1。/
voidGeoToGauss(double jd, double wd, short DH, short DH_width, double *y, double *x, double LP)
{
double t; // t=tgB
double L; // 中央经线的经度
double l0; // 经差
double jd_hd,wd_hd; // 将jd、wd转换成以弧度为单位
double et2; // et2 = (e' ** 2) * (cosB ** 2)
double N; // N = C / sqrt(1 + et2)
double X; // 克拉索夫斯基椭球中子午弧长
double m; // m = cosB * PI/180 * l0
doubletsin,tcos; // sinB,cosB
double PI=3.14159265358979;
double b_e2=0.0067385254147;
doubleb_c=6399698.90178271;
jd_hd = jd / 3600.0 * PI / 180.0 ; // 将以秒为单位的经度转换成弧度wd_hd = wd / 3600.0 * PI / 180.0 ; // 将以秒为单位的纬度转换成弧度
// 如果不设中央经线（缺省参数: -1000），则计算中央经线，
// 否则，使用传入的中央经线，不再使用带号和带宽参数
//L = (DH - 0.5) * DH_width ; // 计算中央经线的经度
if (LP == -1000)
{
L = (DH - 0.5) * DH_width ; // 计算中央经线的经度
}
else
{
L = LP ;
}
l0 = jd / 3600.0 - L ; // 计算经差
tsin = sin(wd_hd); // 计算sinB
tcos = cos(wd_hd); // 计算cosB
// 计算克拉索夫斯基椭球中子午弧长X
X = 111134.8611 / 3600.0 * wd - (32005.7799 * tsin + 133.9238 * pow(tsin,3) + 0.6976 * pow(tsin,5) + 0.0039 * pow(tsin,7) ) * tcos;
et2 = b_e2 * pow(tcos,2) ; // et2 = (e' ** 2) * (cosB ** 2)
N =b_c / sqrt( 1 + et2 ) ; // N = C / sqrt(1 + et2)
t = tan(wd_hd); // t=tgB
m = PI/180 * l0 * tcos; // m = cosB * PI/180 * l0
*x = X + N * t * ( 0.5 * pow(m,2)
+ (5.0 - pow(t,2) + 9.0 * et2 + 4 * pow(et2,2)) * pow(m,4)/24.0
+ (61.0 - 58.0 * pow(t,2) + pow(t,4)) * pow(m,6) / 720.0 ) ;
*y = N * ( m + ( 1.0 - pow(t,2) + et2 ) * pow(m,3) / 6.0
+ ( 5.0 - 18.0 * pow(t,2) + pow(t,4) + 14.0 * et2
- 58.0 * et2 * pow(t,2) ) * pow(m,5) / 120.0 );
}
3度带与6度带
1．我国采用6度分带和3度分带：
1∶2.5万及1∶5万的地形图采用6度分带投影，即经差为6度，从零度子午线开始，自西向东每个经差6度为一投影带，全球共分60个带，用1，2，3，4，5，……表示．即东经0～6度为第一带，其中央经线的经度为东经3度，东经6～12度为第二带，其中央经线的经度为9度。
1∶1万的地形图采用3度分带，从东经1.5度的经线开始，每隔3度为一带，用1，2，3，……表示，全球共划分120个投影带，即东经1.5～4.5度为第1带，其中央经线的经度为东经3度，东经4.5～7.5度为第2带，其中央经线的经度为东经6度．我省位于东经113度－东经120度之间，跨第38、39、40共计3个带，其中东经115.5度以西为第38带，其中央经线为东经114度；东经115.5～118.5度为39带，其中央经线为东经117度；东经118.5度以东到山海关为40带，其中央经线为东经120度。
地形图上公里网横坐标前2位就是带号，例如：1∶5万地形图上的横坐标为20345486，其中20即为带号，345486为横坐标值。
2．当地中央经线经度的计算
六度带中央经线经度的计算：当地中央经线经度＝6°×当地带号－3°，例如：地形图上的横坐标为20345，其所处的六度带的中央经线经度为：6°×20－3°＝117°（适用于1∶2．5万和1∶5万地形图）。
三度带中央经线经度的计算：中央经线经度＝3°×当地带号（适用于1∶1万地形图）。
3、如何计算当地的中央子午线？
当地中央子午线决定于当地的直角坐标系统，首先确定您的直角坐标系统是3度带还是6度带投影公式推算:
6度带中央子午线计算公式：当地经度/6=N；中央子午线L=6 * N （带号）
当没有除尽，N有余数时，中央子午线L=6*N - 3
3度带中央子午线计算公式：当地经度/3=N；中央子午线L=3 X N
我国的经度范围西起73°东至135°，可分成
六度带十一个（13号带—23号带），各带中央经线依次为（75°、81°、……123°、129°、135°）；
三度带二十二个（24号带—45号带）。各带中央经线依次为（72°、75°、……132°、135°）；
六度带可用于中小比例尺（如1：250000）测图，三度带可用于大比例尺（如1：10000）测图，城建坐标多采用三度带的高斯投影
4、如何判断投影坐标是3度带坐标还是6度带坐标
如(4231898,21655933)其中21即为带号，同样所定义的东伪偏移值也需要加上带号，如21带的东伪偏移值为21500000米。假如你的工作区经度在120度至126度范围，则该坐标系为6度带坐标系，该带的中央经度为123度。
如(2949320,36353822)其中36即为带号，已知该地点位于贵阳市附近，而从地图上我们看到贵阳大概的经度是东经108度左右，因此可以36*3=108，所以该坐标系为3度带坐标系，该带的中央经度为108度。而不可能为6度带：36*6=216。