基于matlab的FY-4（全圆盘）影像几何纠正

折腾了三个星期，终于搞定了FY-4的几何纠正问题。期间求人无数，总结起来：还是得靠自己。慢慢摸索，一个个去试，在试错的过程中也让我学会不少东西。话不多说，进入正题。

1.有关FY-4的介绍以及几何校正的两种方法，可以参考官网：http://satellite.nsmc.org.cn/PortalSite/StaticContent/DocumentDownload.aspx?TypeID=3

2.方法选择

         目前有两种方法可以实现FY-4的几何纠正，一种是通过计算方法，一种是通过给定的经纬度文件构建查找表。通过计算的方法经查存在一定的问题，其基本思路是通过计算每个行列号的经纬度坐标文件，然后生成控制点对影像进行校正。具体方法可以自行摸索。这里主要探讨根据官网给出的.raw文件，通过matlab读取，然后生成查找表进行校正。代码如下：

clc,clear;
dim = 10992;
fid=fopen('D:\FY-4\data\FullMask_Grid_1000.raw','rb');
data = fread(fid,dim*dim*2,'double');
b1 = reshape(data,[2,dim,dim]);
b1_lat = reshape(b1(1,:,:),dim,dim);
b1_lon = reshape(b1(2,:,:),dim,dim);
lat = b1_lat';
lon = b1_lon';
fclose(fid);
%%
x = lon(1619:2018,6656:7055);                                              %截取一个子区间
y = lat(1619:2018,6656:7055);
%%
a = [];
b = [];
for i = 1:400                                                              %将图像归一化，主要是为了方便将数值矩阵写成图像
    for j = 1:400                                                          %当然也可以试试mat2gray这个函数
        a(i,j) = (x(i,j) - min(x(:)))/(max(x(:))-min(x(:)));               %由于不清楚上面这个函数用的哪种归一化算法，所以自己写了个函数
        b(i,j) = (y(i,j) - min(y(:)))/(max(y(:))-min(y(:)));
    end
end
imwrite(a,'D:\FY-4\data\xxx_longitude.tif');                               %将查找表文件输出
imwrite(b,'D:\FY-4\data\yyy_latitude.tif');                                %注意输出后的文件在0-255之间，还需将数据还原为原始值，具体方法此处略去

由于全圆盘数据量较大，所以这里只选取了400*400个子区间做实验。得到的X,Y查找表文件如下：

上图分别为经纬度的查找表文件，有了查找表文件后，就可按照如下博客内容，在ENVI里进行基于GLT方法的几何校正工作，博客链接如下：http://blog.sina.com.cn/s/blog_764b1e9d0100qzow.html。

最后附上校正成功后的结果：