云图渲染基础教程——以Python实现MODIS红外与伪VIS云图绘制为例

我不是Carl2

FC-Leo 发表于 2024-4-14 10:14
话说vis and RGB(Ture Colour)也能用以上这几个库完成吗？

vis是可以的
True Color有大气校正，应该需要另外的库

FC-Leo

我不是Carl2 发表于 2024-4-14 12:03
vis是可以的
True Color有大气校正，应该需要另外的库

还有那注册咋办，试了N遍都不行

FC-Leo

还有就是想问一下Carl 用的是3.1几版本？我用3.12貌似很多都适应不了

esseco

FC-Leo 发表于 2024-4-14 12:42
还有就是想问一下Carl 用的是3.1几版本？我用3.12貌似很多都适应不了

3.12很多库都还没适配，建议用3.11或3.10

我不是Carl2

FC-Leo 发表于 2024-4-14 12:42
还有就是想问一下Carl 用的是3.1几版本？我用3.12貌似很多都适应不了

我是3.11

我不是Carl2

c.        HDF文件读取，MODIS波段简介

卫星云图数据存储的格式相当多样，其中不少格式需要较复杂的脚本才能读取（比如Himawari-8/9最广泛使用的DAT格式）。相对容易读取的数据存储方式包括HDF4，HDF5，NetCDF等。MODIS的L1B数据正是用HDF4格式存储的。HDF4有多种读取方式，这里我们选择使用第三方库pyhdf。HDF5和NetCDF我使用的分别是h5py和netCDF4，在此不作详细介绍。

pyhdf的介绍：https://fhs.github.io/pyhdf/modules/SD.html，下图列举了部分会用到的功能

1. from pyhdf.SD import SD, SDC
   
2. import numpy as np

复制代码

首先导入pyhdf与numpy，此处pyhdf只导入了其中的一部分模块，numpy则是导入时命名为了np，后续反复调用的时候少打几个字母。

1. file = SD(r'E:\Data\MOD021KM.A2015296.0515.061.2017322234443.hdf', SDC.READ)

复制代码

在只读模式（SDC.READ）下打开文件，文件路径请自行替换，这里是用了Patricia10.23 0515z的Terra，此处路径前的r是为了直接拷贝过来的绝对路径能够正常读取。

我们先使用datasets()看看文件中放了哪些数据，这一步相当于列出数据集的“目录”，还没有真正将数据取出来：

1. data = file.datasets()
   
2. print(data)
   
3. print(data.keys())

复制代码

直接print(data)会返回每个数据集的名称与一些附带信息，加上keys()能仅返回数据集名称，以下是print(data.keys())的结果：

1. dict_keys(['Latitude', 'Longitude', 'EV_1KM_RefSB', 'EV_1KM_RefSB_Uncert_Indexes', 'EV_1KM_Emissive', 'EV_1KM_Emissive_Uncert_Indexes', 'EV_250_Aggr1km_RefSB', 'EV_250_Aggr1km_RefSB_Uncert_Indexes', 'EV_250_Aggr1km_RefSB_Samples_Used', 'EV_500_Aggr1km_RefSB', 'EV_500_Aggr1km_RefSB_Uncert_Indexes', 'EV_500_Aggr1km_RefSB_Samples_Used', 'Height', 'SensorZenith', 'SensorAzimuth', 'Range', 'SolarZenith', 'SolarAzimuth', 'gflags', 'EV_Band26', 'EV_Band26_Uncert_Indexes', 'Band_250M', 'Band_500M', 'Band_1KM_RefSB', 'Band_1KM_Emissive', 'Noise in Thermal Detectors', 'Change in relative responses of thermal detectors', 'DC Restore Change for Thermal Bands', 'DC Restore Change for Reflective 250m Bands', 'DC Restore Change for Reflective 500m Bands', 'DC Restore Change for Reflective 1km Bands'])

复制代码

其中'Latitude'，'Longitude'是经纬度数据，'EV_1KM_RefSB'是反射波段数据，这里暂时用不到，'EV_1KM_Emissive'是热辐射波段数据，包含了我们所需要的各个红外波段。

使用select将'Latitude'，'Longitude'，'EV_1KM_Emissive'读取出来，并用get()读取数组，用np.shape()获取每个数组的形状（不过其实print(data)也会将数组形状作附带信息返回）：

1. lon = file.select('Longitude').get()
   
2. lat = file.select('Latitude').get()
   
3. Emissive = file.select('EV_1KM_Emissive').get()
   
4. print(lon, lat, Emissive)
   
5. print(np.shape(lon), np.shape(lat), np.shape(Emissive))

复制代码

