Satpy[1]目前支持的卫星数据[2]有50种(MSG, Himawari 8, GOES-R, MODIS, Sentinel- 1/2/3/5, SNPP等)。html
本文以最近朝曦dawn[3]添加的风云4A(FY4A) AGRI L1数据为例。python
Notebook[4]已放在GitHub上，供你们学习。
git
FY4A AGRI L1数据有两种类别：github
1.全圆盘微信
FY4A-_AGRI--_N_DISK_1047E_L1-_FDI-_MULT_NOM_20190807060000_20190807061459_4000M_V0001.HDFapp
2.中国区学习
FY4A-_AGRI--_N_REGC_1047E_L1-_FDI-_MULT_NOM_20190807045334_20190807045750_1000M_V0001.HDFflex
全圆盘标识：DISK , 中国区标识：REGC.ui
风云数据连接spa
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, globfrom 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 = 549y_size = 499area_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
本文分享自微信公众号 - 气象杂货铺（meteogs）。
若有侵权，请联系 support@oschina.cn 删除。
本文参与“OSC源创计划”，欢迎正在阅读的你也加入，一块儿分享。