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

我不是Carl2

陆地叠图：

伪VIS的陆地叠图有两种方式：Blue Marble和Black Marble。这里我们主要讲解流程更复杂的Blue Marble叠图。

Blue Marble下载网址如下：
https://visibleearth.nasa.gov/im ... d-shaded-topography
下面是个人预处理的Blue Marble底图，海面弄蓝了一点：
https://archive.dapiya.top/AI-VIS_Demo/Worldv3.png
（感谢王站长帮忙生成图片链接）

读取图片可以使用第三方库opencv：

1. import cv2
   
2. image = cv2.imread(r'底图文件路径', cv2.IMREAD_COLOR)

复制代码

Blue Marble是覆盖全球的等经纬投影地图。我们希望将其裁到和卫星亮温、经纬度等数据相同的网格上，即2030×1354的扫描网格，以便后续运算。我们将每个像素的经纬度转化为该经纬度在全球等经纬地图中的行列号，并把对应行列号的Blue Marble像素填回扫描网格。

1. indexlon = np.round(newlon*60).astype(int)
   
2. indexlat = np.rint((90-newlat)*60).astype(int)
   
3. BM = image[indexlat, indexlon, :]/255

复制代码

这时候如果我们直接把Blue Marble画出来：

1. plt.pcolormesh(newlon, newlat, BM)

复制代码

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

1. plt.pcolormesh(newlon, newlat, BM[:, :, ::-1])

复制代码

接下来就是将模拟反射率和底图通过合成算法叠起来。此处许多步骤相当经验性，并且色调与代码可能还有不小的优化空间。有能力的读者可以自己尝试优化这部分代码，不过直接用我现在这版应该也大差不差。
将底图分解为rgb三个通道，并去除一列与合成好的模拟反射率对齐

1. b1, g1, r1 = cv2.split(BM[:, 1:])

复制代码

提亮模拟反射率以免合成后海面过暗

1. b = 0.96 * ref + 25/255

复制代码

初步合成模拟反射率与底图

1. blue = b * b1 + np.power(b, 1.6) * (1 - b1)
   
2. green = b * g1 + np.power(b, 1.6) * (1 - g1)
   
3. red = b * r1 + np.power(b, 1.6) * (1 - r1)
   
4. land_brightness = np.clip((g1 + r1)/2, 0, 0.4)

复制代码

后续调色需要海洋和陆地分开调，通过映射底图的r通道做一个海陆分离，对于陆地像素ocean=0，海面像素ocean=1。

1. norm2 = Normalize(vmin=0, vmax=1, clip=True)
   
2. map2 = ScalarMappable(norm2, LinearSegmentedColormap('1', oceandata))
   
3. ocean = map2.to_rgba(r1)[:, :, 0]

复制代码

用于映射的色阶oceandata：

1. oceandata={'green': [(0, 0, 0),
   
2.                   (1, 1, 1)],
   
3.         'red': [(0, 255/255, 255/255),
   
4.                    (15/255, 255/255, 255/255),
   
5.                    (30/255, 0/255, 0/255),
   
6.                    (1, 0, 0)],
   
7.         'blue': [(0, 0, 0),
   
8.                   (1, 1, 1)]}

复制代码

根据亮温，模拟反射率调整图像饱和度。反射率越高/亮温越低，饱和度越低，陆地的饱和度也设置的比海面更高一些。这里饱和度的定义并不正式，和HSV空间中的有些区别。

1. 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
   
2. avg = (blue + green + red) / 3
   
3. blue = (avg + (blue - avg) * Saturation) * 0.98
   
4. green = avg + (green - avg) * Saturation
   
5. red = avg + (red - avg) * Saturation

复制代码

调整一下洋面颜色（这段两年多前写的有点过于繁琐，大概率是可以继续优化的）：

1. recolor = np.clip(green, 0.1, 0.3) - 0.1
   
2. red = red * (0.2 + recolor * 4) * ocean + (1-ocean) * red
   
3. green = green * (0.8 + recolor) * ocean + (1-ocean) * green
   
4. norm3 = Normalize(vmin=0, vmax=1, clip=True)
   
5. map3 = ScalarMappable(norm3, LinearSegmentedColormap('1', recolordata))
   
6. recolor2 = map3.to_rgba(green)[:, :, 0]
   
7. green = green + (- 10/255 + recolor2) * ocean

复制代码

用于调整洋面颜色的色阶：

1. recolordata={'green': [(0, 0, 0),
   
2.                   (1, 1, 1)],
   
3.         'red':   [(0, 0/255, 0/255),
   
4.                   (8/255, 4/255, 4/255),
   
5.                   (20/255, 2/255, 2/255),
   
6.                   (37/255, 10/255, 10/255),
   
7.                   (57/255, 10/255, 10/255),
   
8.                   (74/255, 16/255, 16/255),
   
9.                   (140/255, 15/255, 15/255),
   
10.                   (160/255, 10/255, 10/255),
    
11.                   (1, 10/255, 10/255)],
    
12.         'blue': [(0, 0, 0),
    
13.                   (1, 1, 1)]}

复制代码

将三个颜色通道合并，限制到0~1的值以便pcolormesh画图：

1. rgb = np.stack([red, green, blue], axis=2)
   
2. rgb = np.clip(rgb, 0, 1)

复制代码

出图效果如下：

FC-Leo

话说有无方法做个类似的把ir-bw云图叠个底图（

我不是Carl2

FC-Leo 发表于 2024-4-29 00:39
话说有无方法做个类似的把ir-bw云图叠个底图（

