matlab大地测量高斯投影正反算程序设计实验

1、-1-高斯投影正反算一、实验目的编写高斯投影正反算的程序，并对格式化文件数据进行计算，验证程序。二、实验内容：1、高斯投影正算公式己知大地坐标(B,L)及中央子午线经度，计算高斯平面坐标(x,y),公式如下：x=X+sin(B)cos(B)l2+sin(B)cos3(B)(5-12+9+4)1224+舟sin(B)cos5(B)(61-58t3+t4+270-330干)1y=Ncos(B)l+cos3(B)(l-t2+n2)l3+-cos5(B)(5-18t2+14+14n2-58n2t2)l56120其中：B为纬度，1uL-Lo,单位为弧度，N=,为卯酉圈曲率半/Vl-e-siirBI2-2

2、径，t=taiiB,2=e-cos2B,e=为第二偏心率，a为旋转椭球长半轴，bb为短半轴，X为子午线弧长。根据上式创建以GaussiaiiMapDirect命名的函数，函数输入输出格式为x,y,L0=GaussianMapDirect(B,L,RefEllipsoid)或者x,y.L0=GaussianMapDirect(B,L,a,f)RefEllipsoid为椭球参数RefEllipsoid=a,b,c,f,e2,e2；WGS84椭球参数：长半轴a=6378137扁率f=1/298.257223563b=a(l-Qc=a*a/b;e2=(a*a-b*b)/(a*a);e2_=(a*ab*

3、b)/(b*b);GaussianMapDiiect函数如下：fimctionx,y,L0=GaussianMapDirect(B,L,X)%WGS84椭球参数：a=6378137;%长半轴f=1/298.257223563;%扁率b=a*(l-f);%短半轴e=(sqrt(a八22)/玄;第一偏心率e_=(sqrt(a八2bT)/b;%第二偏心率%当地中央子午线决定于当地的直角坐标系统，首先确定您的直角坐标系统是3度带还是6度带投影，然后再根据如下公式推算。Q=input(iW选择投影带:n);ifQ=6%6度带：M=round(L+3)./6);%带号M=roimd(L+3)/6,即对(L

4、+3)/6的值四舍五入取整数，L为当地经度；L0=6.*M3;%则中央子午线经度L0=6XM-3else%3度带：M=round(L./3);%带号M=round(L/3),即对(L/3)的值四舍五入取整数，L为当地经度；L0=3.*M;%则中央子午线经度L0=3XMendl=(L-L0).*pi/180;N=a./(sqrt(l八2);t=tan(B);p=e_.*cos(B);%p表示yita%计算高斯平面坐标(x,y)x=X+(N.*(siii(B).*(cos(B).*(l.A2)./2+(N.*(siii(B).*(cos(B).A3).*(5-(t).A2)+9.*(p.八2)+4

5、.*(p.八4).*(1.八4)./24.+(N.*sin(B).*(cos(B).A5.*(61-58.*(t.A2)+(tA4)+270.*(p.A2)-330.*(p.A2).*(tA2).*(l.八6)./720;y=N.*cos(B).*l+(N.*(cos(B).A3.*(l-t.A2+p.A2).*(l.A3)./6+(N.*(cos(B).A5).*(5-18.(t.A2)+.t.八4+14严9.八2)5&(p.A2).*(t.A2).*(l.A5)./120;LO=degree3dms(LO);endlatitude2meridiaii函数如下：fimctionx=latit

6、iide2meiidian(B,a,e)%由纬度B求对应的子午线弧长x,计算公式mO=a*(l-e八2);ni2=(3*(e/x2)3|:m0)/2;m4=(5*(2)*1112)/4;m6=(7*(“2)*m4)/6;m8=(9*(“2)*m6)/8;a0=m0+m2/2+3*m4/8+5!*m6/16+35518/128;a2=m2/2+in4/2+15*m6/32+7*m8/16;a4=m4/8+3*m6/l6+7*1118/32;a6=m6/32+m8/16;a8=m8/128;x=aO*B-(a2*sin(2*B)/2+(a4*sin(4*B)/4-(a6*sin(6*B)/6+(a

7、8*sin(8*B)/8;end度分秒转度函数如下：fimctiondegiee=dins2degree(jiaodii)%度分秒(dd.mmss)度-4-4-degree=fix(jiaochi);niimute=fix(jiaodii-degiee)*100);second=(jiaodii-degiee-niimute/100)*10000;clegiee=degiee+mimute/60+second/3600;end度转度分秒函数如下:fimctiondms=degree3dms(du)%度一度分秒(dclmmss)degiee=fix(dii);niin=fix(dii-degie

8、e)*60);second=(dii-degiee)H：60-niin)*60);dins=degree+miii/100+seconc1/10000;end度分秒转弧度dms2rad函数如下:fimctionracl=clins2rad(n)cleg=fix(n);%度取所给数n的整数部分min_tem=(n-deg)*100;%去掉整数部分后小数点后移两位niin=fix(min_tem);%分取整sec=(niin_temmin)*100;%秒是小数点再向后移两位的数字-4-4-4-4-rad=(deg+miii/60.00+sec/3600)*pi/180.0;end2、高斯反算公式己

9、知高斯平面坐标（x,y）及指定中央子午线经度，计算大地坐标（B,L）：式+，Nf=九ZB”，Mf=a(1_列(1*sin也尸=CCS2Bftf=taiiBf,Bf为根据子午线弧长x反算的底点纬度。创建以GaussianMaphiveise命名的函数，函数输入输出格式为B丄=Gaussia11Mapinverse(x,y,RefEllipsoid)或者B,L=GaussianMaphiverse(x,y,a,f)GaussianMapInverse函数如下：fimctionL,B=GaussianMapInverse(x,y,LO)%WGS84椭球参数：a=6378137;%长半轴f=1/298

10、.257223563;%扁率b=a*(l-f);%短半轴e=(sqrt(a八22)/玄;第一偏心率e_=(s(pt(a八2-bA2)/b;%第二偏心率%根据中央子午线弧长x反算底点纬度BfBf=meridian21atitude(x,a,e);N4a./sqrt(l-(e.A2).*(sin(Bf).A2);Mf=a.*(l-e.A2)./scpt(l-(e.A2).*(sin(Bf).A2).A3);pf=e_.*cos(Bf);tf=tan(Bf);%己知高斯平面坐标(x,y)及指定中央子午线经度LO,计算大地坐标(E)B=Bf-tf*(y.A2)./(2.*Mf.*Nf)+tf*(5+3

11、.*(tf.A2)+pfA2-9.*(tfA2).*(pf.A2).*(y.A4)./(2+tf.*(61+90.*(tf.A2)+45.*(tf.A4).*(y.A6)./(720.3*:Mf.*(Nf.A5);L=L04y./(Nf.*cos(Bf)-(l+2.*(tf.A2)+pfA2).*y.A3./(6.*(NfA3).*cos(Bf).+(5+2&*tfA2+24.*tfA44.*pfA2+8.*(tf.A2).*pfA2).*y.A5./(120.*NfA5.*cos(Bf);B=rad2dins(B);L=rad2dms(L);end子午线弧长反算函数meri(Iian21at

12、itii(Ie如下：functionB=meridian21atitude(x,a,e)mO=a*(l-e八2);ni2=(3*(e2)*1110)/2;m4=(5*(2)*1112)/4;m6=(7*(2)*1114)/6;m8=(9*(“2)*m6)/8;a0=m0+m2/2+3水1114/8+5*1116/16+35*1118/128;a2=m2/2+in4/2+15*m6/32+7*m8/16;a4=m4/8+3*m6/l6+7恤8/32;a6=m6/32+m8/16;a8=m8/128;%纬度B的计算B0=x./a0;%B的初始值wliile1F=-(a2.*sin(2.*B0)./

13、2+(a4.*sin(4.*B0)./4-(a6.*sn(6.*B0)./6+(a8.*sin(8.*B0)./8;B=(x-F)./a0;ifabs(B0-B)10.A-6break;enclabs(B0-B)B0=B;endend弧度转度分秒rad2dms函数如下：fimctionchns=rad2chns(azimuth)%弧度转度分秒dgiee_tem=azimuthH：l80/pi;dgree=Jfix(dgree_tem);niin_tem=(clgiee_teindgree)水60);niin=fix(niin_tem);second=(min_tem-niin)*60);dms

14、=dgree+niiii/100+seconcl/10000;end度分秒转弧度dmsliad函数如下：%度取所给数ii的整数部分%去掉整数部分后小数点后移两位%分取整%秒是小数点再向后移两fimctionracl=clins2rad(n)cleg=fix(n);niin_tem=(n-deg)*100;min=fix(min_tem);sec=(niin_tem*100;位的数字racl=(deg+miii/60.00+sec/3600)H：pi/180.0;end3、实例计算验证首先读取文件datal.txt中的数据，计算其在相应六度带高斯投影带内的高斯平面直角坐标，并存贮在文件data2

15、.txt中，datal.txt格式为：经度(dclmmss)纬度(dd.mmss)大地高data2.txt格式为：x(m)y(m)中央子午线经度(dd.mmss)test61data文件程序如下：clear;%文件查找窗口%合并路径文件名clc;filename,patlmame=uigetfile(r*.*r);file=fiillfile(pathname,filename);77A=importclata(file);%将生成的文件导入工作空间，变量名为A%RefEllipsoicl为椭球参数%RefEllipsoid=a,b,c,f,e2,e2_;a=6378137;%WGS84椭球参

16、数：长半轴41/298.257223563;%扁率b=a*(l-f);%WGS84椭球参数：短半轴e=(sqrt(aA22)/玄；X=latitude2meiidian(dins2rad(A.data(:,2),a,e);%X为子午线弧长,有纬度B算岀x,yJLO=GaussiaiiMapDirect(chns2rad(A.data(:,2),dins2degiee(A.clata(:,l),X);B=x,yJLO;%B为重组矩阵fileiiameout,pathname_out=uiputfile(*.txt7请输入文件名*);%文件查找窗口fileout=fiillfile(pathnam

17、e_out,filenameout);%合并路径文件名ficfopeifileoutvt);%新建打开txt文件jprintf(fid,x(m)y(m)中央子午线经度(ddmmss)n);fclose(fid);运行结果显示计算datal.txt中的数据在相应六度带高斯投影带内的高斯平面直角坐标如下:文件(F)佔E)梧弍(0)egfVlW(H)x(m)2468078.0719291149369.3306908956597.5191818160632.482358v(m)-8487&594330166754.858518-10905.223119-26304.060915中央子午线经度(dd.m

18、mss)21.000000177.000000111.00000099.000000data35记事本data35记事本99-8-在相应三度带高斯投影带内的高斯平面直角坐标如下:data35记事本data35记事本99-8-data35记事本data35记事本99-8-文件(F)歸(E)磅云(0)牛辛助IHIx(m)2468078.0719291149345.9267658956597.5191818160632.482358y(m)中央子午线经度(dd.mmss)-84878.59433021.000000-161799.018071180.000000-10905.223119111.000000-26304.06091599.000000然后根据文件data2.txt中的高斯平面直角坐标及其中央子午线经度，计算其经纬度，并将计算结果按照格式存贮在文件data3.txt中，data3.txt格式为：经度(dclmmss)纬度(dd.mmss)test