经纬度坐标与高斯坐标的转换代码
合集下载
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
double l0; //经差
double jd_hd,wd_hd; //将jd、wd转换成以弧度为单位
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
/适用范围：本函数适用于将地球东半球中北半球（即东经0度到东经180度，北纬0度至90度）范围
内所有地理坐标到高斯坐标的转换/
/使用说明：调用本函数后返回的结果应在满足精度的条件下进行四舍五入/
// double jd;输入参数:地理坐标的经度，以秒为单位
// double wd;输入参数:地理坐标的纬度，以秒为单位
第3、4位组成的字符串先转换成数字，再减去30。例如某点在图幅06490701中，该点所在的带号就
是49-30，即19。
如果调用本函数去进行一般的从地理坐标到基于六度带高斯坐标的变换（非邻带转换），则参
数DH的选取按前一段的方法去确定。
如果调用本函数去进行基于六度带邻带转换，则参数DH的选取先按上述方法去确定，然后看是
// short DH;输入参数:三度带或六度带的带号
/六度带(三度带)的带号是这样得到的：从东经0度到东经180度自西向东按每6度(3度)顺序编号
(编号从1开始)，这个顺序编号就称为六度带(三度带)的带号。因此，六度带的带号的范围是1-30，
三度带的带号的范围是1-60。
如果一个点在图号为TH的图幅中，那麽该点所处的六度带的带号就可以这样得到：将该图号的
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)) ;
}
//----------------------------------
*B = (double)(t_B0 * 3600.0) ; //将纬差转成秒，并返回指针
}
/功能说明：（1）将地理坐标(wd,jd)转换成绝对的高斯坐标(y,x)
（2）本函数支持基于六度带（或三度带）、克拉索夫斯基椭球进行转换/
Bf0 = 27.11115372595 + 9.02468257083 * 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) ;
//计算经差，注意这里计算出来的结果是以度为单位的
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) ) ;
往前一个带还是后一个带进行邻带转换再确定是加1还是减1。/
void GeoToGauss(double jd, double wd, short DH, short DH_width, double *y, double *x, double LP)
{
double t; // t=tgB
double L; //中央经线的经度
// double *B;输出参数:指向纬度坐标的指针，其中纬度坐标以秒为单位
void GaussToGeo(double y, double x, short DH, double *L, double *B, double LP)
{
double l0; //经差
double tf; // tf = tg(Bf0),注意要将Bf转换成以弧度为单位
l0 = (t_l0 * 3600.0); //将经差转成秒
if (LP == -1000)
{
*L = (double)((DH * 6 - 3) * 3600.0 + l0); //根据带号计算出以秒为单位的绝对经度，返回ቤተ መጻሕፍቲ ባይዱ针
}
else
{
*L = LP * 3600.0 + l0; //根据带号计算出以秒为单位的绝对经度，返回指针
经纬度坐标与高斯坐标的转换代码
/功能说明：将绝对高斯坐标(y,x)转换成绝对的地理坐标(wd,jd)。/
// double y;输入参数:高斯坐标的横坐标，以米为单位
// double x;输入参数:高斯坐标的纵坐标，以米为单位
// short DH;输入参数:带号，表示上述高斯坐标是哪个带的
// double *L;输出参数:指向经度坐标的指针，其中经度坐标以秒为单位
double X_3 ;
double PI=3.14159265358979;
double b_e2=0.0067385254147;
double b_c=6399698.90178271;
X_3 = x / 1000000.00 - 3 ; //以兆米（1000000）为单位
//对于克拉索夫斯基椭球，计算Bf0
double jd_hd,wd_hd; //将jd、wd转换成以弧度为单位
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
/适用范围：本函数适用于将地球东半球中北半球（即东经0度到东经180度，北纬0度至90度）范围
内所有地理坐标到高斯坐标的转换/
/使用说明：调用本函数后返回的结果应在满足精度的条件下进行四舍五入/
// double jd;输入参数:地理坐标的经度，以秒为单位
// double wd;输入参数:地理坐标的纬度，以秒为单位
第3、4位组成的字符串先转换成数字，再减去30。例如某点在图幅06490701中，该点所在的带号就
是49-30，即19。
如果调用本函数去进行一般的从地理坐标到基于六度带高斯坐标的变换（非邻带转换），则参
数DH的选取按前一段的方法去确定。
如果调用本函数去进行基于六度带邻带转换，则参数DH的选取先按上述方法去确定，然后看是
// short DH;输入参数:三度带或六度带的带号
/六度带(三度带)的带号是这样得到的：从东经0度到东经180度自西向东按每6度(3度)顺序编号
(编号从1开始)，这个顺序编号就称为六度带(三度带)的带号。因此，六度带的带号的范围是1-30，
三度带的带号的范围是1-60。
如果一个点在图号为TH的图幅中，那麽该点所处的六度带的带号就可以这样得到：将该图号的
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)) ;
}
//----------------------------------
*B = (double)(t_B0 * 3600.0) ; //将纬差转成秒，并返回指针
}
/功能说明：（1）将地理坐标(wd,jd)转换成绝对的高斯坐标(y,x)
（2）本函数支持基于六度带（或三度带）、克拉索夫斯基椭球进行转换/
Bf0 = 27.11115372595 + 9.02468257083 * 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) ;
//计算经差，注意这里计算出来的结果是以度为单位的
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) ) ;
往前一个带还是后一个带进行邻带转换再确定是加1还是减1。/
void GeoToGauss(double jd, double wd, short DH, short DH_width, double *y, double *x, double LP)
{
double t; // t=tgB
double L; //中央经线的经度
// double *B;输出参数:指向纬度坐标的指针，其中纬度坐标以秒为单位
void GaussToGeo(double y, double x, short DH, double *L, double *B, double LP)
{
double l0; //经差
double tf; // tf = tg(Bf0),注意要将Bf转换成以弧度为单位
l0 = (t_l0 * 3600.0); //将经差转成秒
if (LP == -1000)
{
*L = (double)((DH * 6 - 3) * 3600.0 + l0); //根据带号计算出以秒为单位的绝对经度，返回ቤተ መጻሕፍቲ ባይዱ针
}
else
{
*L = LP * 3600.0 + l0; //根据带号计算出以秒为单位的绝对经度，返回指针
经纬度坐标与高斯坐标的转换代码
/功能说明：将绝对高斯坐标(y,x)转换成绝对的地理坐标(wd,jd)。/
// double y;输入参数:高斯坐标的横坐标，以米为单位
// double x;输入参数:高斯坐标的纵坐标，以米为单位
// short DH;输入参数:带号，表示上述高斯坐标是哪个带的
// double *L;输出参数:指向经度坐标的指针，其中经度坐标以秒为单位
double X_3 ;
double PI=3.14159265358979;
double b_e2=0.0067385254147;
double b_c=6399698.90178271;
X_3 = x / 1000000.00 - 3 ; //以兆米（1000000）为单位
//对于克拉索夫斯基椭球，计算Bf0