会比叠伪VIS简单很多。伪VIS叠图的主要困难是陆地反射率比较明显，和海面不一样，IR-BW就没这个问题。IR-BW叠图的话你可以先在Photoshop中试试看哪个叠加模式效果好，然后搜索该叠加模式的算法抄下来就行了。

FC-Leo

1. import matplotlib.pyplot as plt
   
2. from skimage import io
   
3. import math
   
4. import numpy as np
   
5. def Screen(img_1, img_2):
   
6.     img = 1- (1-img_1)*(1-img_2)
   
7. 
   
8.     return img

复制代码

我不是Carl2

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将文件名单独提取出来。

1. import os
   
2. path = r'E:\Data\MOD021KM.A2015296.0515.061.2017322234443.hdf'
   
3. dir, filename = os.path.split(path)
   
4. file = SD(path, SDC.READ)

复制代码

使用第三方库datetime将一年中的日期转化为UTC时间，这里相比实际的时间减了50天以体现季节/气候带相对于太阳赤纬的滞后

1. import datetime
   
2. day = int(filename[14:17])
   
3. time = datetime.datetime(2023, 1, 1) - datetime.timedelta(days=50) + datetime.timedelta(days=day)

复制代码

使用第三方库pyorbital计算此时的太阳赤纬。求季节导致的气候带偏移时，考虑实际的偏移幅度比太阳直射点的偏移幅度小，且主要洋区北半球高海温区域更广阔，这里将太阳赤纬±23.5°的变化范围缩放到了±10度，并向北半球略微偏移了一点。
（顺带一提pyorbital.astronomy还是挺常用的，尤其是不方便直接从数据集取角度信息时，比如AI-VIS的太阳/卫星天顶角/方位角都是用这个读的）

1. from pyorbital.astronomy import sun_ra_dec
   
2. lat_shift = np.rad2deg(sun_ra_dec(time))[1] / 23.5 * 9 + 2

复制代码

根据气候带偏移计算纬度参数，该参数越大对应越冷的区域，注意这里带入的纬度需要是已插值的。

1. 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参考值：

1. SST_guess = 28 - lat_factor * 28

复制代码

用简化后的劈窗算法反演SST。官方文件给出的反演系数对于TERRA与AQUA略有差异，因此需要根据文件名写一个卫星判定。

1. Terra_or_Aqua = (filename[1])
   
2. if Terra_or_Aqua == 'O':
   
3.     SST = 5.77141142 + 0.777202129*bt31 + 0.116606250 *(bt31 - bt32) * SST_guess + 2 #TERRA
   
4. else:
   
5.     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这个参数同样也是为了防止热带云系的反射率偏暗设置的。

1. dif2 = bt32 - bt20
   
2. 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)
   
3. correction = 0.9 + np.clip(lat_factor, 0, 0.5) * 0.2
   
4. 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的亮温差，如果替换到其他波段系数需要作一定调整。

1. cirrus2 = np.clip(dif - 2.2, 0, 2.4)/2.4

复制代码

通过幂运算提高卷云覆盖区域的对比度。由于幂运算会让图像变暗，再乘一个系数线性提亮一下即可。

1. ref_corrected = np.power(ref_warm, 1 + cirrus2) * np.power((1 + cirrus2), 1.35)

复制代码

通过权重参数对高云修正的区域和未经修正的区域做一个过渡。

1. cirrus_weight = np.clip(ref_corrected - 0.735, 0, 0.21)/0.21
   
2. ref_warm = ref_corrected * (1 - cirrus_weight) + ref_warm * cirrus_weight

复制代码

例图就不放了，高云修正在这张Patricia上效果并不明显，主要是渲染VIIRS那种风眼/CDO特写时比较好用。

我不是Carl2

总算在4月的最后几个小时内把教程写完了！明天可能会再来写个简短的后记。

FC-Leo

我不是Carl2 发表于 2024-4-16 08:08
c.        HDF文件读取，MODIS波段简介

卫星云图数据存储的格式相当多样，其中不少格式需要较复杂的脚本才 ...

第一行代码import不了，bug了

我不是Carl2

后记

耗时18天断断续续地写完了这份教程。内容上肯定是不完美的。和刚开始写时预想的一样，教程中的代码有诸多不成熟的地方。一方面缺乏模块化，也没写什么函数和类，导致代码没什么重用价值；两三年前逐步调试陆地叠图、海温反演修正等经验性算法时也遗留下来了形式过于繁杂的公式。不过我想这份有瑕疵的教程多多少少还是有帮助到一些风迷的，不论是画MODIS的部分还是伪VIS的部分，都有人在贴吧上给我“交作业”，让我挺有成就感的。

在编写教程的过程中我自己也有不少收获，首先是搜集资料时对第三方库函数的用途有了更全面的了解，其次优化了一些拖慢运算的步骤（比如以前进行陆地叠图需要渲染几张辅助图片再重新读取），再就是编写教程时的灵感让我时隔10个月再对MODIS伪VIS算法做了两处小更新，分别减轻了扫描边缘过度的光影效果以及解决了云层透光问题。

后续有什么需要勘误的地方我会及时跟进，最后再感谢一下当初刚开始做MODIS处理和伪VIS时给过我很多指导的HCl（py）吧，也谢谢各位一直以来对伪VIS项目和别的云图存档的支持！

平顶山气象萌新

注册账户用QQ邮箱是不行的吗，我只有qq邮箱啊别搞

FC-Leo

平顶山气象萌新 发表于 2024-5-5 21:45
注册账户用QQ邮箱是不行的吗，我只有qq邮箱啊别搞

应该可以，不过那个机器人认证必须到国外才能搞吧