合成低亮温区模拟反射率：
先进行亮温差分
dif = bt31 - bt32
我们希望把差分亮温通过一定的映射转化为模拟反射率。3.a介绍的色阶本质上就是一类映射，而这里我们所需的映射也可以通过色阶达成。pcolormesh+色阶的办法只能给出色阶映射后的图像，而我们想要的是映射后的二维数组以方便后续计算，这就需要matplotlib.colors.Normalize和matplotlib.cm.ScalarMappable这两个功能。
from matplotlib.colors import LinearSegmentedColormap, Normalize
from matplotlib.cm import ScalarMappable
定义映射方式。色阶范围为-15~15度的亮温差分值，色阶名为colddata。
norm1 = Normalize(vmin=-15, vmax=15, clip=True)
map1 = ScalarMappable(norm1, LinearSegmentedColormap('1', colddata))
colddata是用于强调高反射率云盖与低反射率卷云的色阶，在映射过程中拉高了两者的对比度，有助于分辨CDO边缘。色阶的具体数值是经验性的，各位也可以按照个人喜好微调。做从亮温差分到模拟反射率的映射中，并不需要多种颜色，因此这里rgb三色中只有red起作用。
colddata={'green': [(0, 0, 0),
(1, 1, 1)],
'red': [(0, 1, 1),
(0.37, 1, 1),
(0.438, 0.98, 0.98),
(0.468, 0.94, 0.94),
(0.498, 0.85, 0.85),
(0.508, 0.81, 0.81),
(0.538, 0.75, 0.75),
(0.638, 0.61, 0.61),
(0.758, 0.51, 0.51),
(1, 0.5, 0.5)],
'blue': [(0, 0, 0),
(1, 1, 1)]}
映射后会得到rgba形式的低温区模拟反射率，是2030×1354×4的三维数组，其中第三个维度分别是r, g, b, a，用索引取出r就得到了低温区的模拟反射率。
ref_cold = map1.to_rgba(dif)[:, :, 0]

对高温区与低温区的模拟反射率做加权合成：
先对高温区模拟反射率进行一下预处理，去除Band20低温区信噪比不够造成的nan值，并将小于0的值全部设为0，大于1的值全部设为1。
ref_warm = np.nan_to_num(ref_warm, nan=1)
ref_warm = np.clip(ref_warm, 0, 1)
低温区反射率由于是用ScalarMappable映射的，已经处在0~1的区间了，因此无需进一步处理。
通过Band31亮温设置两份反射率的权重并合成。这里本质上也是一个从Band31亮温到权重参数的映射，但由于是线性的，不需要ScalarMappable创建映射，用np.clip也可以达到相同效果。
weight = (np.clip(bt31, -70, -20) + 70)/50
ref = ref_warm * weight + ref_cold * (1 - weight)

合成光影效果：
用数据切片将Band31亮温最左侧/最右侧一列数据去除并相减，达到“错位差分”的效果。[:, 0:-1]与[:, 1:]可以互换以达到模拟相反方向光照。
shade = bt31[:, 0:-1] - bt31[:, 1:]
随后将错位差分的亮温值限制在-4~4的范围内并乘上系数0.08，这意味着光影效果最多改变±0.32的反射率。这些数值可以按照个人喜好与不同卫星的分辨率调整，对于分辨率较高的卫星，错位差分的范围可以设置的小一些，系数则可以设置的大一些。
shade = np.clip(shade, -4, 4) * 0.08

（上图为shade的图像，渲染范围改到了-0.5~0.5，并把经纬度也去除了一列以对齐数据）
图像中不是所有区域都应该叠加上如图所示的光影效果，本身反射率偏低的区域应有的光影效果会弱一些，卷云区域和其他非深对流云顶呈现的光影也应稍弱。这里我们定义一个“光影参数”，由先前合成的模拟反射率，Band31-32差分做的卷云识别，以及Band31亮温对深对流的识别共同决定。
cirrus = (np.clip(dif, 1, 3) - 1) / 2
deep_conv = (np.clip(bt31, -70, -50) + 70) / 20
coeff = np.clip(ref - cirrus, 0, 1) * (1 - 0.5 * deep_conv)
将模拟反射率叠上光影效果，注意光影效果的数组少一列，因此反射率和光影系数也要裁掉一列才能运算。
ref = ref[:, 1:] + shade * coeff[:, 1:]

