Featured image of post 利用 Python 处理向日葵卫星 netCDF 数据生成真彩色云图
Python

利用 Python 处理向日葵卫星 netCDF 数据生成真彩色云图

初

我早在去年 8 月份便在 JAXA 申请过向日葵卫星 ftp 账户，但一直到最近才看到 JAXA 给我的回信(居然在垃圾邮箱里)。
因此，我就开始着手对 向日葵卫星 L1 数据进行处理。在此期间，学习到了很多东西，因此以此博客记录处理过程。
源代码放置于 github 上，链接:
https://github.com/storm-1614/Himawari/

申请过程

向日葵卫星的申请网站是:https://www.eorc.jaxa.jp/ptree/registration_top.html用邮箱提交然后填写点个人资料很快就回复了。记得检查垃圾邮箱，不要像我一样傻傻的等了半年都忘掉了

FTP 的使用

登录

JAXA 在回信中会给你提供 ftp 账户。而后，我们便要用这个 ftp 账户进行数据的下载。 ftp 软件我使用的是 FileZilla，是全平台的软件，在 Win、Linux 皆有。在 ArchLinux 中，下载方法是:

1

sudo pacman -S filezilla

在进入软件输入网址、用户名和密码后就可以登陆了。
以及你可以在 File>Site Manager 中对登录信息保存，方便下次登陆。

目录

我使用的是 netCDF 文件，它存储在 /jma/netcdf，该目录下有从 2015 年 7 月份到现在为止的数据，一个文件例子是NC_H09_20240219_0000_R21_FLDK.06001_06001.nc
这里的 H09是卫星编号，即 Himawari 9
20230219_0000 是时间
R21_FLDK是区域，该网站提供两个区域。R21是全圆盘、r14是日本区域。
而06001我不知道代表什么，但是 06001 是最清晰文件体积最大的。

向日葵9号信息

向日葵9号是向日葵8号的备份星，上面的传感器是一样的，在维基百科中文上对向日葵8号卫星的波段介绍较为详细因此在这提供: 向日葵9号波段信息 从上图表格能够看出，1～3频段分别是蓝绿红的可见光波段。刚好与 RGB 相反，这是要注意的。而后4～16是波长逐渐增长的红外波段。我关注第13波段，之后要用该波段补上夜晚可见光的缺漏。

使用的Python库

这里使用到的库是 numpy、netCDF4、PIL 以及用来处理参数的sys。

numpy和netCDF4需要自己下载，这里同样提供的是 ArchLinux 下载方法：

1

sudo pacman -S python-numpy python-netcdf4 python-pillow

下载好后，就可以运行我提供在 github 上的程序了，输出的图片如下:
因为图片过大（39MB)分辨率6000 * 6000。所以就放个截图。

编写过程

处理 nc 源文件

这里也就是开端，但是netCDF这种应用在气象等科学上的文件我没有接触过，因而极为艰难，这里就说明下 netCDF 处理的方法。 需要说明的是，导入模块的时候我用了as进行简化，这里放上代码:

1
2

import netCDF4 as nc
import numpy as np

PIL 导入了三个模块:

1

from PIL import Image, ImageEnhance, ImageOps

1
2
3

data = "./NC_H09_20240131_2240_R21_FLDK.06001_06001.nc"

nc_data = nc.Dataset(data)

这里的 data 变量用来存储文件目录，而后用 nc.Dataset把处理的数据存储到 nc_data 内。 之后这个 data 可以用 sys.argv 套下方便以后的处理。

这个时候的 nc_data就存储了处理过的 .nc 数据，通过输出 nc_data.variables.keys()可以看到如下：

1
2
3
4

