前言
最近好多网友有FY-4的几何校正需求,就来更新一下。
数据准备
L1级的FY-4数据是HDF格式的全圆盘数据,数据集里面没有经纬度信息。每个位置对应的经纬度数据需要单独下载,如下图。这里我下载了4 km的经纬度查找表。
文件下载解压后,是一个.raw格式的二进制文件。
下面用python解析这个文件。
import numpy as np
file = r'F:\Py_project\geoFY-4A\FullMask_Grid_4000.raw'
raw_image = np.fromfile(file)
加载进来是一个长度为15103008的一维数组。但是这里面存储的应该是两个波段的数据,分别是经度和纬度。
通过阅读其他博客可以得知它的存储方式是BIP。
BIP的存储方式如下,其中Band1代表纬度数据,Band2代表经度数据,(0,1)代表行列分别为0,1。也就是说这一个长度为15103008的一维数组里面交错存储着每个像元的纬度和经度信息。
所以reshape
成两列转换成如,两列数据,分别对应每个像元的纬度和经度数据。
raw_image = raw_image.reshape(-1, 2)
根据官方文档,可知4 km的影像数组大小为(2748, 2748).
所以可以继续reshape
提取出纬度和经度数组。因为圆盘之外的经纬度数据为999999,为了显示,我转换成了np.nan
lat = raw_image[:, 0].reshape(2748, 2748)
# lat[lat>999] = np.nan
plt.figure()
plt.imshow(lat, cmap=plt.get_cmap('gray'))
plt.colorbar()
plt.title('Latitude')
lon = raw_image[:, 1].reshape(2748, 2748)
# lon[lon>999] = np.nan
plt.figure()
plt.imshow(lon, cmap=plt.get_cmap('gray'))
plt.colorbar()
plt.title('Longitude')
如下图,可以看出提取的经纬度信息是正确的,即经度自西向东增大,纬度自上向下减小。
然后,为了方便存储,输出为HDF文件。
with h5.File(r'F:\Py_project\geoFY-4A\Geo.hdf', 'w') as f:
f['Lat'] = lat.astype('float16')
f['Lon'] = lon.astype('float16')
f.close
在ENVI中打开输出的经纬度HDF文件。
如果想要获取其他分辨率数据的经纬度数据只需要修改(2748, 2748)为对应的行列数。
代码和我提取出来的4KM 经纬度数据可以在这里下载到提取码:bznc
几何校正
至此,我们已经有了FY-4 L1圆盘数据的经纬度信息了,所以就可以在ENVI中通过构建GLT文件来几何校正。具体GLT构建及校正可自行百度,或者看这篇博客。
存在的问题
我输出HDF格式数据的目的起初是为了用gdal.Warp
来校正的,但是一直校正不成功。然后我就使用ENVI来构建GLT,发现会报错。想到可能是我把圆盘外的坐标信息修改成了NaN的原因,然后我就取消修改,输出了默认的值。发现仍然报错,因为默认的值是999999,ENVI里面识别为INF,即无穷大。如果修改为其他的数值,比如300,仍然报错。进行如上修改后用gdal校正也是不可行的。目前看来,想要通过GLT对圆盘数据进行全部校正是不可行了。所以只能先对经纬度、L1数据进行裁剪,针对某一区域(不能裁剪到圆盘外的点)单独构建GLT,然后来校正了。如果考虑到批量的话,可以使用ENVI的二次开发函数,批量使用GLT进行校正,后续有时间会更新。