【月球Ti/Fe元素反演-基于MATLAB（一）】基于Matlab，对月球上的多光谱（IMG数据）进行Fe与Ti元素的反演

主要步骤

1. 分段读取数据 ，主要基于multibandread函数；
2. 基于反演公式 进行元素反演；
3. 将超大矩阵保存为Tif格式；

代码位置：C:\Users\LEI\Documents\MATLAB\TiO2反演\read_mutiband\分幅\LQ11_FINAL

---

一、分段读取数据（multibandread）

1. 提取hdr文件关键词
通过hdr文件，提取出samples，lines，bands，precision，offset，interleave，byteorder等因子。
注： IMG里面的图像像素数据都需要再乘一个变形因子scaling_factor = 2e-5; 由于多光谱里面的数据都是反射率（即，0-1之间），所以直接读取的数据是几千左右的，是变形过的。

2. multibandread

// 不规定range读数
R_b1_b2 = double(multibandread(filename,[2 samples 1], precision, offset, interleave, byteorder));

// 分别规定了行，列，波段的range范围
Row_450 = double(multibandread(filename,[lines samples bands], precision, offset, interleave, byteorder,{'Row', 'Range', [A(1) A(end)]},{'Column','Range',[1,samples]},{'Band', 'Range', [1 1]}));
Row_750 = double(multibandread(filename,[lines samples bands], precision, offset, interleave, byteorder,{'Row', 'Range', [A(1) A(end)]},{'Column','Range',[1,samples]},{'Band', 'Range', [2 2]}));
Row_950 = double(multibandread(filename,[lines samples bands], precision, offset, interleave, byteorder,{'Row', 'Range', [A(1) A(end)]},{'Column','Range',[1,samples]},{'Band', 'Range', [4 4]}));

运行结果：这里R_b1_b2 是直接按[2 samples 1]这种格式读取IMG数据。又因为是‘bil’的存储方式（即，按band1的低第一行存储，接着是band2的第一行； 所以如果multiband不规定range的话，直接读取的话，默认6326那一行就是band2的第一行），所以R_b1_b2的第一行和Row_450 对应，第二行对应Row_750 ，依次类推。第9+1行对应第一波段（Row_450 ）的第二行。

二、反演公式

1. 反演的时候因为得到的是[0,1]区间的百分比含量，是double类型的。可以统一乘以100（factor），再转化为uint16类型以节省空间；

% 由于原图的反射率在[0,1]区间内像为了方便保存，统一放大了。所以在计算时，需要全部乘scaling_factor ；
% 才能得到反射率的真实值
scaling_factor = 2e-5;
theta_Ti = atan(((Row_450./Row_750)-0.208)./(Row_750*scaling_factor+0.108));
theta_Fe = -atan(((Row_950./Row_750)-1.25)./(Row_750*scaling_factor-0.037));
factor = 100；
% Fe的取值范围[0~10（%）]，保存时乘以了一个factor=100；即像素范围为uint16[0~1000（%）]
Wt_Ti = real(0.72*(theta_Ti.^14.964));
Wt_Ti(Wt_Ti<0)=0;
Wt_Ti(Wt_Ti>10)=10;
Wt_Ti = uint16(Wt_Ti*factor);

% Ti的取值范围[0~20（%）]，保存时乘以了一个factor=100；即像素范围为uint16[0~20000（%）]
Wt_Fe = real(20.527*theta_Fe-12.266);
Wt_Fe(Wt_Fe<0)=0;
Wt_Fe(Wt_Fe>20)=20;
Wt_Fe = uint16(Wt_Fe*factor);

三、将MAT数据，保存Tiff格式（geotiffwrite）

1. 首先每128行保存一次分别3个波段的数据，利用反演公式后：动态保存矩阵为MAT格式：

save('Ti_LQ11_uint16.mat','Wt_Ti','-v7.3');
temp_Ti = matfile('Ti_LQ11_uint16mat','Writable',true);
[nrows, ncols] = size(temp_Ti,'Wt_Ti');
if nrows==0
	temp_Ti.Wt_Ti=real(Wt_Ti);
else
    temp_Ti.Wt_Ti(nrows+1:nrows+128,:)=real(Wt_Ti);
end

2. 将这个temp_Ti.Wt_Ti（>10GB）矩阵MAT格式，write为tiff格式（geotiffwrite）：
由于单纯想保存超大tiff格式，这里的R（spatial reference）随便设置了，具体可以参加改为博主的介绍：
geotiffwrite详细介绍

clc;clear all;
R=[0,-1;1,0;0,6000];
load('Ti_LQ11_uint16.mat')
geotiffwrite('LQ11_Ti_inv.tif', Wt_Ti, R, 'CoordRefSysCode', 21417,'TiffType','bigtiff');