用MODIS L1B数据绘制真彩色云图

灭火器

Aqua 卫星的 MYD021KM 产品 36 个波段的 1 km 分辨率的地表反射率数据，将其中波段 1（0.65 μm）的数值作为红色分量，波段 4（0.56 μm）的数值作为绿色分量，波段 3（0.47 μm）的数值作为蓝色分量，即可合成与人眼感受类似的彩色云图。具体来说需要以下步骤：
（1）从 MYD021KM 配套的 MYD03 产品中读取 1 km 分辨率的地理经纬度数据。
（2）根据反射率变量的元数据将 short 类型的反射率缩放为浮点型。
（3）将三个波段的反射率合并为 RGB 数组，处理缺测和超出 [0, 1] 范围的元素。
（4）将 RGB 数组用 pcolormesh 方法画在地图上，每个四边形像元的颜色通过 color 参数给出。
（5）增强图像的亮度，比较增强前后的效果。

MYD021KM 数据的变量说明和缩放方法请见官方文档：
https://mcst.gsfc.nasa.gov/l1b-documents
https://ftp.ssec.wisc.edu/pub/willemm/Creating_Reprojected_True_Color_MODIS_Images_A_Tutorial_process.pdf
增强图像亮度采用了两种方法：
（1）伽马矫正（http://unidata.github.io/python-gallery/examples/mapping_GOES16_TrueColor.html）
（2）分段映射（http://www.idlcoyote.com/ip_tips/brightmodis.html）
pcolormesh 方法画 RGB 数组的颜色参考了
https://stackoverflow.com/questions/41389335/how-to-plot-geolocated-rgb-data-faster-using-python-basemap
https://github.com/matplotlib/matplotlib/issues/4277

读取文件使用 pyhdf 包，画图使用 cartopy 包。代码如下：

1. '''
   
2. @author : 灭火器
   
3. @time : 2022-04-15
   
4. @use : 利用MYD02产品画出真彩色图, 并比较不同图像处理方法的效果.
   
5. '''
   
6. from pathlib import Path
   
7. from datetime import datetime
   
8. 
   
9. from pyhdf.SD import SD, SDC
   
10. import numpy as np
    
11. from scipy.interpolate import interp1d
    
12. import matplotlib.pyplot as plt
    
13. import matplotlib.ticker as mticker
    
14. import cartopy.crs as ccrs
    
15. import cartopy.feature as cfeature
    
16. from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
    
17. 
    
18. def to_reflectance(sds):
    
19.     '''将反射率的SDS变为浮点型数组.'''
    
20.     ref = sds[:]
    
21.     attrs = sds.attributes()
    
22.     scales = np.array(attrs['reflectance_scales'])
    
23.     offsets = np.array(attrs['reflectance_offsets'])
    
24.     index = np.s_[:, np.newaxis, np.newaxis]
    
25.     ref = np.where(
    
26.         ref > 32767, np.nan,
    
27.         (ref - offsets[index]) * scales[index]
    
28.     )
    
29. 
    
30.     return ref
    
31. 
    
32. def scale_modis(rgb):
    
33.     '''增亮MODIS的彩色图.'''
    
34.     x = np.array([0, 30, 60, 120, 190, 255]) / 255
    
35.     y = np.array([0, 110, 160, 210, 240, 255]) / 255
    
36.     f = interp1d(x, y, bounds_error=False, fill_value=1)
    
37. 
    
38.     return f(rgb)
    
39. 
    
40. if __name__ == '__main__':
    
41.     # 输入和输出目录.
    
42.     dirpath_data = Path('./data')
    
43.     dirpath_pics = Path('./pics')
    
44.     if not dirpath_pics.exists():
    
45.         dirpath_pics.mkdir()
    
46. 
    
47.     # 输入文件.
    
48.     filepath_MYD02 = dirpath_data / 'MYD021KM.A2018120.0500.061.2018120164637.hdf'
    
49.     filepath_MYD03 = dirpath_data / 'MYD03.A2018120.0500.061.2018120162824.hdf'
    
50.     # 提取文件名中的时间.
    
51.     parts = filepath_MYD02.name.split('.')
    
52.     date = datetime.strptime(parts[1] + parts[2], 'A%Y%j%H%M')
    
53. 
    
54.     # 读取1km经纬度.
    
55.     sd = SD(str(filepath_MYD03), SDC.READ)
    
56.     lon = sd.select('Longitude')[:]
    
57.     lat = sd.select('Latitude')[:]
    
58.     sd.end()
    
59.     # 地图范围由granule的范围决定.
    
60.     extents = [lon.min(), lon.max(), lat.min(), lat.max()]
    
61. 
    
62.     # 读取反射率.
    
63.     sd = SD(str(filepath_MYD02), SDC.READ)
    
64.     ref250 = to_reflectance(sd.select('EV_250_Aggr1km_RefSB'))
    
65.     ref500 = to_reflectance(sd.select('EV_500_Aggr1km_RefSB'))
    
66.     sd.end()
    
67.     ref1 = ref250[0, :, :]
    
68.     ref3 = ref500[0, :, :]
    
69.     ref4 = ref500[1, :, :]
    
70. 
    
71.     # 合并RGB数组, 处理缺测部分.
    
72.     rgb = np.dstack((ref1, ref4, ref3))
    
