convertxytobl.cpp

地理信息系统中把坐标转换到经纬度的程序,一个课程设计的课题.

#include <iostream>
#include <math.h>
#define pi 3.1415926;

using namespace std;

/* 
功能：把高斯克吕格XY坐标转化为经纬度BL
参数[IN]：a,b -- 椭球体参数 
x,y -- 高斯克吕格XY坐标
T --  分带情况
参数[OUT]：B,L -- 经纬度BL
n   -- 所在带号
l0  -- 所在带的中央经度
*/
void ConvertXYtoBL(double a,double b,double X,double Y,int T,double& B,double& L,int& n,int& l0);
void ConvertBLtoXY(double a,double b,double B,double L,int T, double& X,double& Y,int& n,int& l0);
void ConvertxytoXY(double x, double y,double inputMax[3][2],double &X, double &Y);
bool nizhen(double *inptr,double (*outptr)[3], int n);
void matrixMulti(double inputMax1[3][3],double inputMax2[3][2],double (*outptr)[2]);
void ConvertxytoBL(double x, double y,double inputMax[3][2],double &B, double &L);

void main()
{   
	
	int nMuchConvert = 1;
	while(1 == nMuchConvert)
	{
// 		double Y0=0.0,Y1=0.0,Y2=0.0;					//控制点在高斯坐标的横坐标，以米为单位
// 		double X0=0.0,X1=0.0,X2=0.0;					//控制点在高斯坐标的纵坐标，以米为单位
// 		int T=0;					//输入参数：选择分带标准
// 		
// 		double x0=0.0,x1=0.0,x2=0.0;        //输入参数：控制点在自由坐标系下的横坐标，以米为单位
// 		double y0=0.0,y1=0.0,y2=0.0;        //输入参数：控制点在自由坐标系下的纵坐标，以米为单位
		
// 		double L0=0.0,L1=0.0,L2=0.0,L=0.0;					//输入参数: 经度坐标，其中经度坐标以度为单位
// 		double B0=0.0,B1=0.0,B2=0.0,B=0.0;					//输入参数: 纬度坐标，其中纬度坐标以度为单位
		double L=0.0;					//输入参数: 经度坐标，其中经度坐标以度为单位
		double B=0.0;	
		//用户输入自由坐标系的三个控制点纬度坐标
// 		cout << "请输入控制点1的纬度,经度，在自由坐标系下的横坐标，纵坐标： ";
// 		cin >> B0 >> L0 >> y0 >> x0;
// 		cout << "\n";
// 		
// 		cout << "请输入控制点2的纬度,经度，在自由坐标系下的横坐标，纵坐标： ";
// 		cin >> B1 >> L1>> y1 >> x1;
// 		cout << "\n";
// 		
// 		cout << "请输入控制点3的纬度,经度，在自由坐标系下的横坐标，纵坐标： ";
// 		cin >> B2 >> L2>> y2 >> x2;
// 		cout << "\n";
// 		
// 		double temp;
// 		cout << "请选择分带(0 -- 6°分带，1 -- 3°分带): ";
// 		cin >> temp;
// 		
// 		if(temp!=0 && temp!=1)
// 		{
// 			cout<< "出错"<<"\n";
// 			return;
// 		}
// 		T = temp;
// 		cout << "\n";
// 		
// 		temp = -1;
// 		cout << "请选择球体坐标系:"
// 			<<"\n"<<"(0 -- 1954北京坐标系"
// 			<<"\n"<<" 1 -- 1980西安坐标系"
// 			<<"\n"<<" 2 --  WGS84坐标系  ): ";
// 		cin >> temp;
// 		
// 		double a=0.0, b=0.0;		//椭球长短半轴
// 		if(0 == temp)
// 		{
// 			a = 6378245;
// 			b = 6356863.02;
// 		}
// 		else if(1 == temp)
// 		{
// 			a = 6378140;
// 			b = 6356755.3;
// 		}
// 		else if(2 == temp)
// 		{
// 			a = 6378137;
// 			b = 6356752.3;
// 		}
// 		else
// 		{
// 			cout<< "出错"<<"\n";
// 			return;
// 		}
// 		cout << "\n";
// 	
// 		int l0=0, n=0;  
		//由控制点的经纬度坐标求出点在某一标准坐标系下的横纵坐标X,Y
	//	ConvertBLtoXY(a, b, B0, L0, T, X0, Y0, n, l0);
	//	ConvertBLtoXY(a, b, B1, L1, T, X1, Y1, n, l0);
	//	ConvertBLtoXY(a, b, B2, L2, T, X2, Y2, n, l0);
		
// 		double pointx0[]={B0,L0,y0,x0};
// 		double pointx1[]={B1,L1,y1,x1};
// 		double pointx2[]={B2,L2,y2,x2};

		double pointx0[]={30.521, 114.351,  3549063.4 , 12729710.7};
		double pointx1[]={30.731, 114.612 , 3576095.6, 12758765.6};
		double pointx2[]={30.894, 114.542 , 3597118.7 ,  12750973.1};
		
		//构造矩阵1
		double matrixOne[3][3];
		matrixOne[0][0]=1;
		matrixOne[0][1]=pointx0[2]+pointx1[2]+pointx2[2];
		matrixOne[0][2]=pointx0[3]+pointx1[3]+pointx2[3];
		matrixOne[1][0]=pointx0[2]+pointx1[2]+pointx2[2];
		matrixOne[1][1]=pow(pointx0[2],2)+pow(pointx1[2],2)+pow(pointx2[2],2);
		matrixOne[1][2]=pointx0[2]*pointx0[3]+pointx1[2]*pointx1[3]+pointx2[2]*pointx2[3];
		matrixOne[2][0]=pointx0[3]+pointx1[3]+pointx2[3];
		matrixOne[2][1]=pointx0[2]*pointx0[3]+pointx1[2]*pointx1[3]+pointx2[2]*pointx2[3];
		matrixOne[2][2]=pow(pointx0[3],2)+pow(pointx1[3],2)+pow(pointx2[3],2);
		
		//构造矩阵2
		double matrixTwo[3][2];
		matrixTwo[0][0]=pointx0[1]+pointx1[1]+pointx2[1];
		matrixTwo[0][1]=pointx0[0]+pointx1[0]+pointx2[0];
		matrixTwo[1][0]=pointx0[2]*pointx0[3]+pointx1[2]*pointx1[3]+pointx2[2]*pointx2[3];
		matrixTwo[1][1]=pointx0[2]*pointx0[1]+pointx1[2]*pointx1[1]+pointx2[2]*pointx2[1];
		matrixTwo[2][0]=pointx0[3]*pointx0[2]+pointx1[3]*pointx1[2]+pointx2[3]*pointx2[2];
		matrixTwo[2][1]=pointx0[3]*pointx0[0]+pointx1[3]*pointx1[0]+pointx2[3]*pointx2[0];
		
		//保存计算得到的matrixOne的逆矩阵
		double matrixOutput[3][3]=
		{
			{0, 0, 0},
			{0, 0, 0},
			{0, 0, 0}
		};
		
		//进行逆矩阵计算
		nizhen(matrixOne[0],matrixOutput,3);
		//int i=0,j=0;
		
        //求矩阵的乘积
		double mutiMatrix[3][2]=
		{
			{0, 0},
			{0, 0},
			{0, 0}
		};
		matrixMulti(matrixOutput, matrixTwo, mutiMatrix);
		
		
		
		//输入自由坐标系下任意点的横，纵坐标x,y
		double x=0.0,y=0.0;
		double X=0.0,Y=0.0;
		
		cout << "请输入自由坐标系下某点的横坐标y: ";
		cin >> y;
		cout << "\n";
		
		cout << "请输入自由坐标系下某点的纵坐标x: ";
		cin >> x;
		cout << "\n";
		//转换成某一标准坐标系下的横，纵坐标X,Y
	//	ConvertxytoXY(x, y, mutiMatrix,X, Y);
		
		//某一标准坐标系下的横纵坐标向经纬度转化
	//	ConvertXYtoBL(a, b, X, Y, T, B, L, n, l0);
	
		ConvertxytoBL(x, y, mutiMatrix,B, L);
		//输出该点的经纬度
// 		cout << "该点所在带号是：" << n << "\n";
// 		cout << "该点所在带中央经线是: " << l0 << "\n";
		cout << "该点的经度是：" << L << "\n";
		cout << "该点的纬度是: " << B << "\n";
		
		cout << "是否继续处理？(1 -- 是，非1 -- 否): " << "\n";
		cin >> nMuchConvert;
	}
}

