MATLAB批量读取航摄相片EXIF信息和GNSS信息以及MATLAB批量经纬度坐标转换空间直角坐标

前言

没有前言，创作不易，希望家人们点点赞，如有错误评论留言，感谢支持！

---

---

一、MATLAB批量读取航摄相片的EXIF信息，并批量提取其GNSS经纬度信息

（1）前述：

所谓航摄相片的EXIF信息，就是 （可交换图像文件格式）（英语：Exchangeable image file format，官方简称Exif），是专门为数码相机的照片设定的，可以记录数码照片的属性信息和拍摄数据。下图分别是航摄图片的详细属性信息以及在matlab中读取的结果。

图1 图片属性信息

图2 MATLAB显示

GNSS经纬度信息就是BLH（纬度，经度，大地高），我们在图片的详细信息可以看见，而在matlab中读取之后字段名分别为BLH（GPSLatitude，GPSLongitude，GPSAltitude）

（2）代码实现和结果：

1: 代码思路： （批量的航摄相片===>matlab读取每张相片EXIF信息===>读取EXIF信息的GNSS经纬度坐标===>经纬度的度分秒转换为度===>导出保存为xls文件。）
2：程序清单：
(1) dir函数----格式：a=dir(‘filename’ , ‘*.txt / *.jpg / *.tiff’ ) 其表示将盘定位到你存放文件的文件里面，会列举出你指定问filename路径下的固定格式的文件，我下面的代码是由于我将相片文件放在和程序运行保存的文件夹里面，所以没有fliename
(2) imfinfo函数----格式：a=imfinfo(‘filename’) 其用于读取图片的文件信息，相同用法的还有exifread函数
(3) xlswrite函数----格式：xlswrite(’ filename+名字.xls ’ , 变量名) 用于将变量里面的内容保存为xls表格输出

clc;
clear all;
data=dir('*.JPG');
m=[];
n=[];
%data.name表示相片的全名，用m集合储存; data.folder表示相片所在文件夹的位置，用n集合储存
for i=1:length(data)
    m_data=data(i).name;
    folder=data(i).folder;
    m=[m;m_data];
    n=[n;folder];
end
%利用imfinfo函数循环读取每一张相片的exif信息
GNSS=[];
for j=1:length(data)
    GNSS_data=imfinfo([n(j,:),'\',m(j,:)]);
    %这里的imfinfo('filename'),这里的filename就是第j张相片的文件夹的位置+   \   +第j张相片的全名
    GNSS=[GNSS;GNSS_data];
end
%将exif信息中的GPS信息赋值给GNSS_data,依次读取其经纬度坐标
E=[];
N=[];
H=[];
for k=1:length(data)
   GNSS_data=GNSS(k,:).GPSInfo;
   E_data=GNSS_data.GPSLatitude;
   N_data=GNSS_data.GPSLongitude;
   H_data=GNSS_data.GPSAltitude;
   E=[E;E_data];
   N=[N;N_data];
   H=[H;H_data];
end
%注意由于相片经纬度坐标的组织形式，因此读取的经纬度的度分秒以矩阵形式存在
%将度分秒转换为度为单位
EAST=E(:,1)+E(:,2)/60+E(:,3)/3600;
NORTH=N(:,1)+N(:,2)/60+N(:,3)/3600;
image_GNSS=[EAST,NORTH,H];
%保存转换后的坐标为xls
xlswrite('C:\Users\稳魂\Desktop\GNSS.xls',image_GNSS);

我们将程序提取保存后的经纬度的表格打开，可以看见明显的BLH格式。

二、地理经纬度坐标转空间直角坐标原理及公式

（1） 转换原理

我国的参心坐标系可以分为参心大地坐标系和参心空间直角坐标系；参心大地坐标系的坐标元素是BLH（也就是我们常见的经纬度坐标）而参心空间直角坐标系的坐标元素是XYZ（注意不是笛卡尔坐标系，xy的方向是反的），他们的原点都是参考椭球的原点。

图1 参心大地坐标系

图2 参心空间直角坐标系

而其转换原理就是将两个坐标系上的同名点投影到点所在的子午平面上，建立对应的数学关系，进行转换。

（2） 转换公式

我们设A点的大地坐标（B，L，H）以及对应的空间直角坐标（X，Y，Z）以及下面的一系列参数：
椭球基准：西安80坐标系
参考椭球长半轴：a=6378140±5（m）
参考椭球短半轴：b=6356755.2882（m）
过A点的卯酉圈曲率半径：

参考椭球的第二偏心率：

公式：

三、利用MATLAB批量转换坐标

2：程序清单：

(1) importdata函数—格式：a=imortdata(‘filename’) 用于读取表格以及文本格式文件
(2) format long g：用于关闭数据科学计数法显示

clc;
clear all;
%利用importdata函数读取上一个程序保存的GNSS经纬度表格
BLH_data=importdata('C:\Users\稳魂\Desktop\GNSS.xls');
B=BLH_data(:,1);
L=BLH_data(:,2);
H=BLH_data(:,3);
a=6378140;
b=6356755.2882;
e2=(a^2-b^2)/(a^2);
NN=[];
format long g;
for i=1:length(BLH_data)
    N=a/((1-e2*(sind(B(i,:)).^2)).^0.5);
    NN=[NN;N];
end
XX=[];
YY=[];
ZZ=[];
for j=1:length(BLH_data)
    X=((NN(j,:)+H(j,:)).*cosd(B(j,:)).*cosd(L(j,:)));
    Y=((NN(j,:)+H(j,:)).*cosd(B(j,:)).*sind(L(j,:)));
    Z=((NN(j,:).*(1-e2)+H(j,:)).*sind(B(j,:)));
    XX=[XX;X];
    YY=[YY;Y];
    ZZ=[ZZ;Z];
end
XYZ=[XX,YY,ZZ];
figure;
xlswrite('C:\Users\稳魂\Desktop\XYZ.xls',XYZ);

四、根据外部工具验证其转换精度

根据上述转换后的坐标数据，绘制其这一组航摄相片的航线图

plot(XX,YY,'--gs',...
    'LineWidth',2,...
    'MarkerSize',3,...
    'MarkerEdgeColor','r',...
    'MarkerFaceColor',[1,0.5,0.5]);
xlim([-1517450 -1517150]);
ylim([5268200 5268340]);
xlabel('Y',"Color",'r');
ylabel('X',"Color",'r');
title('该航摄相片组的航线','color','r');
grid on；

(1)利用pix4D的航线报告与MATLAB转换后绘制的航线比较

需要注意的是我在pix4D中选用的WGS 84 UTM 48N的平面投影，是经纬度与平面坐标之间的转换
而我计算的经纬度与空间直角坐标之间的转换，所以航线会有差别，这里只为了做参考。

图1 MATLAB绘制的航线

图2 pix4D质量报告的航线

(2)利用坐标转换COORD软件转换后与MATLAB转换后的结果比较

我们抽取第一张相片的经纬度利用坐标转换软件转换后与matlab转换后的结果比较，转换精度还可以，至于毫米级之后的精度，是因为西安80椭球参数长半轴a=6378140±5（m）,我直接省略了±5，所以各位可以自行选用。

图1 COORD软件

图2 MATLAB