可以看到热辐射波段数组的形状是(16, 2030, 1354)，代表有16个不同波段堆叠在一起，每个波段的图像有2030×1354个像素，而数组的值都是0~32767或65535的整数，显然不是我们最终需要的亮温，需经过辐射定标与亮温反演（见2.d）；经纬度数组的值是正确的，但数组形状是(406, 271)，像素比波段数据少，这是因为文件中存储的经纬度降过分辨率，需插值后才能使用（见2.e）。

'EV_1KM_Emissive'这个数据集并非只包含上面读取的数组，HDF文件可以为数据集附上attributes，其中往往包含辐射定标的关键信息，查看方式如下：

1. print(file.select('EV_1KM_Emissive').attributes())
   
2. print(file.select('EV_1KM_Emissive').attributes().keys())

复制代码

以下是加上keys()，只获取attributes名称的结果：

1. dict_keys(['long_name', 'units', 'valid_range', '_FillValue', 'band_names', 'radiance_scales', 'radiance_offsets', 'radiance_units'])

复制代码

其中'long_name'是数据集全称， 'band_names'是波段号，这两个attribute可以在看官方文献看晕掉的时候让你知道读的是什么数据；'radiance_scales'，'radiance_offsets'是辐射定标的两个系数（见2.d），'radiance_units'是辐射定标后得到的辐亮度的单位，有时候不同卫星数据会采用不同的辐亮度单位，反演的亮温如果偏差很大则需要检查该单位并调整黑体辐射公式的形式。

获取attributes中的数组同样是使用get()即可：

1. radiance_scales = file.select('EV_1KM_Emissive').attributes().get('radiance_scales')
   
2. radiance_offsets = file.select('EV_1KM_Emissive').attributes().get('radiance_offsets')

复制代码

MODIS波段介绍：https://modis.gsfc.nasa.gov/about/specifications.php
先前查看热辐射波段attributes的时候16个波段号分别为20,21,22,23,24,25,27,28,29,30,31,32,33,34,35,36。细心的读者可能已经发现了，Band26不在其中。（不细心的作者因为这个在伪VIS教程篇里面把中波红外要用的Band20错写成Band21了，前两天刚勘误）Band26是短波红外反射波段，在2.d中读取指定波段数据的时候需要注意到这一点，不然可能会读错波段。

16个热辐射波段中，Band31用于渲染我们最熟知的IR云图，Band27则是标准水汽波段，Band28是响应层面偏低一些的水汽波段，反映中低层大气的暖心和水汽吸收情况。伪VIS所用到的波段则是Band31，Band32与Band20，波段选择的具体原因可以参考《伪VIS算法原理与教程 - 理论篇》。值得一提的是，Band20与Band21,22的波长非常相近，选择Band20的主要原因是Band20接收的波长宽一些，更容易接收到较冷云顶上微弱的辐射，云顶噪点相对而言少一些。

我不是Carl2

d.        辐射定标与亮温反演

辐射定标将卫星仪器测量的原始数据转化为辐亮度这个物理量，而亮温反演通过黑体辐射公式将辐亮度转化为亮温。
从卫星原始数据计算出辐亮度是个很复杂的过程，需要考虑卫星本身的电子元件问题等等，好在L1B数据已经将这些效应处理掉了，只需要我们做辐射定标的最后一步，把波段数据里存储的整数转化为辐亮度。MODIS User guide（https://mcst.gsfc.nasa.gov/sites ... rra_V6.2.3_Aqua.pdf）中的equation 5.8给出了辐射定标的公式：

其中SI（Scaled Integer）就是波段数据存储的整数，L则是辐亮度。

具体到python代码上：

1. radiance = (Emissive[10] - radiance_offsets[10]) * radiance_scales[10]

复制代码

这里计算了Band31的辐亮度。波段数据和辐射定标系数均按照波段序号排列，跳过Band26这个反射波段，Band31是热辐射波段中第11个波段，而python索引从0开始算，因此索引[10]正对应Band31。

黑体辐射的概念在《伪VIS算法原理与教程》理论篇2.1章节有较为详细的介绍，这里只放由此而来的亮温计算公式：

其中T是开氏度单位的亮温，转换到摄氏度需要-273.15，ℎ是普朗克常数，𝑐是光速，𝑘是玻尔兹曼常数，𝜆是波长，ln()是取自然对数，B𝜆是辐亮度。对常数做一下简化：

带入计算可得在SI制单位下C1=0.0143877，C2=1.191043*10^-16，但红外波段的波长一般用微米表示，且查看attributes（见2.c）可知辐亮度也用到了微米作单位，因此调整单位后C1=1.43877*10^4，C2=1.191043*10^8。
带入Band31中心波长（11.030微米）即可将亮温计算出来，其中**代表指数运算，np.log取的就是自然对数：