// 由点的高斯坐标得到点的经纬度
// void ConvertXYtoBL(double a,double b,double X,double Y,int T, double& B,double& L,int& n,int&l0 )
// {
// 	double yy;   
// 	
// 	n = (int)(X / 1000.0 / 1000.0)  ;	// 查找带号
// 	yy =  (n * 1000.0 + 500.0) * 1000.0;	
// 	if(0 == T)
// 	{
// 		l0 = n * 6 - 3;							//点所在带的中央经线的经度
// 	}
// 	if(1 == T)
// 		l0 = n * 3 ;
// 	
// 	double e, e1, e2, alpha, Nf, Mf, Tf, Cf, D, Bf, Rf;
// 	
// 	// 离心率
// 	e = sqrt(1 - b / a * b / a);			// 这里不能写成b * b / a * a,否则会超出double范围
// 	e1 = (1 - b / a) / (1 + b / a);		
// 	e2 = sqrt(a / b * a / b - 1);
// 	
// 	// 计算过程中的参数
// 	Mf =Y;
// 	alpha = Mf / (a * (1 - e * e / 4 - 3 * pow(e,4) / 64 - 5 * pow(e,6) / 256));
// 	Bf = alpha + (3 * e1 / 2 - 27 * pow(e1,3) / 32) * sin(2 * alpha) 
// 		+ (21 * e1 * e1 / 16 - 55 * pow(e1,4) / 32) * sin(4 * alpha) 
// 		+ (151 * pow(e1,3) / 96) * sin(6 * alpha);
// 	Tf = tan(Bf) * tan(Bf);
// 	Cf = e2*e2 * cos(Bf)* cos(Bf);
// 	Nf = a / sqrt(1 - e*e/sin(Bf)*sin(Bf));
// 	D = (X-yy )/Nf;
// 	Rf = a*(1-e*e)/sqrt(pow((1-e*e*sin(Bf)*sin(Bf)),3));
// 	
// 	// 求出BL
// 	B = Bf - Nf*tan(Bf)/Rf * (D*D/2 - (5 + 3*Tf + Cf - 9*Tf*Cf) * pow(D,4)/24 + (61 + 90*Tf+45 *Tf*Tf)*pow(D,6)/720);
// 	B = B * 180.0 / pi;
// 	
// 	L =  (D -(1+2*Tf +Cf)*pow(D,3)/6 +(5+28*Tf+6*Cf+8*Tf*Cf+24*Tf*Tf)*pow(D,5)/120)/cos(Bf);
// 	L = l0 + L * 180.0 /pi;
// }
// 
// void ConvertBLtoXY(double a,double b,double B,double L,int T, double& X,double& Y,int& n,int& l0)
// {
// 	double FE;   
// 	
// 	
// 	
// 	if(T==0)
// 	{
// 		n = (int)(L / 6) +1 ;	// 查找带号
// 		l0 = n * 6 - 3;	//点所在带的中央经线的经度	
// 	}
// 	if(T==1)
// 	{
// 		n = (int)(L / 3) +1 ;	// 查找带号
// 		l0 = n * 3 ;
// 	}
// 	FE =  500000 +n*1000000;
// 	double e, e2, K, C, A, M, N;
// 	
// 	// 离心率
// 	e = sqrt(1 - b / a * b / a);			// 这里不能写成b * b / a * a,否则会超出double范围
// 	
// 	e2 = sqrt(a / b * a / b - 1);
// 	
// 	// 计算过程中的参数
// 	B= B / 180 * pi; 
// 	L = L / 180 * pi;
// 	K=tan(B)*tan(B);
// 	C=e2*e2 * cos(B)*cos(B);
// 	double tempL = l0 / 180.0 * pi;
// 	A=(L- tempL)*cos(B);
// 	
// 	M=a*(1-e*e/4 -3*pow(e,4)/64 -5*pow(e,6)/256)*B-a*(3*pow(e,2)/8 + 3*pow(e,4)/32 +45*pow(e,6)/1024)*sin(2*B) + a*(15*pow(e,4)/256 + 45*pow(e,6)/1024)*sin(4*B)-a*35*pow(e,6)/3072*sin(6*B);
// 	N=a*a/b/sqrt(1 + e2*e2*cos(B)*cos(B));
// 	// 求出XY
// 	Y=M + N*tan(B)*(A*A/2 + (5 - K +9*C+4*C*C)*pow(A,4)/24) +(61-58*K +K*K +270*C-330*K*C)*pow(A,6)/720;
// 	X=FE + N*(A + (1-K+C)*pow(A,3)/6 + (5-18*K +K*K +14*C -58*K*C)*pow(A,5)/120);
// 	
// 	
// 	
// 	}