注入灵魂.png

陆地叠图：
伪VIS的陆地叠图有两种方式：Blue Marble和Black Marble。这里我们主要讲解流程更复杂的Blue Marble叠图。
Blue Marble：
https://visibleearth.nasa.gov/images/57752/blue-marble-land-surface-shallow-water-and-shaded-topography
下面是个人预处理的Blue Marble底图，海面弄蓝了一点：
archive.dapiya.top/AI-VIS_Demo/Worldv3.png
（感谢@小谦.WX 帮忙生成图片链接）
读取图片可以使用第三方库opencv：
import cv2
image = cv2.imread(r'底图文件路径', cv2.IMREAD_COLOR)
Blue Marble是覆盖全球的等经纬投影地图。我们希望将其裁到和卫星亮温、经纬度等数据相同的网格上，即2030×1354的扫描网格，以便后续运算。我们将每个像素的经纬度转化为该经纬度在全球等经纬地图中的行列号，并把对应行列号的Blue Marble像素填回扫描网格。
indexlon = np.round(lon*60).astype(int)
indexlat = np.rint((90-lat)*60).astype(int)
BM = image[indexlat, indexlon, :]/255
这时候如果我们直接把Blue Marble画出来：
plt.pcolormesh(newlon, newlat, BM)

陆地的裁取没什么问题，且底图确实和亮温，经纬度等数据处于同一网格。唯一的问题是颜色不对，这是由于Matplotlib默认按照r, g, b的顺序，而opencv读取的顺序是b, g, r，将BM的第三个维度取倒序就能得到颜色正确的底图了。
plt.pcolormesh(newlon, newlat, BM[:, :, ::-1])

接下来就是将模拟反射率和底图通过合成算法叠起来。此处许多步骤相当经验性，并且色调与代码可能还有不小的优化空间。有能力的读者可以自己尝试优化这部分代码，不过直接用我现在这版应该也大差不差。
将底图分解为rgb三个通道，并去除一列与合成好的模拟反射率对齐
b1, g1, r1 = cv2.split(BM[:, 1:])
提亮模拟反射率以免合成后海面过暗
b = 0.96 * ref + 25/255
初步合成模拟反射率与底图
blue = b * b1 + np.power(b, 1.6) * (1 - b1)
green = b * g1 + np.power(b, 1.6) * (1 - g1)
red = b * r1 + np.power(b, 1.6) * (1 - r1)
land_brightness = np.clip((g1 + r1)/2, 0, 0.4)
后续调色需要海洋和陆地分开调，通过映射底图的r通道做一个海陆分离，对于陆地像素ocean=0，海面像素ocean=1。
norm2 = Normalize(vmin=0, vmax=1, clip=True)
map2 = ScalarMappable(norm2, LinearSegmentedColormap('1', oceandata))
ocean = map2.to_rgba(r1)[:, :, 0]
用于映射的色阶oceandata：
oceandata={'green': [(0, 0, 0),
(1, 1, 1)],
'red': [(0, 255/255, 255/255),
(15/255, 255/255, 255/255),
(30/255, 0/255, 0/255),
(1, 0, 0)],
'blue': [(0, 0, 0),
(1, 1, 1)]}
根据亮温，模拟反射率调整图像饱和度。反射率越高/亮温越低，饱和度越低，陆地的饱和度也设置的比海面更高一些。这里饱和度的定义并不正式，和HSV空间中的有些区别。
Saturation = 1.5 * (1 - np.clip(ref - 0.4, 0, 0.5)/0.5) + (1 - ocean) * 3 * (np.clip(bt31[:, 1:], -60, 10) + 60)/70
avg = (blue + green + red) / 3
blue = (avg + (blue - avg) * Saturation) * 0.98
green = avg + (green - avg) * Saturation
red = avg + (red - avg) * Saturation
调整一下洋面颜色（这段两年多前写的有点过于繁琐，大概率是可以继续优化的）：
recolor = np.clip(green, 0.1, 0.3) - 0.1
red = red * (0.2 + recolor * 4) * ocean + (1-ocean) * red
green = green * (0.8 + recolor) * ocean + (1-ocean) * green
norm3 = Normalize(vmin=0, vmax=1, clip=True)
map3 = ScalarMappable(norm3, LinearSegmentedColormap('1', recolordata))
recolor2 = map3.to_rgba(green)[:, :, 0]
green = green + (- 10/255 + recolor2) * ocean
用于调整洋面颜色的色阶：
recolordata={'green': [(0, 0, 0),
(1, 1, 1)],
'red': [(0, 0/255, 0/255),
(8/255, 4/255, 4/255),
(20/255, 2/255, 2/255),
(37/255, 10/255, 10/255),
(57/255, 10/255, 10/255),
(74/255, 16/255, 16/255),
(140/255, 15/255, 15/255),
(160/255, 10/255, 10/255),
(1, 10/255, 10/255)],
'blue': [(0, 0, 0),
(1, 1, 1)]}
将三个颜色通道合并，限制到0~1的值以便pcolormesh画图：
rgb = np.stack([red, green, blue], axis=2)
rgb = np.clip(rgb, 0, 1)
出图效果如下：

