测绘—matlab—坐标正反算

在测量实验过程中需要进行单个、批量坐标正反算时可以用计算机运算提高效率

1.正算函数

输入的azimuth是弧度值

[x2,y2] = xy_direct(x1,y1,distance, azimuth)

x1、y1、x2、y2是坐标值，distance是两点间距离， azimuth是坐标方位角

function  [x2,y2] = xy_direct(x1,y1,distance, azimuth) %正算函数
x2 = x1+distance.*cos(azimuth);
y2 = y1+distance.*sin(azimuth);
end

1.1度分秒转弧度的函数

function degree = dms2rad(jiaodu)%度分秒（dd.mmss）－>弧度
 degree = fix(jiaodu);
 mimute = fix((jiaodu-degree)*100);
 second = (jiaodu-degree-mimute/100)*10000;
 degree = ((degree+mimute/60+second/3600).*pi)/180;
end

2.反算函数

function [distance,azimuth] = xy_inv(x1,y1,x2,y2)%反算函数，注意该过程三角函数算出来的是弧度
distance = sqrt((x2-x1).^2+(y2-y1).^2);%算出两点间的距离
if (x2-x1)>0 && (y2-y1)>0 %象限角在第一象限时
    azimuth = atan2(abs(y2-y1),abs(x2-x1));%方位角等于象限角
elseif (x2-x1)>0 && (y2-y1)<0%象限角在第四象限时
    azimuth = 2*pi-atan2(abs(y2-y1),abs(x2-x1));%方位角等于“360-象限角”
elseif (x2-x1)<0 && (y2-y1)>0%象限角在第二象限时
    azimuth = pi-atan2(abs(y2-y1),abs(x2-x1));%方位角等于“180-象限角”
elseif (x2-x1)<0 && (y2-y1)<0%象限角在第三象限时
    azimuth = pi+atan2(abs(y2-y1),abs(x2-x1));%方位角等于“180+象限角”
elseif (x2-x1) == 0 && (y2-y1)>0%y轴正方向
    azimuth=pi;
elseif (x2-x1) == 0 && (y2-y1)<0%y轴负方向
    azimuth=1.5*pi;
elseif (x2-x1) > 0 && (y2-y1)==0%x轴正方向
    azimuth=0;
elseif (x2-x1) < 0 && (y2-y1)==0%x轴负方向
    azimuth=pi;
else                            %其他情况
    azimuth=[];
end
end

2.1弧度转度分秒函数（将算出来的azimuth值转换为角度值）

function a = rad2dms(x)%弧度化度分秒函数
b = x*180/pi;%弧度化角度
c = fix(b);%取度数
d = (b-c)*60;%度的小数部分化为分
e = fix(d);%取分的整数
f = e/100;%分改为对应小数的数值
g = d-e;%取分的小数部分
h = g*60/10000;%分的小数部分化为秒
a = f+c+h;%将度分秒合在一起
end

3.总引用实验程序

[filename,pathname] = uigetfile('data1.txt');%找到文件data1.txt的文件名、文件路径。
file = [pathname,filename];%将文件名与文件路径合并成“路径\文件名”形式
A = importdata(file);%打开文件data1.txt中的数据
out1 = fopen('data2.txt','wt+');%新建文件夹data2.txt，并以可读可写的.txt格式打开
out2 = fopen('data3.txt','wt+');%新建文件夹data2.txt，并以可读可写的.txt格式打开
A.data(:,4) = dms2degree(A.data(:,4));%将第四行中“度分秒-->>>度”
[A.data(:,3),A.data(:,4)] = xy_direct(A.data(:,1),A.data(:,2),A.data(:,3), A.data(:,4));%将坐标x1,x2,距离，方位角带入正算函数[x2,y2] = xy_direct(x1,y1,distance, azimuth)进行计算并将算得的值赋值给A中的3、4列
fprintf(out1,'%s\n','  x1      y1     x2       y2');%在out1文件中写入表头
fprintf(out1,'%.2f %.2f %.2f %.2f\n',A.data(:,[1,2,3,4])');%将算后的矩阵A.data写入文件out1中
[m,~] = size(A.data);%找出A.data的行数m进行循环运算
for i = 1:m   %循环开始于1，到m结束
[A.data(i,5),A.data(i,6)] = xy_inv(A.data(i,1),A.data(i,2),A.data(i,3),A.data(i,4));%将A.data中的数据带入到反算函数[distance,azimuth] = xy_inv(x1,y1,x2,y2)中运算并将所得结果填入到A.data的5、6列
end
A.data(:,6) = rad2dms(A.data(:,6));%将函数算出的弧度角值转换为度分秒形式
fprintf(out2,'%s\n','    x1       y1      x2       y2     距离  方位角');%在out2文件中写入表头
fprintf(out2,'%.4f %.4f %.4f %.4f %.2f %.4f\n',A.data(:,[1,2,3,4,5,6])');%将算后的矩阵A.data写入文件out2中
fclose('all');%关闭所有打开的文件

3.1此次实验中的源数据文件data1.txt是下面形式。如若数据文件形式不同、需要改代码读取数据如：A.data(:,4)中的“数字”或可以改一下源数据形式也可。