高斯投影坐标转换（c语言）
- 一、高斯投影简介
- 1.1 高斯投影基本概念
- 1.2 高斯投影的条件
- 1.3 高斯投影分带
- 二、高斯投影转换
- 2.1 几个概念
- 2.2 计算已知量
- 2.3 计算步骤
- 三、具体问题求解
- 3.1 问题重述
- 3.2 实现代码
- 3.3 输出结果
一、高斯投影简介
1.1 高斯投影基本概念
高斯投影又称横轴椭圆柱等角投影，属于正形投影。世界上最先采用高斯投影的国家是奥地利和德国，我国于1952年正式决定采用高斯投影。
1.2 高斯投影的条件
（1）投影后角度不产生变形，满足正形投影要求；
（2）中央子午线投影后是一条直线；
（3）中央子午线投影后长度不变，其投影长度比恒等于1。
高斯投影除了在中央子午线上没有长度变形外，不在中央子午线上的各点，其长度比都大于1，且离开中央子午线愈远，长度变形愈大。
1.3 高斯投影分带
为限制长度投影变形，投影分带有6度分带和3度分带两种方法。
- 6°带带号N和中央子午线经度 LN的关系式： LN=6N-3
- 3°带带号n和中央子午线经度 Ln的关系式： Ln=3n
- 6°带与3°带带号之间的关系为：n=2N-1
二、高斯投影转换
2.1 几个概念
高斯平面坐标系：指的是以中央子午线与赤道的交点作为坐标原点，以中央子午线的投影为纵坐标轴X，规定X轴向北为正，以赤道的投影为横坐标轴Y，Y轴向东为正，形成的坐标系。
通用坐标：在同一投影带内横坐标有正值、有负值，这对坐标的计算和使用不方便。为了使Y值都为正，将纵坐标x轴西移500km，并在Y坐标前面冠以带号，称为通用坐标。
大地坐标（Geodetic coordinate）：是大地测量中以参考椭球面为基准面的坐标，地面点P的位置用大地经度L、大地纬度B和大地高H表示。
高斯坐标正算：大地坐标（B、L、H）转换为高斯投影坐标（X、Y、Z）
高斯坐标反算：高斯投影坐标（X、Y、Z）转换为大地坐标（B、L、H）
注意：本文高斯坐标投影转换指的是高斯投影的 6° 带转换成 3° 带 或者 3° 带转换成 6° 带。
2.2 计算已知量
地球属性 | 符号表示 | 值 |
长半轴 | a | 6378245 |
短半轴 | b | 6356863.01877 |
c | 6399698.9018 | |
第一偏心率 | e1^2 | 0.0066934216 |
第二偏心率 | e2^2 | 0.0067385254 |
2.3 计算步骤
以6° 带转换成 3° 带为例， 3° 带转换成 6° 带同。
- 通用值转换成自然值
- 高斯投影坐标反算
- 求解6°带中央经线
- 求解3°带的经差
- 高斯投影坐标正算
- 自然值转换成通用值
三、具体问题求解
3.1 问题重述
已知：某点P在1954年北京坐标系6°带的平面直角坐标为
x1 = 3589644.286m
y1 = 20679136.438m
求：P点在3°带第40带的平面坐标（x2，y2）
3.2 实现代码
/*-------------------------------------------
高斯坐标代号转换代码：6°带转换成3°带
时间：2019年12月1日
编译：鸽子
编程语言：C语言
---------------------------------------------*/
#include <stdio.h>
#include <math.h>
/* 定义已知值 */
double a = 6378245; //长半轴
double b = 6356863.01877; //短半轴
double c = 6399698.9018; //c
double e_square = 0.0066934216; //第一偏心率
double e1_square = 0.0067385254; //第二偏心率
double e = sqrt(e_square);
double e1 = sqrt(e1_square);
double beita0 = 1 - 3.0 / 4.0*e1_square + 45 / 65 * pow(e1_square, 2) - 175 / 256 * pow(e1_square, 3) + 11025 / 16384 * pow(e1_square, 4); //β0，下同
double beita2 = beita0 - 1;
double beita4 = 15 / 32 * pow(e1_square, 2) - 175 / 384 * pow(e1_square, 3) + 3675 / 8192 * pow(e1_square, 4);
double beita6 = -35 / 96 * pow(e1_square, 3) + 735 / 2048 * pow(e1_square, 4);
double beita8 = 315 / 1024 * pow(e1_square, 4);
/*------------------------------高斯坐标反算-------------------------------------*/
double Fx(double B)
{
return (c*beita2 + (c*beita4 + (c*beita6 + c*beita8*pow(cos(B), 2))*pow(cos(B), 2))*pow(cos(B), 2))*sin(B)*cos(B);
}
double Fx(double B, double l)
{
double N = a / sqrt(1 - e_square*pow(sin(B), 2)); //卯酉圈曲率半径
double t = tan(B);
double fai = e1*cos(B);
double a2 = 1.0 / 2.0 * N*sin(B)*cos(B);
double a4 = 1.0 / 24.0 * N*sin(B)*pow(cos(B), 3)*(5 - pow(t, 2) + 9 * pow(fai, 2) + 4 * pow(fai, 4));
double a6 = 1.0 / 720.0 * N*sin(B)*pow(cos(B), 5)*(61 - 58 * pow(t, 2) + pow(t, 4));
return (a2*pow(l, 2) + a4*pow(l, 4) + a6*pow(l, 6));
}
double Fy(double B, double l)
{
double N = a / sqrt(1 - e_square*pow(sin(B), 2)); //卯酉圈曲率半径
double t = tan(B);
double fai = e1*cos(B);
double a3 = 1.0 / 6.0 * N*pow(cos(B), 3)*(1 - pow(t, 2) + pow(fai, 2));
double a5 = 1.0 / 120.0 * N*pow(cos(B), 5)*(5 - 18 * pow(t, 2) + pow(t, 4) + 14 * pow(fai, 2) - 58 * pow(fai, 2)*pow(t, 2));
return a3*pow(l, 3) + a5*pow(l, 5);
}
void InverseCalculation(double x, double y, double &B, double &a1, double &l)
{
B = (x - Fx(B) - Fx(B, l)) / c*beita0;
a1 = a*cos(B) / sqrt(1 - e_square*pow(sin(B), 2));
l = (y - Fy(B, l)) / a1;
}
//坐标反算
void InverseCal(double x, double y, double &B, double&l)
{
//初始值设置
double B0 = (x / c*beita0);
double a10 = a*cos(B0) / sqrt(1 - e_square*pow(sin(B0), 2));
double l0 = (y / a10);
printf("设置初始值为：\n B0 = %lf l0 = %lf\n", B0, l0);
B = B0;
double a1 = a10;
l = l0;
int n = 1;
InverseCalculation(x, y, B, a1, l);
printf("第%d次迭代结果为： B = %lf l = %lf\n", n, B, l);
while (fabs(B0 - B) > 0.0001 || fabs(l0 - l) > 0.0001)
{
n++; //次数
B0 = B;
a10 = a1;
l0 = l;
InverseCalculation(x, y, B, a1, l);
printf("第%d次迭代结果为： B = %lf l = %lf\n", n, B, l);
}
}
/*------------------------------------------------------------------------------*/
/*带号转换*/
void DegreeTrans(double &B, double &l, int &n, double &L, int degree)
{
switch (degree)
{
case 6:
{
int centerLongitude = n * 6 - 3; //中央经度
B = B * 206265 / 3600;
l = l * 206265 / 3600;
L = l + centerLongitude;
break;
}
case 3:
{
n = (1.5 + L) / 3; //计算带号
int centerLongitude = n * 3; //中心经度
l = L - centerLongitude; //经差
l = l * 3600 / 206265;
B = B * 3600 / 206265;
break;
}
}
}
/*坐标正算：大地坐标转化为高斯投影坐标*/
void Calculate(double B, double l, double &x, double &y)
{
double N = a / sqrt(1 - e_square*pow(sin(B), 2)); //卯酉圈曲率半径
double t = tan(B);
double fai = e1*cos(B);
double X = c*(beita0*B + (beita2*cos(B) + beita4*pow(cos(B), 3) + beita6*pow(cos(B), 5) + beita8*pow(cos(B), 7))*sin(B)); //赤道至纬度B的子午线弧长
x = X + pow(l, 2)*N*sin(B)*cos(B) + pow(l, 4) / 24 * N*sin(B)*pow(cos(B), 3)*(5 - pow(t, 2) + 9 * pow(fai, 2) + 4 * pow(fai, 4)) +
pow(l, 6) / 720 * N*sin(B)*pow(cos(B), 5)*(61 - 58 * pow(t, 2) + pow(t, 4));
y = l*N*cos(B) + pow(l, 3) / 6 * N*pow(cos(B), 3)*(1 - pow(t, 2) + pow(fai, 2)) +
pow(l, 5) / 120 * N*pow(cos(B), 5)*(5 - 18 * pow(t, 2) + 14 * pow(fai, 2) - 58 * pow(fai, 2)*pow(t, 2));
}
int main()
{
printf("*********高斯投影带转换（6°——3°）************\n\n");
double x = 3589544.286; //初始高斯坐标xy
double y = 20679136.438;
printf("需要转换的6°带坐标为：\n x = %lf y = %lf\n\n", x, y);
int n = 20;
y = y - n * 1000000 - 500000; //转换成自然值
printf("计算过程：\n（一）通用值转换成自然值\n x = %lf y = %lf\n", x, y);
printf("\n（二）高斯投影坐标反算\n");
double B = 0;
double l = 0;
InverseCal(x, y, B, l); //坐标反算
printf("迭代结果为：\n B = %lf l = %lf\n", B, l);
double L = 0; //定义大地经度
DegreeTrans(B, l, n, L, 6); //求中央经线及大地坐标
printf("\n（三）求解6°带中央经线\n转换后的大地坐标为：\n B = %lf L = %lf\n", B, L);
printf("\n（四）求解3°带的经差\n");
int n1 = 0; //3°带代号
double l1 = 0; //经差
DegreeTrans(B, l1, n1, L, 3); //求中央经线
printf("转换后的3°带的经度纬度和代号为：\n B = %lf l = %lf n=%d\n", B, l1, n1);
printf("\n（五）高斯投影坐标正算\n");
double x1 = 0, y1 = 0;
Calculate(B, l1, x1, y1); //坐标正算
printf("正算结果为：\n x = %lf y = %lf\n", x1, y1);
y1 = y1 + n1 * 1000000 + 500000; //高斯投影坐标y
printf("\n（六）自然值转换成通用值\n x = %lf y = %lf\n", x1, y1);
return 0;
}注：代码是根据中国矿业大学出版社出版的《应用大地测量学》高斯投影一章的知识编写的，为了理解书本内容，精度没经过验证。
3.3 输出结果
注：以下是编译器输出的内容，未经修改
需要转换的6°带坐标为：
x = 3589544.286000 y = 20679136.438000
计算过程：
（一）通用值转换成自然值
x = 3589544.286000 y = 179136.438000
（二）高斯投影坐标反算
设置初始值为：
B0 = 0.558058 l0 = 0.033078
第1次迭代结果为： B = 0.560073 l = 0.033116
第2次迭代结果为： B = 0.560076 l = 0.033116
迭代结果为：
B = 0.560076 l = 0.033116
（三）求解6°带中央经线
转换后的大地坐标为：
B = 32.090028 L = 118.897435
（四）求解3°带的经差
转换后的3°带的经度纬度和代号为：
B = 0.560076 l = -0.019243 n=40
（五）高斯投影坐标正算
正算结果为：
x = 3552710.678780 y = -104087.421312
（六）自然值转换成通用值
x = 3552710.678780 y = 40395912.578688
总结：测量学与其他学科的不同就是对结果的精度近乎苛刻（尤其是大地测量学），在实际应用中还是需要精度验证的。