1. C1 = 1.4387685*10**4
   
2. C2 = 1.19104356*10**8
   
3. bt31 = C1/11.030 / (np.log(1 + C2/11.030**5/radiance)) - 273.15

复制代码

至此我们已经获得了MODIS Band31的准确摄氏度亮温。

我不是Carl2

e.        经纬度插值

先前提到过MODIS L1B文件中的经纬度数值降过分辨率，云图像素数为2030×1354，而经纬度像素数为406×271。为了正常投影云图，我们需要将经纬度数组插值到与云图像素一一对应。我们可以用 scipy.interpolate.RectBivariateSpline类 实现插值。

RectBivariateSpline类，提供样条插值功能：


导入scipy.interpolate

1. import scipy.interpolate

复制代码

首先提取出经度的形状406×271

1. ny, nx = lon.shape

复制代码

用np.arange创建两个数值为0~405，0~270，间隔为1的一维数组，这代表经纬度原始的矩形网格

1. x = np.arange(nx)
   
2. y = np.arange(ny)

复制代码

再创建数值为0~405，0~270，间隔为271/1354，1/5的一维数组，数组长度为2030，1354，这代表插值后的矩形网格

1. newx = np.arange(0, nx, 271/1354)
   
2. newy = np.arange(0, ny, 1/5)

复制代码

应用RectBivariateSpline类，首先明确原网格上每个点的经度数值，随后用刚生成的数组将其插值到2030×1354的网格上

1. splinelon = scipy.interpolate.RectBivariateSpline(y, x, lon)
   
2. newlon = splinelon(newy, newx)

复制代码

此处newlon就是插值后的经度，要获取插值后的纬度重复最后一个步骤并用lat代替即可。有了一一对应的亮温、经度、纬度三个数组，我们就可以开始绘制投影后的红外云图了。

我不是Carl2

红外云图渲染

a.        色阶

相信各位风迷对色阶并不陌生，色阶将红外云图的每个像素按照其亮温上色。色阶颜色随亮温的变化可以分为渐变（如WV色阶0度以下的部分）和跳变（如IR-BD OW以下的部分），当然许多色阶两者兼有（如IR-RAMMB），黑白的IR和VIS云图本质上也可以看做使用了一类从黑到白的渐变色阶。

左侧是Patricia10.23 0515z Terra Band31未经投影的IR-Rammb图像，右侧是IR-Rammb，IR-BD，WV，IR-BW的色阶条。

色阶编写需要用到matplotlib.colors.LinearSegmentedColormap类（https://matplotlib.org/stable/ap ... mentedColormap.html），最简单的IR-BW色阶长这样：

1. bwdata={'green': [(0, 1, 1),
   
2.                   (1, 0, 0)],
   
3.         'red': [(0, 1, 1),
   
4.                   (1, 0, 0)],
   
5.         'blue': [(0, 1, 1),
   
6.                   (1, 0, 0)]}

复制代码

卫星云图每个像素的颜色可以用RGB值表示，如果想要自己编写色阶，可以拿读取像素RGB值的软件找找感觉。LinearSegmentedColormap的色阶格式也有'green' 'red' 'blue'独立的三组，对应G，R，B的亮度值。IR-BW是黑白色阶，因此三组都是一样的，彩色色阶则会有所不同。

每一行色阶由三个数字组成（称为x, y0, y1），并指定某一亮温节点对应的该颜色（RGB其中之一）亮度。x代表该亮温节点的位置，当色阶渲染范围是-100到50度时，x为0代表-100度，为0.5代表-25度，以此类推。y0, y1代表该颜色亮度，其0~1的范围对应0~255的RGB值。为什么需要用两个数字代表同一个亮温节点上的亮度呢？这是因为每个亮温节点两侧的亮度可以不一样，出现颜色的跳变。y0代表更低亮温那一侧的亮度，y1代表更高亮温那一侧的亮度。在两个亮温节点的间隙中，LinearSegmentedColormap会根据上一行的y1与下一行的y0自动进行渐变，若这两个值相同则无渐变。对于上面的IR-BW色阶来说，就是在-100度 亮度为1 和50度 亮度为0的区间内进行渐变。


下面展示几个更复杂的色阶实例：