73.     rgb[np.any(np.isnan(rgb), axis=-1)] = 1
    
74.     rgb = np.clip(rgb, 0, 1)
    
75. 
    
76.     # 三种处理的图像.
    
77.     rgb1 = rgb
    
78.     rgb2 = np.power(rgb, 1 / 2.2)
    
79.     rgb3 = scale_modis(rgb)
    
80. 
    
81.     # 组图形状为(1, 3).
    
82.     proj = ccrs.PlateCarree()
    
83.     subplot_kw = {'projection': proj}
    
84.     fig, axes = plt.subplots(1, 3, figsize=(15, 5), subplot_kw=subplot_kw)
    
85.     fig.subplots_adjust(wspace=0.2)
    
86.     # 画出地图.
    
87.     for ax in axes:
    
88.         ax.coastlines(resolution='50m', lw=0.5)
    
89.         ax.add_feature(cfeature.BORDERS, lw=0.5)
    
90.         ax.set_xticks(np.arange(-180, 181, 10), crs=proj)
    
91.         ax.set_yticks(np.arange(-90, 91, 10), crs=proj)
    
92.         ax.xaxis.set_minor_locator(mticker.AutoMinorLocator(2))
    
93.         ax.yaxis.set_minor_locator(mticker.AutoMinorLocator(2))
    
94.         ax.xaxis.set_major_formatter(LongitudeFormatter())
    
95.         ax.yaxis.set_major_formatter(LatitudeFormatter())
    
96.         ax.set_extent(extents, crs=proj)
    
97.         ax.tick_params(labelsize='small')
    
98. 
    
99.     # 画出三种图像.
    
100.     rgbs = [rgb1, rgb2, rgb3]
     
101.     titles = ['Original', 'Gamma Correction', 'Enhancement']
     
102.     for i, ax in enumerate(axes):
     
103.         colors = rgbs[i][:-1, :-1].reshape(-1, 3)
     
104.         ax.pcolormesh(
     
105.             lon, lat, lon[:-1, :-1], color=colors,
     
106.             shading='flat', transform=proj
     
107.         )
     
108.         ax.set_title(titles[i], fontsize='medium')
     
109. 
     
110.     # 设置图片标题.
     
111.     t = date.strftime('%Y-%m-%d %H:%M')
     
112.     fig.suptitle(t, y=0.9, fontsize='large')
     
113. 
     
114.     # 保存图片.
     
115.     filepath_output = dirpath_pics / 'true_color.png'
     
116.     fig.savefig(str(filepath_output), dpi=300, bbox_inches='tight')
     
117.     plt.close(fig)

复制代码

效果如下图，可以看到如果不增强亮度的话，无云的陆面和洋面会显得黑糊糊的，难以识别；而增强后云又会显得过亮，损失了一些纹理细节。


下图是 https://worldview.earthdata.nasa.gov/ 网站的在线图片，与上图的差别是陆面植被的绿色更加醒目，可能是官方对反射率进行了其它订正操作。但总的来说效果还是比较接近的。


代码中使用 pcolormesh 画图是为了保留 MODIS granule 非规则的网格结构。本来 pcolor 和 pcolormesh 只能画伪彩色图，默认使用 viridis 色标填色，这里通过 color 参数用 RGB 数组的颜色覆盖了默认的填色。其实 Axes 的 pcolorfast 方法能直接画 RGB(A) 数组，但因为这个方法仍处于实验状态，且 GeoAxes 并没有对其进行改写，所以不保证能在别的地图投影下正常工作。相比之下 Meteoinfo 直接就有 imshowm 函数……

对于大小为 (2030, 1354) 的数组来说，这个程序每张图大概要画上十几秒，效率堪忧。这里给出另一种处理思路：利用 pyresample 包的 resample_nearest 函数对 swath 数据进行最邻近插值，得到网格化的 RGB 数组，然后通过 imshow 方法直接画在地图上。计算时间会稍微延长，但画图时间会快十倍。不过 pyresample 包对区域的定义我还看不懂，所以这里就只提供简单的思路。

我在网上搜到的教程大多是画同步卫星的真彩色图（使用 Satpy 包），而关于 MODIS 产品的教程则非常少，所以这里的程序也只是一种尝试，也可能存在错误或需要改进的地方，还请熟悉的朋友多多交流。

我是sswwkk

大神 怎么修改绘图范围呢

灭火器

我是sswwkk 发表于 2022-4-28 10:48
大神 怎么修改绘图范围呢

你是指什么，地图范围还是云图范围

忘川思源

请问楼主为什么我会出现这个情况呢？应该如何解决

坎坷

imshow应该是最快的，特别对云图这种高分辨率数据，其他绘图方式都太慢了

灭火器

忘川思源 发表于 2022-7-2 11:26
请问楼主为什么我会出现这个情况呢？应该如何解决

看不太懂是哪里出了问题。另外建议阅读这篇：公众号
帖子中的方法非常慢

忘川思源

灭火器 发表于 2022-7-4 13:27
看不太懂是哪里出了问题。另外建议阅读这篇：公众号
帖子中的方法非常慢

好嘞 谢谢大佬 我现在也开始学您尝试处理GPM资料

Franny

lon, lat, lon[:-1, :-1], color=colors
请问pcolormesh的第三个变量是怎么弄的？我自己没理解