3月15日近10年来最强的沙尘暴袭击了北京。关于此次沙尘暴的天气分析已经非常多了，本文不会分析相关的天气背景，主要从技术方面讲一下如何利用satpy处理卫星数据，从卫星视角看一下此次沙尘过程的演变。
本文的数据为 Himawari-8 静止卫星L1b产品：
代码语言：javascript复制from glob import glob
from datetime import datetime
import dask
import satpy
from satpy import Scene
from satpy.resample import get_area_def
from pyresample import create_area_def
from pyresample.geometry import AreaDefinition
以下是日期解析和数据处理函数：
代码语言：javascript复制def process(f):
scn = Scene(f, reader='ahi_l1b_nc')
scn.load([composite_name])
area_def = create_area_def('china',
' proj=eqc pm=180',
area_extent=area_extent,
units='degrees',
resolution=res
)
lcn = scn.resample(area_def)
#lcn.show(composite_name)
lcn.save_dataset(composite_name, f'dust/{composite_name}_{parse_date(f[0]):%Y%m%d%H%M}.tiff')
parse_date = lambda f: datetime.strptime(''.join(f.split('_')[2:4]), '%Y%m%d%H%M')
因为所处理的数据比较多，上述 process 函数直接将处理后的结果保存为 .tiff 格式文件。如果你想单独看某一个时刻的卫星图像，可以在 save_dataset 之前使用 show 显示图片。
注意：satpy 官方暂时不支持 Himawari-8 卫星的 L1b 数据接口，需要单独添加接口。
真彩色图
利用satpy绘制真彩色图非常方便，给定 composite 参数即可，同时给定经纬度范围限制图片显示范围。
fp = '/data/phd/satellite/himawari8/work/202103/15/'
fn = 'NC_H08_20210315_0*_R21_FLDK.06001_06001.nc'
files = glob(fp fn)
res = 0.01
area_extent = (80, 20, 140, 70)
composite_name = 'true_color'
由于需要绘制的图形比较多，为了加快绘图速度，使用 dask 并行绘图：
%%time
tasks = dask.delayed(process([f]) for f in files)
tasks.compute()
3月15日0500UTC Himawari-8真彩色图
沙尘
其实从真彩色图上已经能够看出沙尘的发展了。此外，satpy 还提供了 dust 合成产品，主要是通过不同的光谱通道对粒子反射特性进行计算。
composite_name = 'dust'
tasks = dask.delayed(process([f]) for f in files)
tasks.compute()
3月15日0500UTC 沙尘合成产品
多说几句，除了 真彩色图 和 dust 的合成产品之外，satpy 还支持很多合成产品，比如 fog 和 convection 等，处理方法是类似，只需要更改 composite 参数即可。
下面是此次沙尘暴发展视频 http://mpvideo.qpic.cn/0bf2keabgaaakiaajuksyrqfauodcniqaeya.f10002.mp4?
—END—
