satpy 处理卫星 FY4A 数据

1. 读取数据并画图

   import os
   import glob
   from datetime import datetime, timedelta
   from satpy.scene import Scene
   from pyresample import get_area_def
   import warnings
   
   warnings.filterwarnings('ignore')
   
   
   def fileNameSplit(fileNs):
       fSplit = []
       for fpn in fileNs:
           fp = os.path.dirname(fpn)
           fn = os.path.basename(fpn).split('_P_')[-1]
           nfpn = os.path.join(fp, fn)
           if not os.path.exists(nfpn):
               os.rename(fpn, nfpn)
           fSplit.append(nfpn)
       return fSplit
   
   
   def time_format(info, format="%Y%m%d%H%M%S", format2="%Y%m%d%H%M00"):
       dt = datetime.strptime(info, format)
       strdt = dt + timedelta(hours=8)
       return datetime.strftime(strdt, format2)
   
   
   def get_fy4a(file_paths, channel_list):
       fileTime = time_format(os.path.basename(file_paths[0]).split('_')[-3])
       # 创建scene对象
       scn = Scene(file_paths, reader='agri_l1')
       # 查看可使用的通道
       if 'C03' in channel_list:
           print(scn.available_dataset_names())
       if 'true_color' in channel_list:
           # 查看可用的合成选项
           print(scn.available_composite_names())
       # 读取指定通道数据
       scn.load(channel_list)
       # 打印加载通道后 xarray.Dataset 的数据集
       # print(scn)
   
       """
           >python coord2area_def.py china merc 5 54 60 139 3
           ### +proj=merc +lat_0=29.5 +lon_0=99.5 +ellps=WGS84
       
           china:
             description: china
             projection:
               proj: merc
               ellps: WGS84
               lat_0: 29.5
               lon_0: 99.5
             shape:
               height: 2194
               width: 2931
             area_extent:
               lower_left_xy: [-4397119.886334, 553583.846816]
               upper_right_xy: [4397119.886334, 7135562.567523]
       """
       # 画自定义区域图片
       area_id = 'china'
       area_name = "myarea"
       proj_id = 'china_99.5_29.5'
       proj4_args = '+proj=merc +lat_0=29.5 +lon_0=99.5 +ellps=WGS84'
       height = 2194
       width = 2931
       area_extent = [-4397119.886334, 553583.846816, 4397119.886334, 7135562.567523]
       # 定义地图投影和区域，然后将数据投影到该区域上
       areadef = get_area_def(area_id, area_name, proj_id, proj4_args, width, height, area_extent)
       china_scene = scn.resample(areadef)
   
       channel_ele = {'C03': 'vis', 'C09': 'vap', 'C12': 'ir', 'true_color': 'thr'}
       for cele in channel_list:
           pic_dir = r'./'
           pic_name = f'fy4a_{channel_ele[cele]}_{fileTime}.png'
           pic_path = os.path.join(pic_dir, pic_name)
           if not os.path.exists(pic_path):
               print(f'生成图片 {pic_name}')
               # 保存图片
               china_scene.save_dataset(cele, filename=pic_path)
   
   
   if __name__ == '__main__':
       t1 = datetime.now()
   
       grayFilePathList = ['./FY4A-_AGRI--_N_REGC_1047E_L1-_FDI-_MULT_NOM_20211020022336_20211020022753_4000M_V0001.HDF',
                           './FY4A-_AGRI--_N_REGC_1047E_L1-_GEO-_MULT_NOM_20211020022336_20211020022753_4000M_V0001.HDF']
       print(grayFilePathList)
       channelList = ['C03', 'C09', 'C12']
       get_fy4a(grayFilePathList[:1], channelList)
       t2 = datetime.now()
       print(t2 - t1)
       print('------')
       thrFilePathList = ['./FY4A-_AGRI--_N_REGC_1047E_L1-_FDI-_MULT_NOM_20211020022336_20211020022753_4000M_V0001.HDF',
                           './FY4A-_AGRI--_N_REGC_1047E_L1-_GEO-_MULT_NOM_20211020022336_20211020022753_4000M_V0001.HDF']
   
       print(thrFilePathList)
       channelList1 = ['true_color']
       get_fy4a(thrFilePathList, channelList1)
       t3 = datetime.now()
       print(t3 - t2)

2. reader='agri_l1'   

   版本<=0.36.0  'agri_l1'     >=0.37.0   'agri_fy4a_l1'   其他请看网址

   https://satpy.readthedocs.io/en/stable/index.html#reader-table

3. HTTPSConnectionPool(host='zenodo.org', port=443): Max retries exceeded with url: /record/1288441/files/pyspectral_atm_correction_luts_no_aerosol.tgz

   # 原因下载依赖文件失败
   # 手动下载地址

   https://zenodo.org/record/3824535/files/pyspectral_rsr_data.tgz
   
   https://zenodo.org/record/1288441/files/pyspectral_atm_correction_luts_no_aerosol.tgz

   
   win
   pyspectral_rsr_data.tgz  解压放到下面地址
   C:\Users\admin\AppData\Local\pytroll\pyspectral
   
   pyspectral_atm_correction_luts_no_aerosol.tgz  解压放到下面地址
   C:\Users\admin\AppData\Local\pytroll\pyspectral\rayleigh_only
   
   linux
   pyspectral_rsr_data.tgz  解压放到下面地址
   /root/.local/share/pyspectral
   
   pyspectral_atm_correction_luts_no_aerosol.tgz  解压放到下面地址
   /root/.local/share/pyspectral/rayleigh_only

    

4. Don't know how to open the following files: {'./FY4A.......HDF'}

   # 按我猜测，应该是satpy 和其依赖包版本的问题
   h5py==3.7.0
   numpy==1.21.6
   satpy==0.36.0
   pyresample==1.22.0
   scikit-image==0.19.3
   pyspectral==0.12.0

    

5. 参考地址

   官网：https://satpy.readthedocs.io/en/stable/readers.html#available-readers
   博客：https://cloud.tencent.com/developer/article/1584969
   牛人：https://github.com/aiwei169/FY4A/blob/main/codes/fy4a.py