matlab大地测量高斯投影正反算程序设计实验.pdf
- 文档编号:3431604
- 上传时间:2023-05-05
- 格式:PDF
- 页数:10
- 大小:194.03KB
matlab大地测量高斯投影正反算程序设计实验.pdf
《matlab大地测量高斯投影正反算程序设计实验.pdf》由会员分享,可在线阅读,更多相关《matlab大地测量高斯投影正反算程序设计实验.pdf(10页珍藏版)》请在冰点文库上搜索。
高斯投影正反算高斯投影正反算一、实验目的一、实验目的编写高斯投影正反算的程序,并对格式化文件数据进行计算,验证程序。
二、实验内容:
二、实验内容:
1、高斯投影正算公式高斯投影正算公式已知大地坐标(B,L)及中央子午线经度L0,计算高斯平面坐标(x,y),公式如下:
6222425442232)3302705861)(cos)sin(720)495)(cos)sin(24)cos()sin(2ltttBBNltBBNlBBNXx52224253223)5814185)(cos120)1)(cos6)cos(ltttBNltBNlBNy其中:
B为纬度,0LLl,单位为弧度,BeaN22sin1,为卯酉圈曲率半径,Bttan,Be222cos,bbae22为第二偏心率,a为旋转椭球长半轴,b为短半轴,X为子午线弧长。
根据上式创建以GaussianMapDirect命名的函数,函数输入输出格式为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(1-f)c=a*a/b;e2=(a*a-b*b)/(a*a);-1-e2_=(a*a-b*b)/(b*b);GaussianMapDirect函数函数如下如下:
functionx,y,L0=GaussianMapDirect(B,L,X)%WGS84椭球参数:
a=6378137;%长半轴f=1/298.257223563;%扁率b=a*(1-f);%短半轴e=(sqrt(a2-b2)/a;%第一偏心率e_=(sqrt(a2-b2)/b;%第二偏心率%当地中央子午线决定于当地的直角坐标系统,首先确定您的直角坐标系统是3度带还是6度带投影,然后再根据如下公式推算。
Q=input(请选择投影带:
n);ifQ=6%6度带:
M=round(L+3)./6);%带号M=round(L+3)/6,即对(L+3)/6的值四舍五入取整数,L为当地经度;L0=6.*M-3;%则中央子午线经度L0=6M-3else%3度带:
M=round(L./3);%带号M=round(L/3),即对(L/3)的值四舍五入取整数,L为当地经度;L0=3.*M;%则中央子午线经度L0=3Mendl=(L-L0).*pi/180;N=a./(sqrt(1-(e2).*(sin(B).2);t=tan(B);p=e_.*cos(B);%p表示yita-2-%计算高斯平面坐标(x,y)x=X+(N.*(sin(B).*(cos(B).*(l.2)./2+(N.*(sin(B).*(cos(B).3).*(5-(t).2)+9.*(p.2)+4.*(p.4).*(l.4)./24.+(N.*sin(B).*(cos(B).5.*(61-58.*(t.2)+(t.4)+270.*(p.2)-330.*(p.2).*(t.2).*(l.6)./720;y=N.*cos(B).*l+(N.*(cos(B).3.*(1-t.2+p.2).*(l.3)./6+(N.*(cos(B).5).*(5-18.*(t.2)+.t.4+14.*(p.2)-58.*(p.2).*(t.2).*(l.5)./120;L0=degree3dms(L0);endlatitude2meridian函数如下:
函数如下:
functionx=latitude2meridian(B,a,e)%由纬度B求对应的子午线弧长x,计算公式m0=a*(1-e2);m2=(3*(e2)*m0)/2;m4=(5*(e2)*m2)/4;m6=(7*(e2)*m4)/6;m8=(9*(e2)*m6)/8;a0=m0+m2/2+3*m4/8+5*m6/16+35*m8/128;a2=m2/2+m4/2+15*m6/32+7*m8/16;a4=m4/8+3*m6/16+7*m8/32;a6=m6/32+m8/16;a8=m8/128;x=a0*B-(a2*sin(2*B)/2+(a4*sin(4*B)/4-(a6*sin(6*B)/6+(a8*sin(8*B)/8;end度分秒转度函数如下:
度分秒转度函数如下:
functiondegree=dms2degree(jiaodu)%度分秒(dd.mmss)度-3-degree=fix(jiaodu);mimute=fix(jiaodu-degree)*100);second=(jiaodu-degree-mimute/100)*10000;degree=degree+mimute/60+second/3600;end度转度分秒函数如下:
度转度分秒函数如下:
functiondms=degree3dms(du)%度度分秒(dd.mmss)degree=fix(du);min=fix(du-degree)*60);second=(du-degree)*60-min)*60);dms=degree+min/100+second/10000;end度分秒度分秒转转弧度弧度dms2rad函数如下函数如下:
functionrad=dms2rad(n)deg=fix(n);%度取所给数n的整数部分min_tem=(n-deg)*100;%去掉整数部分后小数点后移两位min=fix(min_tem);%分取整sec=(min_tem-min)*100;%秒是小数点再向后移两位的数字rad=(deg+min/60.00+sec/3600)*pi/180.0;end2、高斯反算公式高斯反算公式已知高斯平面坐标(x,y)及指定中央子午线经度L0,计算大地坐标(B,L):
64254222232)459061(720)935(242yttNMtyttNMtyNMtBBffffffffffffffff522242532230)8624285(cos1201)21(cos61cos1ytttBNytBNyBNLLfffffffffffff-4-式中,ffBeaN22sin1,3222)sin1()1(ffBeeaM,ffBe222cos,ffBttan,fB为根据子午线弧长x反算的底点纬度。
创建以GaussianMapInverse命名的函数,函数输入输出格式为B,L=GaussianMapInverse(x,y,RefEllipsoid)或者B,L=GaussianMapInverse(x,y,a,f)GaussianMapInverse函数如下:
函数如下:
functionL,B=GaussianMapInverse(x,y,L0)%WGS84椭球参数:
a=6378137;%长半轴f=1/298.257223563;%扁率b=a*(1-f);%短半轴e=(sqrt(a2-b2)/a;%第一偏心率e_=(sqrt(a2-b2)/b;%第二偏心率%根据中央子午线弧长x反算底点纬度BfBf=meridian2latitude(x,a,e);Nf=a./sqrt(1-(e.2).*(sin(Bf).2);Mf=a.*(1-e.2)./sqrt(1-(e.2).*(sin(Bf).2).3);pf=e_.*cos(Bf);tf=tan(Bf);%已知高斯平面坐标(x,y)及指定中央子午线经度L0,计算大地坐标(B,L)B=Bf-tf.*(y.2)./(2.*Mf.*Nf)+tf.*(5+3.*(tf.2)+pf.2-9.*(tf.2).*(pf.2).*(y.4)./(24.*Mf.*(Nf.3).+tf.*(61+90.*(tf.2)+45.*(tf.4).*(y.6)./(720.*Mf.*(Nf.5);L=L0+y./(Nf.*cos(Bf)-(1+2.*(tf.2)+pf.2).*y.3./(6.*(Nf.3).*cos(Bf).+(5+28.*tf.2+24.*tf.4+6.*pf.2+8.*(tf.2).*pf.2).*y.5./(120.*Nf.5.*cos(Bf);B=rad2dms(B);L=rad2dms(L);-5-end子午线弧长反算函数子午线弧长反算函数meridian2latitude如下:
如下:
functionB=meridian2latitude(x,a,e)m0=a*(1-e2);m2=(3*(e2)*m0)/2;m4=(5*(e2)*m2)/4;m6=(7*(e2)*m4)/6;m8=(9*(e2)*m6)/8;a0=m0+m2/2+3*m4/8+5*m6/16+35*m8/128;a2=m2/2+m4/2+15*m6/32+7*m8/16;a4=m4/8+3*m6/16+7*m8/32;a6=m6/32+m8/16;a8=m8/128;%纬度B的计算B0=x./a0;%B的初始值while1F=-(a2.*sin(2.*B0)./2+(a4.*sin(4.*B0)./4-(a6.*sin(6.*B0)./6+(a8.*sin(8.*B0)./8;B=(x-F)./a0;ifabs(B0-B)10.-6break;endabs(B0-B)B0=B;endend弧度转度分秒弧度转度分秒rad2dms函数如下函数如下:
functiondms=rad2dms(azimuth)%弧度转度分秒-6-dgree_tem=azimuth*180/pi;dgree=fix(dgree_tem);min_tem=(dgree_tem-dgree)*60);min=fix(min_tem);second=(min_tem-min)*60);dms=dgree+min/100+second/10000;end度分秒度分秒转转弧度弧度dms2rad函数如下函数如下:
functionrad=dms2rad(n)deg=fix(n);%度取所给数n的整数部分min_tem=(n-deg)*100;%去掉整数部分后小数点后移两位min=fix(min_tem);%分取整sec=(min_tem-min)*100;%秒是小数点再向后移两位的数字rad=(deg+min/60.00+sec/3600)*pi/180.0;end3、实例计算验证实例计算验证首先读取文件data1.txt中的数据,计算其在相应六度带高斯投影带内的高斯平面直角坐标,并存贮在文件data2.txt中,data1.txt格式为:
经度(dd.mmss)纬度(dd.mmss)大地高data2.txt格式为:
x(m)y(m)中央子午线经度(dd.mmss)test61data文件程序如下文件程序如下:
clear;clc;filename,pathname=uigetfile(*.*);%文件查找窗口file=fullfile(pathname,filename);%合并路径文件名-7-A=importdata(file);%将生成的文件导入工作空间,变量名为A%RefEllipsoid为椭球参数%RefEllipsoid=a,b,c,f,e2,e2_;a=6378137;%WGS84椭球参数:
长半轴f=1/298.257223563;%扁率b=a*(1-f);%WGS84椭球参数:
短半轴e=(sqrt(a2-b2)/a;X=latitude2meridian(dms2rad(A.data(:
2),a,e);%X为子午线弧长,有纬度B算出x,y,L0=GaussianMapDirect(dms2rad(A.data(:
2),dms2degree(A.data(:
1),X);B=x,y,L0;%B为重组矩阵filename_out,pathname_out=uiputfile(*.txt,请输入文件名);%文件查找窗口fileout=fullfile(pathname_out,filename_out);%合并路径文件名fid=fopen(fileout,wt);%新建打开txt文件fprintf(fid,x(m)y(m)中央子午线经度(dd.mmss)n);fprintf(fid,%ft%ft%ftn,B(:
:
);fclose(fid);运行结果显示运行结果显示计算计算data1.txt中的数据在相应六度带高斯投影带内的高斯平面直角坐标如下:
中的数据在相应六度带高斯投影带内的高斯平面直角坐标如下:
-8-在相应三在相应三度带高斯投影带内的高斯平面直角坐标如下:
度带高斯投影带内的高斯平面直角坐标如下:
然后根据文件data2.txt中的高斯平面直角坐标及其中央子午线经度,计算其经纬度,并将计算结果按照格式存贮在文件data3.txt中,data3.txt格式为:
经度(dd.mmss)纬度(dd.mmss)test62data文件程序如下:
文件程序如下:
clear;clc;filename,pathname=uigetfile(*.*);%文件查找窗口file=fullfile(pathname,filename);%合并路径文件名A=importdata(file);%将生成的文件导入工作空间,变量名为AL,B=GaussianMapInverse(A.data(:
1),A.data(:
2),dms2rad(A.data(:
3);C=L,B;%C为重组矩阵filename_out,pathname_out=uiputfile(*.txt,请输入文件名);%文件查找窗口fileout=fullfile(pathname_out,filename_out);%合并路径文件名fid=fopen(fileout,wt);%新建打开txt文件fprintf(fid,经度(dd.mmss)纬度(dd.mmss)n);fprintf(fid,%ft%ftn,C(:
:
);fclose(fid);运行结果显示运行结果显示将三度带或六度带的高斯平面坐标转化为大地坐标的结果如下将三度带或六度带的高斯平面坐标转化为大地坐标的结果如下:
-9-
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- matlab 大地测量 投影 正反 程序设计 实验
![提示](https://static.bingdoc.com/images/bang_tan.gif)