Satpy基礎系列教程(1)-FY4A AGRI L1數據處理
- 2020 年 2 月 17 日
- 筆記
以下文章來源於氣象雜貨鋪 ，作者朝曦dawn
Satpy[1]目前支援的衛星數據[2]有50種(MSG, Himawari 8, GOES-R, MODIS, Sentinel- 1/2/3/5, SNPP等)。
本文以最近朝曦dawn[3]添加的風雲4A(FY4A) AGRI L1數據為例。
Notebook[4]已放在GitHub上，供大家學習。
FY4A AGRI L1數據有兩種類別：
1.全圓盤
FY4A-_AGRI–_N_DISK_1047E_L1-_FDI-_MULT_NOM_20190807060000_20190807061459_4000M_V0001.HDF
2.中國區
FY4A-_AGRI–_N_REGC_1047E_L1-_FDI-_MULT_NOM_20190807045334_20190807045750_1000M_V0001.HDF
全圓盤標識：DISK , 中國區標識：REGC.
風雲數據鏈接
1.實時數據[5](30天內)及文檔
2.歷史數據[6] (2018-03-12 — )
3.FY4A官方應用平台[7]
安裝
推薦通過conda安裝，簡單快捷。
$ conda install -c conda-forge satpy
定標
FY4A AGRI L1有以下三種定標方式:
1.原始數字量化值 (所有通道)
2.反射率 (C01 – C06)
3.輻射率和亮溫 (C07 – C14)
讀取數據
只需幾行即可讀取數據：
import os, glob from satpy.scene import Scene # 載入FY4A文件 filenames = glob.glob('/xin/data/FY4A/20190807/FY4A-_AGRI*4000M_V0001.HDF') # 創建scene對象 scn = Scene(filenames, reader='agri_l1') # 查看可用的通道 scn.available_dataset_names()
輸出結果:
['C01', 'C02', 'C03', 'C04', 'C05', 'C06', 'C07', 'C08', 'C09', 'C10', 'C11', 'C12', 'C13', 'C14', 'satellite_azimuth_angle', 'satellite_zenith_angle', 'solar_azimuth_angle', 'solar_glint_angle', 'solar_zenith_angle']
讀取指定通道數據：
# 以紅外通道為例 ir_channel = 'C12' scn.load([ir_channel])
繪圖
全圓盤圖
一行命令即可保存或顯示圖片：
# 保存到文件 # scn.save_dataset(ir_channel, filename='{sensor}_{name}.png') # 在notebook中顯示 scn.show(ir_channel)
全圓盤真彩色圖
首先查看可用的合成選項：
scn.available_composite_names()
輸出結果顯示，支援真彩圖：
['ash', 'dust', 'fog', 'green', 'green_snow', 'ir108_3d', 'ir_cloud_day', 'natural_color', 'natural_color_sun', 'night_background', 'night_background_hires', 'overview', 'overview_sun', 'true_color']
於是我們可以直接載入真彩色數據並繪圖：
# 註：這步需要大記憶體 (取決於cpu核數) # 可查看FAQ關於記憶體的討論： # https://satpy.readthedocs.io/en/latest/faq.html composite = 'true_color' scn.load([composite]) scn.show(composite) # scn.save_dataset(composite, filename='{sensor}_{name}.png')
特定區域圖
以下以颱風利奇馬為例。
首先，我們需要定義地圖投影和區域，然後將數據投影到該區域上。
方式一
用Pyresample[8]來定義區域。
from pyresample import get_area_def area_id = 'lekima' x_size = 549 y_size = 499 area_extent = (-1098006.560556, -967317.140452, 1098006.560556, 1026777.426728) projection = '+proj=laea +lat_0=19.0 +lon_0=128.0 +ellps=WGS84' description = "Typhoon Lekima" proj_id = 'laea_128.0_19.0' areadef = get_area_def(area_id, description, proj_id, projection,x_size, y_size, area_extent)
方式二
用coord2area_def.py[9]程式來直接生成區域定義。
比如用python coord2area_def.py lekima_4km laea 10 28 118 138 4，即可得到之前定義的利奇馬區域：
lekima_4km: description: lekima_4km projection: proj: laea ellps: WGS84 lat_0: 19.0 lon_0: 128.0 shape: height: 499 width: 549 area_extent: lower_left_xy: [-1098006.560556, -967317.140452] upper_right_xy: [1098006.560556, 1026777.426728]
然後將該定義拷貝到$PPP_CONFIG_DIR/areas.yaml中，即可直接調用。
# 如果你已經添加區域到areas.yaml中,可直接調用: os.environ['PPP_CONFIG_DIR'] = '/yin_raid/xin/satpy_config/' lekima_scene = scn.resample('lekima_4km') # 否則需要使用之前定義的區域: # lekima_scene = scn.resample(areadef)
lekima_scene.show(composite) # lekima_scene.save_dataset(composite, filename='{sensor}_{name}_resampled.png')
如果想利用自定義的colormap來生成影像（如下圖），請參閱關於enhancement的notebook（正在忍餓趕稿中）。
References
[1] Satpy: https://satpy.readthedocs.io/en/latest/ [2] 支援的衛星數據: https://satpy.readthedocs.io/en/latest/index.html#reader-table [3] 朝曦dawn: https://dreambooker.site/ [4] Notebook: https://github.com/zxdawn/FY-4/tree/master/satpy/examples [5] 實時數據: https://fy4.nsmc.org.cn/data/en/data/realtime.html [6] 歷史數據: http://satellite.nsmc.org.cn/PortalSite/Data/Satellite.aspx [7] FY4A官方應用平台: http://rsapp.nsmc.org.cn/geofy/ [8] Pyresample: http://pyresample.readthedocs.org/ [9] coord2area_def.py: https://github.com/pytroll/satpy/blob/master/utils/coord2area_def.py
