matlab大地測量高斯投影正反算程序設(shè)計(jì)實(shí)驗(yàn)

1、-1-高斯投影正反算一、實(shí)驗(yàn)?zāi)康木帉懜咚雇队罢此愕某绦?，并對格式化文件?shù)據(jù)進(jìn)行計(jì)算，驗(yàn)證程序。二、實(shí)驗(yàn)內(nèi)容：1、高斯投影正算公式己知大地坐標(biāo)(B,L)及中央子午線經(jīng)度，計(jì)算高斯平面坐標(biāo)(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為旋轉(zhuǎn)橢球長半軸，bb為短半軸，X為子午線弧長。根據(jù)上式創(chuàng)建以GaussiaiiMapDirect命名的函數(shù)，函數(shù)輸入輸出格式為x,y,L0=GaussianMapDirect(B,L,RefEllipsoid)或者x,y.L0=GaussianMapDirect(B,L,a,f)RefEllipsoid為橢球參數(shù)RefEllipsoid=a,b,c,f,e2,e2；WGS84橢球參數(shù)：長半軸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函數(shù)如下：fimctionx,y,L0=GaussianMapDirect(B,L,X)%WGS84橢球參數(shù)：a=6378137;%長半軸f=1/298.257223563;%扁率b=a*(l-f);%短半軸e=(sqrt(a八22)/玄;第一偏心率e_=(sqrt(a八2bT)/b;%第二偏心率%當(dāng)?shù)刂醒胱游缇€決定于當(dāng)?shù)氐闹苯亲鴺?biāo)系統(tǒng)，首先確定您的直角坐標(biāo)系統(tǒng)是3度帶還是6度帶投影，然后再根據(jù)如下公式推算。Q=input(iW選擇投影帶:n);ifQ=6%6度帶：M=round(L+3)./6);%帶號M=roimd(L+3)/6,即對(L

4、+3)/6的值四舍五入取整數(shù)，L為當(dāng)?shù)亟?jīng)度；L0=6.*M3;%則中央子午線經(jīng)度L0=6XM-3else%3度帶：M=round(L./3);%帶號M=round(L/3),即對(L/3)的值四舍五入取整數(shù)，L為當(dāng)?shù)亟?jīng)度；L0=3.*M;%則中央子午線經(jīng)度L0=3XMendl=(L-L0).*pi/180;N=a./(sqrt(l八2);t=tan(B);p=e_.*cos(B);%p表示yita%計(jì)算高斯平面坐標(biāo)(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嚴(yán)9.八2)5&(p.A2).*(t.A2).*(l.A5)./120;LO=degree3dms(LO);endlatitude2meridiaii函數(shù)如下：fimctionx=latit

6、iide2meiidian(B,a,e)%由緯度B求對應(yīng)的子午線弧長x,計(jì)算公式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度分秒轉(zhuǎn)度函數(shù)如下：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度轉(zhuǎn)度分秒函數(shù)如下: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度分秒轉(zhuǎn)弧度dms2rad函數(shù)如下:fimctionracl=clins2rad(n)cleg=fix(n);%度取所給數(shù)n的整數(shù)部分min_tem=(n-deg)*100;%去掉整數(shù)部分后小數(shù)點(diǎn)后移兩位niin=fix(min_tem);%分取整sec=(niin_temmin)*100;%秒是小數(shù)點(diǎn)再向后移兩位的數(shù)字-4-4-4-4-rad=(deg+miii/60.00+sec/3600)*pi/180.0;end2、高斯反算公式己

9、知高斯平面坐標(biāo)（x,y）及指定中央子午線經(jīng)度，計(jì)算大地坐標(biāo)（B,L）：式+，Nf=九ZB”，Mf=a(1_列(1*sin也尸=CCS2Bftf=taiiBf,Bf為根據(jù)子午線弧長x反算的底點(diǎn)緯度。創(chuàng)建以GaussianMaphiveise命名的函數(shù)，函數(shù)輸入輸出格式為B丄=Gaussia11Mapinverse(x,y,RefEllipsoid)或者B,L=GaussianMaphiverse(x,y,a,f)GaussianMapInverse函數(shù)如下：fimctionL,B=GaussianMapInverse(x,y,LO)%WGS84橢球參數(shù)：a=6378137;%長半軸f=1/298

10、.257223563;%扁率b=a*(l-f);%短半軸e=(sqrt(a八22)/玄;第一偏心率e_=(s(pt(a八2-bA2)/b;%第二偏心率%根據(jù)中央子午線弧長x反算底點(diǎn)緯度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);%己知高斯平面坐標(biāo)(x,y)及指定中央子午線經(jīng)度LO,計(jì)算大地坐標(biāo)(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子午線弧長反算函數(shù)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的計(jì)算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弧度轉(zhuǎn)度分秒rad2dms函數(shù)如下：fimctionchns=rad2chns(azimuth)%弧度轉(zhuǎn)度分秒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度分秒轉(zhuǎn)弧度dmsliad函數(shù)如下：%度取所給數(shù)ii的整數(shù)部分%去掉整數(shù)部分后小數(shù)點(diǎn)后移兩位%分取整%秒是小數(shù)點(diǎn)再向后移兩fimctionracl=clins2rad(n)cleg=fix(n);niin_tem=(n-deg)*100;min=fix(min_tem);sec=(niin_tem*100;位的數(shù)字racl=(deg+miii/60.00+sec/3600)H：pi/180.0;end3、實(shí)例計(jì)算驗(yàn)證首先讀取文件datal.txt中的數(shù)據(jù)，計(jì)算其在相應(yīng)六度帶高斯投影帶內(nèi)的高斯平面直角坐標(biāo)，并存貯在文件data2

15、.txt中，datal.txt格式為：經(jīng)度(dclmmss)緯度(dd.mmss)大地高data2.txt格式為：x(m)y(m)中央子午線經(jīng)度(dd.mmss)test61data文件程序如下：clear;%文件查找窗口%合并路徑文件名clc;filename,patlmame=uigetfile(r*.*r);file=fiillfile(pathname,filename);77A=importclata(file);%將生成的文件導(dǎo)入工作空間，變量名為A%RefEllipsoicl為橢球參數(shù)%RefEllipsoid=a,b,c,f,e2,e2_;a=6378137;%WGS84橢球參

16、數(shù)：長半軸41/298.257223563;%扁率b=a*(l-f);%WGS84橢球參數(shù)：短半軸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)中央子午線經(jīng)度(ddmmss)n);fclose(fid);運(yùn)行結(jié)果顯示計(jì)算datal.txt中的數(shù)據(jù)在相應(yīng)六度帶高斯投影帶內(nèi)的高斯平面直角坐標(biāo)如下:文件(F)佔(zhàn)E)梧弍(0)egfVlW(H)x(m)2468078.0719291149369.3306908956597.5191818160632.482358v(m)-8487&594330166754.858518-10905.223119-26304.060915中央子午線經(jīng)度(dd.m

18、mss)21.000000177.000000111.00000099.000000data35記事本data35記事本99-8-在相應(yīng)三度帶高斯投影帶內(nèi)的高斯平面直角坐標(biāo)如下:data35記事本data35記事本99-8-data35記事本data35記事本99-8-文件(F)歸(E)磅云(0)牛辛助IHIx(m)2468078.0719291149345.9267658956597.5191818160632.482358y(m)中央子午線經(jīng)度(dd.mmss)-84878.59433021.000000-161799.018071180.000000-10905.223119111.000000-26304.06091599.000000然后根據(jù)文件data2.txt中的高斯平面直角坐標(biāo)及其中央子午線經(jīng)度，計(jì)算其經(jīng)緯度，并將計(jì)算結(jié)果按照格式存貯在文件data3.txt中，data3.txt格式為：經(jīng)度(dclmmss)緯度(dd.mmss)test