dict_keys(['latitude', 'longitude', 'band_id', 'start_time', 'end_time', 'geometry_parameters', 'alb
edo_01', 'albedo_02', 'albedo_03', 'sd_albedo_03', 'albedo_04', 'albedo_05', 'albedo_06', 'tbb_07', 
'tbb_08', 'tbb_09', 'tbb_10', 'tbb_11', 'tbb_12', 'tbb_13', 'tbb_14', 'tbb_15', 'tbb_16', 'SAZ', 'SA
A', 'SOZ', 'SOA', 'Hour'])

这里主要关注 albedo_01~albedo_06和 tbb_07~tbb_16这16个键存储的就是 16 波段的数据。对应波段信息在这： 向日葵9号信息

可见光部分

刚开始我所想的是处理个可见光的云图，网上的文章也不少。但这其中因为数学和编程经验的不足，花费了太多时间在学习上，最终我实现的可见光处理的函数是:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21

def vis(nc):
	# 获取可见光波段数据
    blue = nc.variables['albedo_01'][:]
    green = nc.variables['albedo_02'][:]
    red = nc.variables['albedo_03'][:]
    
	# 分别归一化
    for i in [blue, green, red]:
        i = i / np.max(i)
        
	# 处理成一个三位数组
    rgb = np.dstack((red, green, blue))
    
	# 转成 UTF-8 编码
    rgb = (rgb * 255).astype(np.uint8)
    
	# 处理成图像
    imgVis = Image.fromarray(rgb)
    print("可见光处理完成")
    
    return imgVis

输入参数是 nc 文件，而返回的是imgVis存储了图像数据
从之前便知道nc.variables中存储了 16 个波段的数据，然后我用 [albedo_01][:]分别提取 1~3 波段对应的数据。用颜色存储。变量类型是 <class 'numpy.ma.core.MaskedArray'>，粗略理解是个数组，至于 MaskedArray这个中文翻译是掩码数组，不太理解。
之后很大程度依赖于网上的博客和 ChatGPT 对此理解。
这里简单讲解下:

1
2
3

# 分别归一化
for i in [blue, green, red]:
	i = i / np.max(i)

这里是把一个数组除以这个数组里面的最大值，以达到归一化，这是一个算法，现在理解程度有限。

1
2

# 处理成一个三位数组
rgb = np.dstack((red, green, blue))

这里是把三个波段数据组合成一个三位数组。
要注意的是 这里应该按 红绿蓝排列，向日葵卫星的颜色排列方式刚好相反。

1
2

	# 转成 UTF-8 编码
    rgb = (rgb * 255).astype(np.uint8)

这是因为我发现直接原数据给 PIL 它会报错，网上查询在知道是编码问题，这是 ChatGPT 给出的解决方案。

1

imgVis = Image.fromarray(rgb)

这是用来把上面存储成三维数组的 rgb 变成图片，用的是 PIL 的Image.fromarray 函数。

之后便是返回这个变量。

红外部分

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

def ir(nc):
    wave = nc.variables['tbb_13'][:]
    wave = ((wave - np.min(wave)) * (255 - 0)) / (np.max(wave) - np.min(wave)) + 255

    wave = wave.astype(np.uint8)

    imgIR = Image.fromarray(wave)
    imgIR = ImageOps.invert(imgIR)

    print("红外处理完成")
    return imgIR

这是红外的处理函数，内容差不多，不过多介绍。
用的是 13 波段，也就是向日葵卫星用来做 IR1 的波段。
这里的归一化是用 ChatGPT 提供的方法归一。
要注意的是，这个出来的图片云的黑色的，用反色出来才是正常的，用的是 ImageOps.invert函数 返回一致，是 处理的 IR 图片信息

混合图片

混合图片是为了达到夜间也能看到云图，这里提供函数代码:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

def blend(vis, ir):
    if ir.size != vis.size:
        ir = ir.resize(vis.size)

    if ir.mode != vis.mode:
        ir = ir.convert(vis.mode)

    img = Image.blend(ir, vis, alpha=0.8)
    
    print("合并完成")
    return img

前两个 if 判断是要保证 vis 和 ir 是一致的尺寸和模式，不然混合出来是奇奇怪怪的。

后面就是用 Image.blend 函数进行混合。出来的样子已经很不错了。

图片色彩

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15


def color(img):
    #亮度
    img_1 = ImageEnhance.Brightness(img)
    img_1 = img_1.enhance(1.3)

    #对比度
    img_2 = ImageEnhance.Contrast(img_1)
    img_2 = img_2.enhance(1.1)
    
    #饱和度
    img_3 = ImageEnhance.Color(img_2)
    img = img_3.enhance(2)
    print("色彩处理完成")
    return img

这里的代码写得很烂，随便看看吧（真的像答辩一样烂）
ImageEnhance模块是现找现学的，很多思路没理明白，以后可以写篇博客介绍介绍，然后顺便改改这垃圾代码。

保存图片

1

color(blend(vis(nc), ir(nc))).save('./1.png')

是的，就是这么干脆利落。一行调用函数直接完事了！ 理解起来其实挺好理解的，只需要知道 PIL的save 函数以及我上面写的各个函数就能知道。

参考资料

python处理葵花8 netCDF4（nc格式）数据
向日葵8号维基百科
Pillow Docs
JAXA Himawari Monitor
Himawari-8/9 Operational Information

附记

北京时间2016年9月14日7时 向日葵8号西北太平洋云图 storm1614制作 北京时间2023年7月27日14时20分 向日葵9号西北太平洋卫星云图 storm1614 制作

上面两张图片，饱含了太多意义了……

当向日葵 8 号照下第一张照片的时候，他四年级。天空的跑马云、学校的停课通知与即将到来的中秋假期一同构成了那个带着些许不安的早晨。

那时的他还在用手笔着床边的中国地图，通过昨晚从父亲手机上得知的莫兰蒂卫星云图笔画着它在哪呢？狭窄的海峡，模糊的童年。一切还是懵懂，谁知道呢？

第二张，是去年杜苏芮登陆泉州前一天的图片。彼时的他将升高二，废寝忘食的一周……

愿未来如浩瀚大气长河奔涌不息，愿少年终不负少年梦。

萌ICP备20241614号