FY4A数据读取与插值

先说结果：

已实现FY4A总共57种产品的数据读取，并且进行了：

经纬度矩阵的计算、中国区域全圆盘的转换、有效区域与无效区域的标定，根据指定经纬网格进行插值等功能，后续将会更新定标、定位、几何校正、数据反演等功能

代码在github: https://github.com/aiwei169/FY4Ahttps://github.com/aiwei169/FY4A

下载地址：

https://download.csdn.net/download/aiwei169/25622062https://download.csdn.net/download/aiwei169/25622062

使用方法如下，现任选一个FY4A文件（用户使用时，只需要输入路径即可，程序会自动识别是哪一个FY4A产品），比如输入CLT产品，并插值到中国区域

import utils
import components
import logging
from fy4a import FY4A


arg = components.Arg()
arg_parsed = arg.arg_parse('-d'.split())  # 代码输入参数
# arg_parsed = arg.arg_parse()  # 命令行输入参数

fy4a = FY4A(arg_parsed)


# FY4A文件路径(任写一种产品即可)
file_path = r'/home/developer_13/FY4A/codes/test_FY4A_data/' \
            r'FY4A-_AGRI--_N_DISK_1047E_L2-_CLT-_MULT_NOM_20210923010000_20210923011459_4000M_V0001.NC'
# 存储nc的路径
sv_path = r'/home/developer_13/FY4A/codes/CLT.nc'
# 需要插值的经纬度
new_lon, new_lat = utils.gen_lat_lon(70, 140, 60, 0, 0.04)
# 得到数据块，并存储为nc
fy4a.get_data_from_data_name(file_path, 'CLT', new_lon, new_lat, sv_path=sv_path)
# 数据对象
logging.debug(fy4a.read_FY4A(file_path))

输出一个dataset结果：

输出文件：

 METEINFO画图查看：

再试试HDF文件，比如我们随机选择一个2000m的L1的产品：FY4A-_AGRI--_N_DISK_1047E_L1-_FDI-_MULT_NOM_20210922084500_20210922085959_2000M_V0001.HDF，也是插值到中国区（从全圆盘插值到中国区），选择波段NOMChannel01（如果你想多个波段都做，写个循环就可以了）

代码还是一样的，只需要修改一下文件名，别的不变，

输出结果：

输出nc文件，在meteinfo展示：

 所有文件读取之后，都是生成一个xarray下的dataset对象，方便用户进行其他的操作

经纬度转换的方法，代码实现了从经纬度转行列、从行列转经纬度、读取raw文件等

需要注意几个问题：

1.官方的解释raw文件的说明有问题，如：

这里应该是前8字节是纬度值，后8字节是经度值，我检查了所有的raw文件，都是这样的

 REGC官方没有给出读取的办法，但是都在文件中给了一个属性：

# L2 NC文件
ds = xr.open_dataset(file_path)
l_bias = ds['geospatial_lat_lon_extent'].attrs['begin_line_number']
c_bias = ds['geospatial_lat_lon_extent'].attrs['begin_pixel_number']

# L1 HDF文件（GEO除外）
ds = xr.open_dataset(file_path, engine='netcdf4')
l_bias = ds.attrs['Begin Line Number']
c_bias = ds.attrs['Begin Pixel Number']

l_bias和c_bias分别表示中国区数据块在标称投影中的行列起始位置，

参考代码里的定义：

    def l_c_to_Lat_lon(self, data: ndarray, l_bias: int, c_bias: int, **kwargs) -> Tuple[ndarray, ndarray]:
        """
        根据data的行列，计算得到data每个数据点的经纬度，并插值到指定经纬度网格，存储为nc数据

        Args:
            data:       读取到的行列数据
            l_bias:     行偏移，即：(data的第l行，相当于标称投影中的 (l + l_bias) 行)
            c_bias:     列偏移，即：(data的第c列，相当于标称投影中的 (c + c_bias) 列)
        Returns:
            原始数据块对应的经纬度矩阵
        """

或者：

    def lat_lon_to_l_c(self, new_lat: ndarray, new_lon: ndarray, l_bias: int, c_bias: int, res: str) -> Tuple[ndarray, ndarray]:
        """
        根据指定的经纬度网格，获取数据

        Args:
            data:       原始数据块
            new_lat:    要插值到的纬度数列，m维向量
            new_lon:    要插值到的经度数列，n维向量
            l_bias:     行偏移，即：(标称投影中的第l行，相当于data中的 (l - l_bias) 行)
            c_bias:     列偏移，即：(标称投影中的第c列，相当于data中的 (c - c_bias) 列)
            res:        原始数据的分辨率
        Returns:
            指定经纬度网格下的行列矩阵
        """

l_bias:行偏移，假设读取了数据块data
标称投影中的第m行，相当于data中的 (m - l_bias) 行
data的第n行，相当于标称投影中的 (n + l_bias) 行

c_bias:列偏移，假设读取了数据块data
标称投影中的第p列，相当于data中的 (p - c_bias) 列)
data的第q列，相当于标称投影中的 (q + c_bias) 列)

如果读取到了DISK数据，那么这里的l_bias和c_bias，读取到的都会是0

至此，我们可以得到数据块中每一个点对应的经纬度，通过插值就可以得到研究区的数据了

向dataset对象中添加了经度矩阵、纬度矩阵和flag标定，举例说明如下：

取产品：FY4A-_AGRI--_N_REGC_1047E_L1-_FDI-_MULT_NOM_20210922083836_20210922084253_2000M_V0001.HDF

是一款中国区2KM的产品

读取之后输出如下：

注意这里：

代码计算并添加了这三个数据块，他们的大小与主产品一致（NOMChannel01、NOMChannel02等）

比如说现在想知道这个数据的纬度矩阵，那么：

ds = fy4a.read_FY4A(file_path)
lat = ds['lat'].values
logging.debug(lat)

即可得到lat矩阵

原始HDF里是没有lat lon信息的，这里计算并添加了进去