【月球Ti/Fe元素反演

月球Ti/Fe元素反演"/>

【月球Ti/Fe元素反演

一. 裁剪IMG格式的图像：参考链接1

1.首先打开要裁剪的图像：file->open image file
2.找到ROI工具，在快捷栏中直接找到roi工具，单击；或者通过File—new—Region of Interest。
3.描出感兴趣区域，在Region of interest (ROI) Tool面板中，选择File-> Save as，保存绘制的多边形ROI，选择保存的路径和文件名。


取名为.xml的后缀
4. 在右边Toolbox搜索框中搜subset,下面会出现SubsetData from ROIs,双击，在出现的框中选择你要裁剪的文件名，点击OK。

5.在Subset Data from ROIs Parameters面板中，设置以下参数：
Select Input ROIs：选择刚才生成的矢量文件roi，选择输出路径和文件名，单击OK执行图像裁剪。
保存为.IMG格式，并自动生成头文件.hdr

二. 将IMG格式文件转化为TIF：参考链接2

1. 在TOC面板，选择导出：
   
2. 并输入保存的TIF文件名字：
   

三. 基于MATLAB赋值地理信息

1. 分别读取.img和.hdr文件的像素值和行列波段参数等:
   注: 精度对应的类型如下
   case 1：
   precision=‘uint8=>uint8’;%头文件中datatype=1对应ENVI中数据类型为Byte，对应MATLAB中数据类型为uint8
   case 2：
   precision=‘int16=>int16’;%头文件中datatype=2对应ENVI中数据类型为Integer，对应MATLAB中数据类型为int16
   case 12：
   precision=‘uint16=>uint16’;%头文件中datatype=12对应ENVI中数据类型为Unsighed Int，对应MATLAB中数据类型为uint16
   case 3：
   precision=‘int32=>int32’;%头文件中datatype=3对应ENVI中数据类型为Long Integer，对应MATLAB中数据类型为int32
   case 13：
   % precision=‘uint32=>uint32’;%头文件中datatype=13对应ENVI中数据类型为Unsighed Long，对应MATLAB中数据类型为uint32
   case 4：
   % precision=‘float32=>float32’;%头文件中datatype=4对应ENVI中数据类型为Floating Point，对应MATLAB中数据类型为float32
   case 5：
   % precision=‘double=>double’;%头文件中datatype=5对应ENVI中数据类型为Double Precision，对应MATLAB中数据类型为double
   otherwise：
   % error(‘invalid datatype’);%除以上几种常见数据类型之外的数据类型视为无效的数据类型
   % end

% 文件地址：C:\Users\LEI\Documents\MATLAB\TiO2反演\read_mutiband\分幅\LQ11_90分之1
file = ['Z:\Data\Moon_data\KAGUYA\MI\img2tif\l90_1ROI_LQ11'];
filename = [file,'.img'];
fid = fopen([file,'.hdr'],'r');
% fid = fopen('Z:\Data\Moon_data\KAGUYA\MI\MI_Multi_band_Imager\MI_MAP_03_N00E160S01E161SC\MI_MAP_03_N00E160S01E161SC.ctg','r');
info = fread(fid,'char=>char');
info = info'%% 并提取.hdr如下参数：
a=strfind(info,'samples = ');
b=length('samples = ');
c=strfind(info,'lines');
samples=[];
for i=a+b:c-1samples=[samples,info(i)];
end
samples=str2num(samples);
lines   = 51200;
bands   = 9;
precision=['int16=>int16']; %因为data type = 2
offset = 0;
interleave = 'bil';
byteorder = 'ieee-le';

2. 反演-计算Fe，Ti含量

%% for 全部,记得每次改save的名字！！
% !!!!!记得最开始把Fe_LQ11里面的Wt_Fe=[]然后保存。
S = [1:lines];
step = 2; %两行两行的读取 
s = reshape(S,step,[]);
s1 = s(:,2);
Wt_Ti=[];
Wt_Fe = [];
save('Ti_LQ11_uint16.mat','Wt_Ti','-v7.3');% 先保存两个空的MAT变量
save('Fe_LQ11_uint16.mat','Wt_Fe','-v7.3');% 先保存两个空的MAT变量
factor=100; %统一将计算的含量放大100倍保存；
for i=1:length(s)A =  s(:,i);AA = [A(1),A(end)];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]}));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));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);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);temp_Ti = matfile('Ti_LQ11_uint16.mat','Writable',true);temp_Fe = matfile('Fe_LQ11_uint16.mat','Writable',true);[nrows, ncols] = size(temp_Ti,'Wt_Ti');if nrows==0temp_Ti.Wt_Ti=real(Wt_Ti);temp_Fe.Wt_Fe=real(Wt_Fe);elsetemp_Ti.Wt_Ti(nrows+1:nrows+step,:)=real(Wt_Ti);temp_Fe.Wt_Fe(nrows+1:nrows+step,:)=real(Wt_Fe);enddisp(['第',num2str(i),'次'])
end

3. 将Ti含量的矩阵保存为tif文件，并附上地理坐标：

[A, R] = geotiffread('Z:\Data\Moon_data\KAGUYA\MI\img2tif\LQ11_90per_arcgis.TIF');
info1 = geotiffinfo('Z:\Data\Moon_data\KAGUYA\MI\img2tif\LQ11_90per_arcgis.TIF');
geoTags = info1.GeoTIFFTags.GeoKeyDirectoryTag;
% load('Ti_LQ11_uint16.mat') or temp_Ti.Wt_Ti
geotiffwrite('Z:\Data\Moon_data\KAGUYA\MI\Ti_Fe反演结果\LQ11_Ti_inv_90per_geowrite.tif', temp_Ti.Wt_Ti, R, 'GeoKeyDirectoryTag', geoTags ,'TiffType','bigtiff');

注意：
这里先读取第二步另存为的tif文件的R和geoTags；
再用geotiffwrite函数，将R和geoTags赋值给上一步计算得到的temp_Ti.Wt_Ti；
如果直接用geotiffwrite(filename,A,R)；则会报错：参考链接3

R 如下：

info1 如下：

geoTags 如下：