葵花八号AHI真彩图可视化

用户11172986

发布于 2024-07-18 13:47:40

840

发布于 2024-07-18 13:47:40

文章被收录于专栏：气python风雨

葵花八号真彩图绘制

🎓 作者：酷炫用户名

前言

需求：大家看到诸多文献使用卫星云图作为天气形势系统介绍时想必也想自己也为文章中加一张 面向群体：需要使用卫星云图进行天气学分析或天气系统阐释的小伙伴，当然你喜欢卫星云图做壁纸也可以画着玩 应用场景：汇报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分钟的日本地区观测图像。

以上文字摘自[这里]

(www.satdatafresh.com/readme/H8%E5%8F%B0%E9%A3%8Ereadme%20file3.0.pdf)

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'][0,:,:]
b2 = ds2['channel_0002_scaled_radiance'][0,:,:]
b3 = ds2['channel_0003_scaled_radiance'][0,:,:]

x = ds2.variables['x'][:]
y = ds2.variables['y'][:]

plt.figure(figsize=(6,6))
plt.imshow(b1, extent=[x[0], x[-1], y[-1], y[0]])
plt.show()

三色可视化

RGB: red green blue

fig, ([ax1, ax2, ax3]) = 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')
plt.show()

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[0], x[-1], y[-1], y[0]])
# Add labels
plt.title('Himawari 8 \n', fontsize=20)

# Adjust tick mark size
plt.tick_params(labelsize=16)
plt.show()

/tmp/ipykernel_58/2589970593.py:10: RuntimeWarning: invalid value encountered in power
  rgb = np.clip((rgb/0.7)**0.7, 0, 1)
/opt/conda/lib/python3.9/site-packages/matplotlib/cm.py:494: RuntimeWarning: invalid value encountered in cast
  xx = (xx * 255).astype(np.uint8)

两种对比度增强的方式有各自的优缺点。

第一种方式（第一段代码）使用了非线性的幂函数转换来增强对比度，然后通过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()

/opt/conda/lib/python3.9/site-packages/matplotlib/cm.py:494: RuntimeWarning: invalid value encountered in cast
  xx = (xx * 255).astype(np.uint8)

增加经纬度刻度

import cartopy.crs as ccrs
import matplotlib.pyplot as plt

dat = ds2.metpy.parse_cf('channel_0001_scaled_radiance')
geos = dat.metpy.cartopy_crs
x = dat.x
y = dat.y

fig = plt.figure(figsize=(15, 12))
ax = fig.add_subplot(1, 1, 1, projection=geos)
ax.imshow(rgb, origin='upper', extent=(x.min(), x.max(), y.min(), y.max()), transform=geos)

# 添加经纬度网格线
ax.gridlines(draw_labels=True, dms=True, x_inline=False, y_inline=False)

# 添加海岸线和边界
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()

/opt/conda/lib/python3.9/site-packages/matplotlib/cm.py:494: RuntimeWarning: invalid value encountered in cast
  xx = (xx * 255).astype(np.uint8)
/opt/conda/lib/python3.9/site-packages/matplotlib/cm.py:494: RuntimeWarning: invalid value encountered in cast
  xx = (xx * 255).astype(np.uint8)

显示局部地区并叠加地图

#获取投影
dat = ds2.metpy.parse_cf('channel_0001_scaled_radiance')

geos = dat.metpy.cartopy_crs
#转换RGB
vmin = 0
vmax = .5
B1 = b1[1000:4000,1000:3500].clip(vmin, vmax) / vmax * 255
B2 = b2[1000:4000,1000:3500].clip(vmin, vmax) / vmax * 255
B3 = b3[1000:4000,1000:3500].clip(vmin, vmax) / vmax * 255

X = x[1000:4000]
Y = y[1000:3500]
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[0], X[-1], Y[-1], Y[0]))

ax2.imshow(rgb, origin='upper',extent=(X[0], X[-1], Y[-1], Y[0]), transform=geos)

# Add Coastlines and States
ax2.coastlines(resolution='50m', color='white', linewidth=1)
ax2.add_feature(ccrs.cartopy.feature.STATES, linewidth=1)
# # 添加经纬度网格线
# ax1.gridlines(draw_labels=True, dms=True, x_inline=False, y_inline=False)
# ax2.gridlines(draw_labels=True, dms=True, x_inline=False, y_inline=False)
# Add labels
plt.title('Himawari 8 \n', fontsize=20)

# Adjust tick mark size
plt.tick_params(labelsize=16)

显示局部地区并叠加地图,经纬度

#获取投影
dat = ds2.metpy.parse_cf('channel_0001_scaled_radiance')

geos = dat.metpy.cartopy_crs
#转换RGB
vmin = 0
vmax = .5
B1 = b1[1000:4000,1000:3500].clip(vmin, vmax) / vmax * 255
B2 = b2[1000:4000,1000:3500].clip(vmin, vmax) / vmax * 255
B3 = b3[1000:4000,1000:3500].clip(vmin, vmax) / vmax * 255

X = x[1000:4000]
Y = y[1000:3500]
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[0], X[-1], Y[-1], Y[0]))

ax2.imshow(rgb, origin='upper',extent=(X[0], X[-1], Y[-1], Y[0]), transform=geos)

# Add Coastlines and States
ax2.coastlines(resolution='50m', color='white', linewidth=1)
ax2.add_feature(ccrs.cartopy.feature.STATES, linewidth=1)
gl = ax1.gridlines(draw_labels=True, alpha=0.5, linewidth=1, color='w', linestyle='--')

gl.xlabels_top = False  # 关闭顶端标签
gl.ylabels_right = False  # 关闭右侧标签
gl.xformatter = LONGITUDE_FORMATTER  # x轴设为经度格式
gl.yformatter = LATITUDE_FORMATTER  # y轴设为纬度格式
plt.title('Himawari 8 \n', fontsize=20)

# Adjust tick mark size
plt.tick_params(labelsize=16)