教程还剩最后一节没写，这个月也只剩一天半了
怎么我做啥都得赶着ddl交啊

b. 海温反演修正与高云修正
相对于前面介绍的伪VIS基础算法，海温反演修正和高云修正更像是锦上添花的步骤，并且更具实验性也更不成熟。但既然目录中列都列出来了，还是把相关代码贴一下供参考。
目前版本的两项修正互相不兼容，一颗卫星只能应用一项修正。MODIS和H8伪VIS应用的是海温反演修正以改善高纬度的成像效果，而VIIRS，SLSTR等卫星一般绘制的范围较小，因此应用的是高云修正以改善TC卷云下方的细节。
海温反演修正调整了高温区模拟反射率（ref_warm）算法，伪VIS的其他部分则保持不变。海温反演公式中除了两个劈窗算法用到的波段Band31和Band32，还需要输入一个“SST参考值” Tsfc。我们这里不作精确的海温反演，仅用于修正伪VIS，这个参考值用季节+纬度估算一下即可。

文件名中包含的一年中的日期可以用于估算季节，例如MOD021KM.A2015296.0515.061.2017322234443.hdf这个文件就是2015年第296天拍的。我们稍微修改一下读取HDF时路径的输入方式即可用第三方库os将文件名单独提取出来。
import os
path = r'E:\Data\MOD021KM.A2015296.0515.061.2017322234443.hdf'
dir, filename = os.path.split(path)
file = SD(path, SDC.READ)
使用第三方库datetime将一年中的日期转化为UTC时间，这里相比实际的时间减了50天以体现季节/气候带相对于太阳赤纬的滞后
import datetime
day = int(filename[14:17])
time = datetime.datetime(2023, 1, 1) - datetime.timedelta(days=50) + datetime.timedelta(days=day)
使用第三方库pyorbital计算此时的太阳赤纬。求季节导致的气候带偏移时，考虑实际的偏移幅度比太阳直射点的偏移幅度小，且主要洋区北半球高海温区域更广阔，这里将太阳赤纬±23.5°的变化范围缩放到了±10度，并向北半球略微偏移了一点。
（顺带一提pyorbital.astronomy还是挺常用的，尤其是不方便直接从数据集取角度信息时，比如AI-VIS的太阳/卫星天顶角/方位角都是用这个读的）
from pyorbital.astronomy import sun_ra_dec
lat_shift = np.rad2deg(sun_ra_dec(time))[1] / 23.5 * 9 + 2
根据气候带偏移计算纬度参数，该参数越大对应越冷的区域，注意这里带入的纬度需要是已插值的。
lat_factor = np.clip(np.absolute(newlat - lat_shift) - 15, 0, 45)/50 * 1.22 * (1.06 - lat_shift/100) + 0.01 - 9/250
给出（非常）粗略的SST参考值：
SST_guess = 28 - lat_factor * 28
用简化后的劈窗算法反演SST。官方文件给出的反演系数对于TERRA与AQUA略有差异，因此需要根据文件名写一个卫星判定。
Terra_or_Aqua = (filename[1])
if Terra_or_Aqua == 'O':
SST = 5.77141142 + 0.777202129*bt31 + 0.116606250 *(bt31 - bt32) * SST_guess + 2 #TERRA
else:
SST = 4.98995447 + 0.806239545*bt31 + 0.106239545 * (bt31 - bt32) * SST_guess + 2 #AQUA
劈窗算法的SST对于有云遮盖的区域严重偏低，直接带入高温区模拟反射率的计算公式会导致部分云系明显偏暗。这里我们用Band32和Band20的差分将这些云识别出来，化为一个权重（SST_weight），对于有云的区域该权重偏低，倾向于使用原本的公式(bt20 - 28) / (bt31 - 30)，而对于无云区域该权重接近1，会采用带入反演SST的方式计算反射率。correction这个参数同样也是为了防止热带云系的反射率偏暗设置的。
dif2 = bt32 - bt20
SST_weight = np.clip(dif2 + 13, 0, 13) / 13 * (0.6 + 1.1 * lat_factor) * (0.3 + np.clip(SST + 20, 0, 20)/20 * 0.7)
correction = 0.9 + np.clip(lat_factor, 0, 0.5) * 0.2
ref_warm = (bt20 - (SST-30) * correction - 28)/(bt31 - (SST - 30) * correction - 30) * SST_weight + (bt20 - 28) / (bt31 - 30) * (1 - SST_weight)
（当时公式写这么长不分段也不写模块化真的是很糟糕的习惯，大家不要学我）
这里得到的ref_warm即可替代原有的高温区反射率。海温修正前与修正后效果对比如下：