//将点的坐标由自由坐标系转换到标准坐标系下
// void ConvertxytoXY(double x, double y,double inputMax[3][2],double &X, double &Y)
// {
// 	//构造参数矩阵
// 	double A=0.0, B=0.0, C=0.0, D=0.0, E=0.0, F=0.0;
// 	
// 	A = inputMax[1][0];
// 	B = inputMax[2][0];
// 	C = inputMax[0][0];
// 	D = inputMax[1][1];
// 	E = inputMax[2][1];
// 	F = inputMax[0][1];
// 	
// 	X=D*y+ E*x+ F;
// 	Y=A*y+ B*x+ C;
// 	
// }

void ConvertxytoBL(double x, double y,double inputMax[3][2],double &B, double &L)
{
	//构造参数矩阵
	double A=0.0, M=0.0, C=0.0, D=0.0, E=0.0, F=0.0;
	
	A = inputMax[1][0];
	M = inputMax[2][0];
	C = inputMax[0][0];
	D = inputMax[1][1];
	E = inputMax[2][1];
	F = inputMax[0][1];
	
	B=D*y+ E*x+ F;
	L=A*y+ M*x+ C;
	
}

double jshls(double *inptr, int n);//计算方阵行列式  inptr方阵第一元素指针 n阶数 
double dsyzs(double *inptr, int n, int i, int j);//求代数余子式的值 i,j为1~N符合数学书的习惯
void yzs(double *inptr, double *outptr, int n, int i, int j);//余子式
//计算逆矩阵
bool nizhen(double *inptr,double (*outptr)[3], int n)
{
    int i=0,j=0,k=0;
    double tmpa;//存储行列式|A|的值 
    tmpa=jshls(inptr, n);
    if(tmpa==0)
    {
        //printf("充要条件（|A|<>0）不成立");
        //交给上层显示，为了以后移植用，因为在GUI编程中printf会造成出错 
        return false;
    }
	for(i=0;i<n;i++)
	{
		for(j=0;j<n;j++)
		{
			//求逆阵各元素的值，包含了伴随方阵的运算 
			outptr[i][j]=dsyzs(inptr, n, i+1, j+1)/tmpa; //交换ij，看伴随方阵的定义 
		}
	}
	return true;
}
double jshls(double *inptr, int n)
{
    double tmp=0.0,tmp1=0.0,tmp2=0.0;
    int j;
    tmp=0;
    if(n==1) return *inptr;
    for(j=0; j<n; j++) tmp += (*(inptr+j))*dsyzs(inptr, n, 1, j+1);//按第一行展开 
    return tmp;
}
double dsyzs(double *inptr, int n, int i, int j)
{
    double *yzsptr, yzsjg;//余子式的 指针 临时指针 计算结果     
    yzsptr = (double*)malloc(((n-1)*(n-1)*sizeof(double)));//分配n-1阶方阵的内存 
    if( yzsptr==NULL )
    {
        printf("无法分配内存，代号001，任意键退出");//GUI的时候用MESSAGEBOX对话框替换掉 
        exit(-1);
    }
    if(n<2)
    {
        printf("参数：余子式错误。n大于或等于2。");//GUI的时候用MESSAGEBOX对话框替换掉 
        exit(-1);
    }
    yzs(inptr, yzsptr, n, i, j);
    yzsjg=jshls(yzsptr, n-1);
    free(yzsptr);//释放n-1阶方阵的内存
    return (((i+j)%2)?(-1):1)*yzsjg;//返回代数余子式的结果，代数余子式和余子式是不同的哦 
}
void yzs(double *inptr, double *outptr, int n, int i, int j)
{
    int k;
    double *tmptr;
    tmptr =outptr;
    for(k=0; k<(n*n); k++)
    {
        if( (k%n == (j-1)) || ((k-k%n)/n == (i-1))) continue;
        *tmptr = *(inptr+k);
        tmptr++;
    }
}


// 计算矩阵的乘积
void matrixMulti(double inputMax1[3][3],double inputMax2[3][2],double (*outptr)[2])
{
	double t=0.0;
	int i,j,k,temp=0;
	for(i=0;i<3;i++)
		for(j=0;j<2;j++){
			for(k=0;k<3;k++)
				outptr[i][j] +=inputMax1[i][k]*inputMax2[k][j];
		}
}