批处理就是文件形式读取数据,四个子程序可以不同函数放一起
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设计资料,需要的话私信可以发你邮箱