看上去效果差不多，这是因为Patricia本身海温就很高，不需要什么海温修正，甚至原本因为高海温过暗的洋面还稍微变亮了一点。如果画一下5分钟后的下一片MODIS数据对比就很明显了：


海温修正前的海陆颜色均有明显偏白的现象，修正后基本正常了。部分陆地云系和冷流云的模拟还是不够准确，不过伪VIS作为主打热带气旋的算法，高纬区域效果做到大致能看即可。

高云修正同样是调整了高温区模拟反射率（ref_warm）的算法。
首先将卷云识别出来，这里的系数对于不同的长波波段组合会不同，以下dif代表VIIRS M15 M16的亮温差，如果替换到其他波段系数需要作一定调整。
cirrus2 = np.clip(dif - 2.2, 0, 2.4)/2.4
通过幂运算提高卷云覆盖区域的对比度。由于幂运算会让图像变暗，再乘一个系数线性提亮一下即可。
ref_corrected = np.power(ref_warm, 1 + cirrus2) * np.power((1 + cirrus2), 1.35)
通过cirrus_weight对高云修正的区域和未经修正的区域做一个过渡。
cirrus_weight = np.clip(ref_corrected - 0.735, 0, 0.21)/0.21
ref_warm = ref_corrected * (1 - cirrus_weight) + ref_warm * cirrus_weight
例图就不放了，高云修正在这张Patricia上效果并不明显，主要是渲染VIIRS那种风眼/CDO特写时比较好用。

总算在四月最后几个小时把教程写完了！不过最后一小段看起来还在被吞
明天可能会再来写个后记之类的

后记
耗时18天断断续续地写完了这份教程。内容上肯定是不完美的。和刚开始写时预想的一样，教程中的代码有诸多不成熟的地方。一方面缺乏模块化，也没写什么函数和类，导致代码没什么重用价值；两三年前逐步调试陆地叠图、海温反演修正等经验性算法时也遗留下来了形式过于繁杂的公式。不过我想这份有瑕疵的教程多多少少还是有帮助到一些风迷的，不论是画MODIS的部分还是伪VIS的部分，都有人在贴吧上给我“交作业”，让我挺有成就感的。
在编写教程的过程中我自己也有不少收获，首先是搜集资料时对第三方库函数的用途有了更全面的了解，其次优化了一些拖慢运算的步骤（比如以前进行陆地叠图需要渲染几张辅助图片再重新读取），再就是编写教程时的灵感让我时隔10个月再对MODIS伪VIS算法做了两处小更新，分别减轻了扫描边缘过度的光影效果以及解决了云层透光问题。
后续有什么需要勘误的地方我会及时跟进，最后再感谢一下当初刚开始做MODIS处理和伪VIS时给过我很多指导的HCl（py）吧，也谢谢各位一直以来对伪VIS项目和别的云图存档的支持！