matlab构造圆盘,基于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。

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