葵花八号卫星AHI云图绘制
前言
需求：大家看到诸多文献使用卫星云图作天气形势系统分析时，想必也想自己也为文章中加一张
面向群体：需要使用卫星云图进行天气学分析或天气系统阐释的小伙伴，当然你喜欢卫星云图做壁纸也可以画着玩
应用场景：汇报or写作
himawari 8是什么
隶属于日本气象厅的（JMA）Himawari-8 卫星于 2014 年 10 月 7 日从种子岛航天中心发射升空，同年 10 月 16 日正式开始在地球同步轨道上运行，
是世界气象组织（WMO）主导的全球观测系统（Global Observing System）的重要组成部分，可为西太平洋的 30 多个国家和地区提供气象服务。
Himawari-8 为三轴姿态控制的静止卫星，运行倾角 0.03 度，搭载了 AdvancedHimawari Imager（AHI）和
Space Environment Data Acquisition Monitor(SDAM)两个主要传感器,分别对地球和宇宙空间进行观测，
其同步轨道位于东经 140.7 度，高度约 35786km。
AHI是什么
Advanced Himawari Imager（AHI）是搭载在 Himawari-8 卫星上对地球系统进行观测的仪器，是一种多通道光谱成像仪，用于捕捉亚太地区的可见光和红外图像。
AHI 传感器在某些波段上的空间分辨率高达 500m，大部分波段的空间分辨率为 2km,
可以提供时间分辨率为 10 分钟的全盘观测图像和时间分辨率为 2.5分钟的日本地区观测图像。
RGB是什么
影像数据根据波段的多少可以分为单波段影像和多波段影像两种，单波段影像一般用黑白色的灰度图来描述，多波段常用RGB 合成象素值的彩色图来描述，就是将三个波段的数据分别通过红、绿、蓝三个通道加载，从而渲染出。
将多波段影像数据添加到地图中之后，可使用多波段栅格数据集中的任意三个可用波段的组合来创建 RGB 合成图。与仅处理一个波段相比，通过将多个波段共同显示为RGB 合成图通常可从数据集收集到更多信息。
下面开始绘图
import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import numpy as np
import metpy
import glob
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
files = glob.glob('/home/mw/input/himawari7828/0800_20160814080000-P1S-ABOM_OBS_B03-PRJ_GEOS141_2000-HIMAWARI8-AHI/*.nc')
ds2 = xr.open_mfdataset(files)
b1 = ds2['channel_0001_scaled_radiance']
b2 = ds2['channel_0002_scaled_radiance']
b3 = ds2['channel_0003_scaled_radiance']
x = ds2.variables['x'][:]
y = ds2.variables['y'][:]
plt.figure(figsize=(6,6))
plt.imshow(b1, extent=, x[-1], y[-1], y])<matplotlib.image.AxesImage at 0x7f41442ff650>
三色可视化
RGB: red green blue
In :fig, () = plt.subplots(1, 3, figsize=(16, 3))
ax1.imshow(b2, cmap='Reds', vmax=1, vmin=0)
ax1.set_title('Red', fontweight='bold')
ax1.axis('off')
ax2.imshow(b3, cmap='Greens', vmax=1, vmin=0)
ax2.set_title('Veggie', fontweight='bold')
ax2.axis('off')
ax3.imshow(b1, cmap='Blues', vmax=1, vmin=0)
ax3.set_title('Blue', fontweight='bold')
ax3.axis('off')
(-0.5, 5499.5, 5499.5, -0.5)
RGB真彩填色
大家可以试试用不同的方法设置RGB，达到自己想要的效果
#此处为第一种方法
##提取3个波段为RGB
rgb = np.stack((
b1.to_numpy(),
b2.to_numpy(),
b3.to_numpy(),
), axis=-1)
#扩充RGB值范围使对比度增强
rgb = np.clip((rgb/0.7)**0.7, 0, 1)
#此处注释的代码是第二种RGB图像增强方法
#vmin = 0
#vmax = 0.5
#B1 = b1.clip(vmin, vmax) / vmax * 255
#B2 = b2.clip(vmin, vmax) / vmax * 255
#B3 = b3.clip(vmin, vmax) / vmax * 255
#rgb = np.stack((B3, B2, B1), axis=2).astype('uint8')
# Plot image
plt.figure(figsize=(12,12))
plt.imshow(rgb, extent=, x[-1], y[-1], y])
# Add labels
plt.title('Himawari 8 \n', fontsize=20)
# Adjust tick mark size
plt.tick_params(labelsize=16)
plt.show()/opt/conda/lib/python3.7/site-packages/ipykernel_launcher.py:17: RuntimeWarning: invalid value encountered in power
两种对比度增强的方式有各自的优缺点。
第一种方式（第一段代码）使用了非线性的幂函数转换来增强对比度，然后通过np.clip函数将值限制在0和1之间。优点如下：
对比度增强更强烈：幂函数转换可以显著增强图像的对比度，使得细节更加突出。
缺点如下：
数据类型变化：最终生成的RGB图像是浮点类型的数组，相比于无符号整数类型，需要更多的内存存储和计算资源。
第二种方式（第二段代码）使用了线性缩放的方式进行对比度增强，并将值限制在0和255之间。优点如下：
数据类型保持不变：最终生成的RGB图像是无符号整数类型(uint8)的数组，节省了内存存储和计算资源。
缺点如下：
对比度增强相对较弱：线性缩放对比度增强的效果相对较弱，可能无法显著改善图像细节的可见性。
综上所述，第一种方式以非线性幂函数转换为基础的对比度增强效果更强烈，但需要更多的计算资源和内存存储；
而第二种方式以线性缩放为基础的对比度增强效果较弱，但在数据类型保持不变的情况下，节省了计算资源和内存存储。
选择适合的方式取决于具体的需求和应用场景
叠加地图
metpy是非常强大的气象python库，在前几篇推文也有它的活跃场景。
#获取投影
dat = ds2.metpy.parse_cf('channel_0001_scaled_radiance')
geos = dat.metpy.cartopy_crs
# We also need the x (north/south) and y (east/west) axis sweep of the ABI data
x = dat.x
y = dat.y
fig = plt.figure(figsize=(15, 12))
# Create axis with Geostationary projection
ax = fig.add_subplot(1, 1, 1, projection=geos)
# Add the RGB image to the figure. The data is in the same projection as the
# axis we just created.
ax.imshow(rgb, origin='upper',extent=(x.min(), x.max(), y.min(), y.max()), transform=geos)
# Add Coastlines and States
ax.coastlines(resolution='50m', color='white', linewidth=1)
ax.add_feature(ccrs.cartopy.feature.BORDERS, linewidth=1)
plt.title('True Color', loc='left', fontweight='bold', fontsize=15)
plt.show()
显示局部地区并叠加地图
#获取投影
dat = ds2.metpy.parse_cf('channel_0001_scaled_radiance')
geos = dat.metpy.cartopy_crs
#转换RGB
vmin = 0
vmax = .5
B1 = b1.clip(vmin, vmax) / vmax * 255
B2 = b2.clip(vmin, vmax) / vmax * 255
B3 = b3.clip(vmin, vmax) / vmax * 255
X = x
Y = y
rgb = np.stack((B3, B2, B1), axis=2).astype('uint8')
fig = plt.figure(figsize=(20,12))
ax1 = fig.add_subplot(121, projection=geos)
ax2 = fig.add_subplot(122, projection=geos)
# Plot image
ax1.imshow(rgb,extent=(X, X[-1], Y[-1], Y))
ax2.imshow(rgb, origin='upper',extent=(X, X[-1], Y[-1], Y), transform=geos)
# Add Coastlines and States
ax2.coastlines(resolution='50m', color='white', linewidth=1)
ax2.add_feature(ccrs.cartopy.feature.STATES, linewidth=1)
# Add labels
plt.title('Himawari 8 \n', fontsize=20)
# Adjust tick mark size
plt.tick_params(labelsize=16)
本文参考了GOES-16: True Color Recipe — Unidata Python Gallery 完整文件与代码可后台回复葵花获取
文章来源于公众号：气ython风雨
页:
[1]