1. wvdata={'green': [(0, 0/255, 0/255),
   
2.                   (30/150, 0/255, 0/255),
   
3.                   (45.5/150, 128/255, 128/255),
   
4.                   (58.5/150, 255/255, 255/255),
   
5.                   (65/150, 255/255, 255/255),
   
6.                   (72.5/150, 255/255, 255/255),
   
7.                   (86/150, 128/255, 128/255),
   
8.                   (100/150, 20/255, 255/255),
   
9.                   (1, 1, 1)],
   
10.         'red':   [(0, 128/255, 128/255),
    
11.                   (30/150, 128/255, 128/255),
    
12.                   (45.5/150, 255/255, 255/255),
    
13.                   (58.5/150, 255/255, 255/255),
    
14.                   (65/150, 128/255, 128/255),
    
15.                   (72.5/150, 128/255, 128/255),
    
16.                   (86/150, 0/255, 0/255),
    
17.                   (100/150, 100/255, 255/255),
    
18.                   (1, 1, 1)],
    
19.         'blue':  [(0, 0/255, 0/255),
    
20.                   (30/150, 0/255, 0/255),
    
21.                   (45.5/150, 0/255, 0/255),
    
22.                   (58.5/150, 128/255, 128/255),
    
23.                   (65/150, 128/255, 128/255),
    
24.                   (72.5/150, 255/255, 255/255),
    
25.                   (86/150, 255/255, 255/255),
    
26.                   (100/150, 100/255, 255/255),
    
27.                   (1, 1, 1)]}

复制代码

WV色阶应用于-100~50度的渲染范围，因此x0位置上的N/150就对应-100+N度。WV色阶以渐变为主，因此前几行色阶y0y1是相同的，代表没有颜色的跳变，上一行y1与下一行y0不同，代表在亮温节点之间进行渐变。注意在x0=100/150也就是0度的时候y0y1不同，代表此处出现了从紫色到白色的跳变，这一点可以通过本节第一张图确认。

1. bddata =  {'blue': [(0.0, 0.0, 0.3333333333333333),
   
2.                     (0.12666666666666668, 0.3333333333333333, 0.5294117647058824),
   
3.                     (0.16, 0.5294117647058824, 1.0),
   
4.                     (0.2, 1.0, 0.0),
   
5.                     (0.24, 0.0, 0.6274509803921569),
   
6.                     (0.30666666666666664, 0.6274509803921569, 0.43137254901960786),
   
7.                     (0.38666666666666666, 0.43137254901960786, 0.23529411764705882),
   
8.                     (0.46, 0.23529411764705882, 0.792156862745098),
   
9.                     (0.7266666666666667, 0.42745098039215684, 1.0),
   
10.                     (0.8533333333333334, 0.0, 0.0), (1.0, 0.0, 0.0)],
    
11.           'green': [(0.0, 0.0, 0.3333333333333333),
    
12.                     (0.12666666666666668, 0.3333333333333333, 0.5294117647058824),
    
13.                     (0.16, 0.5294117647058824, 1.0),
    
14.                     (0.2, 1.0, 0.0),
    
15.                     (0.24,0.0, 0.6274509803921569),
    
16.                     (0.30666666666666664, 0.6274509803921569, 0.43137254901960786),
    
17.                     (0.38666666666666666, 0.43137254901960786, 0.23529411764705882),
    
18.                     (0.46, 0.23529411764705882, 0.792156862745098),
    
19.                     (0.7266666666666667, 0.42745098039215684, 1.0),
    
20.                     (0.8533333333333334, 0.0, 0.0),
    
21.                     (1.0, 0.0, 0.0)],
    
22.             'red': [(0.0, 0.0,0.3333333333333333),
    
23.                     (0.12666666666666668, 0.3333333333333333, 0.5294117647058824),
    
24.                     (0.16, 0.5294117647058824, 1.0),
    
25.                     (0.2, 1.0, 0.0),
    
26.                     (0.24, 0.0, 0.6274509803921569),
    
27.                     (0.30666666666666664, 0.6274509803921569, 0.43137254901960786),
    
28.                     (0.38666666666666666, 0.43137254901960786, 0.23529411764705882),
    
29.                     (0.46, 0.23529411764705882, 0.792156862745098),
    
30.                     (0.7266666666666667, 0.42745098039215684, 1.0),
    
31.                     (0.8533333333333334, 0.0, 0.0),
    
32.                     (1.0, 0.0, 0.0)]}

复制代码

BD色阶中有些小数没化为分数，不影响使用，但确实不太美观（
BD色阶同样应用于-100~50度的渲染范围。仔细观察前几行的y0，y1不难发现，同一行的y0，y1不同，代表亮温节点上会出现颜色的跳变，例如CDG在-81跳变到CMG；上一行y1与下一行y0相同，代表亮温节点之间不存在渐变，对应着从-81~-76区间的任意亮温都会被赋予CMG的颜色。

FC-Leo

理论上除了pip install之外，是不是其他数据都可以复制，然后随便改点就行了