求高斯坐标正反算，要求如图
批处理就是文件形式读取数据，四个子程序可以不同函数放一起
ssuididjhshjjdidudh
没人的话我来做
以下是Matlab实现高斯正反算的两个方法：
方法1，教程上提供的的高斯正反算例程，函数1：高斯正算
function PRO = GKZS(Lat,Lon,L0,a,f)
% Lat: Latitude(rad)
% Lon: longitude(rad)
% REF//程鹏飞,成英燕,文汉江,等.2000国家大地坐标系实用宝典[M].
% //北京:测绘出版社,2008:144-148.
Lat = Lat*pi/180;
Lon = Lon*pi/180;
MedLon = L0*pi/180; %中央子午线经度
Eth.R0 = a;%长半轴
Eth.f = f; % 扁率
Eth.e12 = 2*Eth.f - Eth.f*Eth.f; % 第一偏心率额e^2
Eth.e22 = Eth.e12/((1 - Eth.f)*(1 - Eth.f));% 第二偏心率额e'^2
%% 高斯投影正算公式
RN = Eth.R0/sqrt(1 - Eth.e12*sin(Lat)*sin(Lat));
Lon = Lon - MedLon;
Lon2 = Lon*Lon;
Lon4 = Lon2*Lon2;
tnLat = tan(Lat);
tn2Lat = tnLat*tnLat;
tn4Lat = tn2Lat*tn2Lat;
csLat = cos(Lat);
cs2Lat = csLat*csLat;
cs4Lat = cs2Lat*cs2Lat;
Eta2 = Eth.e22*cs2Lat;
NTBLP = RN*tnLat*cs2Lat*Lon2;
coe1 = (5 - tn2Lat + 9*Eta2 + 4*Eta2*Eta2)*cs2Lat*Lon2/24;
coe2 = (61 - 58*tn2Lat + tn4Lat)*cs4Lat*Lon4/720;
x = Merdian(Eth,Lat) + NTBLP*(0.5 + coe1 + coe2);
NBLP = RN*csLat*Lon;
coe3 = (1 - tn2Lat + Eta2)*cs2Lat*Lon2/6;
coe4 = (5 - 18*tn2Lat + tn4Lat + 14*Eta2 - 58*tn2Lat*Eta2)*cs4Lat*Lon4/120;
y = NBLP*(1 + coe3 + coe4) + 500000;
PRO=[x y]
end
function X0 = Merdian(Eth,Lat)
% REF//过家春.子午线弧长公式的简化及其泰勒级数解释[J].测绘学报,2014,43(2):125-130.
S0 = Eth.R0*(1 - Eth.e12);
e2 = Eth.e12;
e4 = e2*e2;
e6 = e4*e2;
e8 = e6*e2;
e10 = e8*e2;
e12 = e10*e2;
A1 = 1 + 3*e2/4 + 45*e4/64 + 175*e6/256 + 11025*e8/16384 + 43659*e10/65536 + 693693*e12/1048576;
B1 = 3*e2/8 + 15*e4/32 + 525*e6/1024 + 2205*e8/4096 + 72765*e10/131072 + 297297*e12/524288;
C1 = 15*e4/256 + 105*e6/1024 + 2205*e8/16384 + 10395*e10/65536 + 1486485*e12/8388608;
D1 = 35*e6/3072 + 105*e8/4096 + 10395*e10/262144 + 55055*e12/1048576;
E1 = 315*e8/131072 + 3465*e10/524288 + 99099*e12/8388608;
F1 = 693*e10/1310720 + 9009*e12/5242880;
G1 = 1001*e12/8388608;
X0 = S0*(A1*Lat - B1*sin(2*Lat) + C1*sin(4*Lat) - D1*sin(6*Lat) +...
E1*sin(8*Lat) - F1*sin(10*Lat) + G1*sin(12*Lat));
end
函数2：高斯反算
function BL = GKFS(x,y,L0,a,f)
iPI = 0.0174532925199433; %3.1415926535898/180.0;
Eth.a = a ;%长半轴
Eth.f = f;% 扁率
e2 = 2*Eth.f-Eth.f*Eth.f;
e1 = (1.0-sqrt(1-e2))/(1.0+sqrt(1-e2));
ee = e2/(1-e2);
X = y
Y = x
longitude0 = L0*pi/180; %中央子午线经度
%
ZoneWide = 6; %6度带宽
ProjNo = fix(X/1000000) ; %查找带号
X0 = ProjNo*1000000+500000;
Y0 = 0;
xval = X-X0; yval = Y-Y0;%带内大地坐标
M = yval;
u = M/(Eth.a*(1-e2/4-3*e2*e2/64-5*e2*e2*e2/256));
fai = u+(3*e1/2-27*e1*e1*e1/32)*sin(2*u)+(21*e1*e1/16-55*e1*e1*e1*e1/32)*sin(4*u)+(151*e1*e1*e1/96)*sin(6*u)+(1097*e1*e1*e1*e1/512)*sin(8*u);
C = ee*cos(fai)*cos(fai);
T = tan(fai)*tan(fai);
NN = Eth.a/sqrt(1.0-e2*sin(fai)*sin(fai));
R = Eth.a*(1-e2)/sqrt((1-e2*sin(fai)*sin(fai))*(1-e2*sin(fai)*sin(fai))*(1-e2*sin(fai)*sin(fai)));
D = xval/NN;
%计算经度(Longitude) 纬度(Latitude)
longitude1 =longitude0+(D-(1+2*T+C)*D*D*D/6+(5-2*C+28*T-3*C*C+8*ee+24*T*T)*D*D*D*D*D/120)/cos(fai);
latitude1 = fai -(NN*tan(fai)/R)*(D*D/2-(5+3*T+10*C-4*C*C-9*ee)*D*D*D*D/24+(61+90*T+298*C+45*T*T-256*ee-3*C*C)*D*D*D*D*D*D/720);
%转换为度 DD
L = longitude1 / iPI
B = latitude1 / iPI;
L = fix(L)+fix((L-fix(L))*60)/100+((L-fix(L))*60-fix((L-fix(L))*60))*60/10000
B = fix(B)+fix((B-fix(B))*60)/100+((B-fix(B))*60-fix((B-fix(B))*60))*60/10000
BL = [B L]
end
方法2的高斯正算函数：
function [ x,y ] = BL2XY( B,L,a,e,k )
ee=e/sqrt(1-e^2);
% 计算中央子午线经度
if k==6
L0=6*(fix(L/6)+1)-3;
elseif k==3
L0=3*(fix((L-1.5)/3)+1);
end
LP=(L-L0)/180*pi;
B=B/180*pi;
% 计算子午圈弧长
ap=1+3/4*e^2+45/64*e^4+175/256*e^6+11025/16384*e^8+43659/65536*e^10;
bp=3/4*e^2+15/16*e^4+525/512*e^6+2205/2048*e^8+72765/65536*e^10;
cp=15/64*e^4+105/256*e^6+2205/4096*e^8+10395/16384*e^10;
dp=35/512*e^6+315/2048*e^8+31185/131072*e^10;
ep=315/16384*e^8+3465/65536*e^10;
fp=639/131072*e^10;
X=a*(1-e^2)*(ap*B-bp/2*sin(2*B)+cp/4*sin(4*B)-dp/6*sin(6*B)+ep/8*sin(8*B)-fp/10*sin(10*B));
% 计算高斯平面坐标
N=a/sqrt(1-e^2*sin(B)*sin(B));
eta=ee*cos(B);
t=tan(B);
m=cos(B)*LP;
x=X+N*t*(0.5*m^2+1/24*(5-t^2+9*eta^2+4*eta^4)*m^4+1/720*(61-58*t^2+t^4)*m^6);
y=N*(m+1/6*(1-t^2+eta^2)*m^3+1/120*(5-18*t^2+t^4+14*eta^2-58*eta^2*t^2)*m^5);
end
高斯反算函数：
function [ B,L ] = XY2BL( x,y,a,e,L0 )
% CGCS2000椭球参数
% a=6378137;
% e=0.081819191042811;
% ee=0.0820944381518626;
% 1975国际椭球参数
% a=6378140;
% e=0.0818192214555235;
ee=e/sqrt(1-e^2);
% 计算底点的经度
ap=1+3/4*e^2+45/64*e^4+175/256*e^6+11025/16384*e^8+43659/65536*e^10;
bp=3/4*e^2+15/16*e^4+525/512*e^6+2205/2048*e^8+72765/65536*e^10;
cp=15/64*e^4+105/256*e^6+2205/4096*e^8+10395/16384*e^10;
dp=35/512*e^6+315/2048*e^8+31185/131072*e^10;
ep=315/16384*e^8+3465/65536*e^10;
fp=639/131072*e^10;
B=x/(a*(1-e^2)*ap);
deltaB=1000;
while deltaB>10^(-20)
Bi=1/ap*(x/(a*(1-e^2))+bp/2*sin(2*B)-cp/4*sin(4*B)+dp/6*sin(6*B)-ep/8*sin(8*B)+fp/10*sin(10*B));
deltaB=abs(Bi-B);
B=Bi;
end
Bf=B;
% 计算大地经纬度BL
tf=tan(Bf);
Vf=sqrt(1+ee^2*cos(Bf)*cos(Bf));
Wf=sqrt(1-e^2*sin(Bf)*sin(Bf));
Nf=a/Wf;
etaf=ee*cos(Bf);
B=Bf-1/2*Vf^2*tf*((y/Nf)*(y/Nf)-1/12*(5+3*tf^2+etaf^2-9*etaf^2*tf^2)*(y/Nf)^4+1/360*(61+90*tf^2+45*tf^4)*(y/Nf)^6);
LP=1/cos(Bf)*((y/Nf)-1/6*(1+2*tf^2+etaf^2)*(y/Nf)^3+1/120*(5+28*tf^2+24*tf^4+6*etaf^2+8*etaf^2*tf^2)*(y/Nf)^5);
B=B/pi*180;
LP=LP/pi*180;
L=L0+LP;
end
高斯正反算批量计算及GUI设计资料，需要的话私信可以发你邮箱
//高斯坐标的正反算
#include<iostream>
#include<cmath>
#define PI 3.1415
using namespace std;
//正算 x=3380330.773 y=320089.9761
//B=30° 30′ L=114°20′
void zhensuan()
{
double d, f, m;
cout << "输入B(°′″): " << endl;
cin >> d >> f >> m;
double dl, fl, ml;
cout << "输入L(°′″): " << endl;
cin >> dl >> fl >> ml;
double _b = d + f / 60 + m / 3600;
double B = (PI/180)*(_b);//角度要变为弧度 PI/180*B
double L = dl + fl / 60 + ml / 3600;
double c = 6399698.9017827110;
double e = 0.006738525414683;//e"^2
double m_e = 0.006693421622966;//e^2
double V = sqrt(1 + e * cos(B)*cos(B));//v=根号下1+e'^2*cosB*cos(B)
double p = 206264.806247096355;//p"
double t = tan(B);
double n = e * cos(B) * cos(B);//n^2
double n6;
if ((int)L % 6 == 0) {
n6 = floor(L / 6) ;//按6度带划分
}
else n6 = floor(L / 6)+1;
double lm = n6 * 6 - 3;//中央子午线度数
double l = fabs(L - lm)*3600;//l"以秒为单位 经差（p点经度和中央子午线之差）
double a = 6378245;//单位m
double b = 6356863.02;//单位m
double m_p = 1 / p;//p的倒数
double N = c / V;
//先求X
double m1 = N * cos(B);
double m2 = (N / 2) * sin(B) * cos(B);
double m3 = (N / 6)* cos(B) * cos(B) * cos(B) * (1 - t * t + n );
double m4 = (N / 24) * sin(B) * cos(B) * cos(B) * cos(B) * (5 - t * t + 9 * n + 4 * n * n);
double m5 = (N / 120) * cos(B) * cos(B) * cos(B) * cos(B) * cos(B) * (5 - 18 * t * t + t * t * t * t + 14 * n - 58 * n * t * t);
double m6 = (N / 720) * sin(B) * cos(B) * cos(B) * cos(B) * cos(B) * cos(B) * (61 - 58 * t * t + t * t * t * t);
double X = 111134.861 * (_b)-32005.780 * sin(B) * cos(B) - 133.929 * sin(B) * sin(B) * sin(B) * cos(B) - 0.698 * sin(B) * sin(B) * sin(B) * sin(B) * sin(B)*cos(B);
//_b是度数哦
double x = X + (m2 * m_p * m_p * l * l) + (m4 * m_p * m_p * m_p * m_p * l * l * l * l) + (m6 * m_p * m_p * m_p * m_p * m_p * m_p * l * l * l * l * l * l);
double y = (N / p) * cos(B) * l + (N / (6 * p * p*p)) * cos(B) * cos(B) * cos(B) * (1 - t * t + n ) * l * l * l + (N / (120 * p * p * p * p * p)) * cos(B) * cos(B) * cos(B) * cos(B) * cos(B) * (5 - 18 * t * t + t * t * t * t + 14 * n - 58 * n * t * t) * l * l * l * l * l;
//输出
cout << "正算结果：" << endl;
cout << "x= " << x << endl;
cout << "y= " << y << endl;
}
//坐标反算 x=3380330.773 y=320089.9761
//B=30° 30′ L=114°20′
//由x,y求B，L
void fansuan()
{
cout << "输入坐标x,y: " << endl;
double x, y;
cin >> x >> y;
//迭代解法计算底点纬度
double p = x / 6367588.4969;//弧度
double Bd = p + (50221746 + (293622 + (2350 + 22 * cos(p) * cos(p)) * cos(p)) * cos(p)) * (1e-10) * sin(p) * cos(p);//底点坐标
double c = 6399698.9017827110;
double e = 0.006738525414683;//e"^2
double m_e = 0.006693421622966;//e^2
double Vd = sqrt(1 + e * cos(Bd));//v=根号下1+e'^2*cosB
double td = tan(Bd);
double nd = e * cos(Bd) * cos(Bd);
double Nd = c / Vd;
double Md = c / (Vd * Vd * Vd);
double n1 = 1 / (Nd * cos(Bd));
double n2 = (td / (2 * Md * Nd));
double n3 = (1 + 2 * td * td + nd * nd) / (6 * Nd * Nd * Nd * cos(Bd));
double n4 = (td * (5 + 3 * td * td + nd * nd - 9 * nd * nd * td * td)) / (24 * Md * Nd * Nd * Nd);
double n5 = (5 + 28 * td * td + 24 * td * td * td * td + 6 * nd * nd + 8 * nd * nd * td * td) / (120 * Nd * Nd * Nd * Nd * Nd * cos(Bd));
double n6 = (td * (61 + 90 * td * td + 45 * td * td * td * td)) / (720 * Md * Nd * Nd * Nd * Nd * Nd);
double B = Bd - n2 * y * y + n4 * y * y * y * y-n6* y * y * y * y*y * y;
double L = n1 * y - n3 * y * y * y + n5 * y * y * y * y * y;
cout << "B= " << B << endl;//按弧度计算
cout << "L= " << L << endl;
cout << "转换成角度后：B=" << (180 / PI) * B << endl;
cout << "转换成角度后：L= " << (180 / PI) * L << endl;//这里是经差，因为高斯坐标和国家统一坐标不同，不知其带号
}
int main() {
//克拉索夫斯基椭球
//a=6378245m b=6356863.0187730473m
//c=6399698.9017827110m a=1/298.3
//e^2=0.006693421622966 e'^2=0.006738525414683
zhensuan();
fansuan();
}