python解析风云4B,生成红外云图、可见光云图、水汽云图

概要

提示：python实现数据读取跟png图片生成

整体流程

提示：了解hdf文件格式,使用pyhton库解析文件生成想要的卫星云图
风云4B测试数据文件请到官网下载,如果想知道hdf文件格式可下载HDFView 直接查看,但是当你打开这个文章你就不用管,直接代码就能用,

红外云图为通道Channel14 所以需要颜色,可根据实际需求自己设置cmap颜色,目前不做演示,可见光云图为通道Channel02,水汽云图为通道通道Channel11,真彩色云图需要根据通道Channel01,通道Channel02,通道Channel03进行通道融合处理,本文章不做演示,跳转这里
python解析风云4B生成真彩云图
https://blog.csdn.net/qq_38197010/article/details/147305867
大致思路:三个通道对于RGB三个颜色管道，然后合并成一个三通道图像

直接干货拿着就用!

from netCDF4 import Dataset
import matplotlib.pyplot as plt
import math
from numpy import deg2rad, rad2deg, arctan, arcsin, tan, sqrt, cos, sin
import numpy as np
from mpl_toolkits.basemap import Basemap
import h5py
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import cartopy.io.shapereader as shpreader
import matplotlib.colors as mcolors

data_file = '\\Z_SATE_C_BAWX_20250306061714_P_FY4B-_AGRI--_N_DISK_1050E_L1-_FDI-_MULT_NOM_20250306060000_20250306061459_4000M_V0001.HDF'
# 读取省级边界的 GeoJSON 或 Shapefile
china_province_shapefile = 'china.json'  # 替换为实际的 Shapefile 路径
reader = shpreader.Reader(china_province_shapefile)
ea = 6378.137  # 地球的半长轴[km]
eb = 6356.7523  # 地球的短半轴[km]
h = 42164  # 地心到卫星质心的距离[km]
λD = deg2rad(104.7)  # 卫星星下点所在经度
# 列偏移
COFF = {"0500M": 10991.5,
        "1000M": 5495.5,
        "2000M": 2747.5,
        "4000M": 1373.5}
# 列比例因子
CFAC = {"0500M": 81865099,
        "1000M": 40932549,
        "2000M": 20466274,
        "4000M": 10233137}
LOFF = COFF  # 行偏移
LFAC = CFAC  # 行比例因子


def latlon2linecolumn(lat, lon, resolution):
    """
    经纬度转行列
    (lat, lon) → (line, column)
    resolution：文件名中的分辨率{'0500M', '1000M', '2000M', '4000M'}
    line, column
    """
    # Step1.检查地理经纬度
    # Step2.将地理经纬度的角度表示转化为弧度表示
    lat = deg2rad(lat)
    lon = deg2rad(lon)
    # Step3.将地理经纬度转化成地心经纬度
    eb2_ea2 = eb ** 2 / ea ** 2
    λe = lon
    φe = arctan(eb2_ea2 * tan(lat))
    # Step4.求Re
    cosφe = cos(φe)
    re = eb / sqrt(1 - (1 - eb2_ea2) * cosφe ** 2)
    # Step5.求r1,r2,r3
    λe_λD = λe - λD
    r1 = h - re * cosφe * cos(λe_λD)
    r2 = -re * cosφe * sin(λe_λD)
    r3 = re * sin(φe)
    # Step6.求rn,x,y
    rn = sqrt(r1 ** 2 + r2 ** 2 + r3 ** 2)
    x = rad2deg(arctan(-r2 / r1))
    y = rad2deg(arcsin(-r3 / rn))
    # Step7.求c,l
    column = COFF[resolution] + x * 2 ** -16 * CFAC[resolution]
    line = LOFF[resolution] + y * 2 ** -16 * LFAC[resolution]
    return np.rint(line).astype(np.uint16), np.rint(column).astype(np.uint16)


# 中国范围
x_min = 10
x_max = 60
y_min = 70
y_max = 140

column = math.ceil((x_max - x_min) / 0.04)
row = math.ceil((y_max - y_min) / 0.04)
print(row, column)
ynew = np.linspace(y_min, y_max, row)  # 获取网格y
xnew = np.linspace(x_min, x_max, column)  # 获取网格x
xnew, ynew = np.meshgrid(xnew, ynew)  # 生成xy二维数组
data_grid = np.zeros((row, column))  # 声明一个二维数组

keyword = "NOMChannel"


cmap = mcolors.LinearSegmentedColormap.from_list("smooth_gradient", colors, N=256)
# cmap = mcolors.LinearSegmentedColormap.from_list("custom_IR", colors, N=256)
# cmap = plt.cm.viridis

# 1. 读取观测数据文件
with h5py.File(data_file, 'r') as f_data:
    type = f_data['Data'].keys()
    nc_obj=f_data['Data'];
    index = {}
    r_data = {}
    for k in type:
        if str(k)!= 'NOMChannel14':
            continue
        if str(k).find(keyword) == 0:

            value = nc_obj[k][:]
            for i in range(row):
                for j in range(column):
                    lat = xnew[i][j]
                    lon = ynew[i][j]
                    fy_line = 0
                    fy_column = 0
                    if index.get((lat, lon)) == None:
                        # 查找行列并记录下来下次循环使用
                        fy_line, fy_column = latlon2linecolumn(lat, lon, "4000M")
                        index[(lat, lon)] = fy_line, fy_column
                    else:
                        fy_line, fy_column = index.get((lat, lon))
                    data_grid[i][j] = value[fy_line, fy_column]
            r_data[k] = data_grid
            img = plt.figure()
            ax = img.add_subplot(111)
            m = Basemap(llcrnrlon=y_min, llcrnrlat=x_min, urcrnrlon=y_max, urcrnrlat=x_max)
            # 4. 创建地图画布
            plt.figure(figsize=(12, 8))
            ax = plt.axes(projection=ccrs.PlateCarree())  # 确保使用平面投影
            # 关键修正：移除转置操作
            if 'Channel14' in k:
                img = ax.pcolormesh(ynew, xnew, data_grid,  # 直接使用bt_data
                                    cmap='viridis', # 更改为彩色映射
                                    vmin=data_grid.min(),  # 动态范围
                                    vmax=data_grid.max(),
                                    transform=ccrs.PlateCarree())
            else:
                img = ax.pcolormesh(ynew, xnew, data_grid,  # 直接使用bt_data
                                    cmap='gray',  # 更改为彩色映射
                                    # vmin=data_grid.min(),  # 动态范围
                                    # vmax=data_grid.max(),
                                    transform=ccrs.PlateCarree())
            # 4. 使用Cartopy内置方法添加边界
            ax.add_geometries(
                reader.geometries(),  # 自动处理所有几何类型
                crs=ccrs.PlateCarree(),  # 数据坐标系（WGS84）
                edgecolor='white',
                facecolor='none',  # 不填充颜色
                linewidth=0.6
            )
            # 添加地理要素
            ax.add_feature(cfeature.COASTLINE.with_scale('50m'), linestyle=':',linewidth=0.5,color='white')
            ax.add_feature(cfeature.BORDERS.with_scale('50m'),  linewidth=0.5,color='white')

            # 保存图像
            plt.savefig(f'./images/FY4B_China_IR_Color {k}.png', dpi=150, bbox_inches='tight')
            print("通道" + k + "绘制完成")
    plt.close()

小结

提示：这里可以添加总结

很多东西一般我都喜欢用Java处理,处理很麻烦最后才会选择python,针对部分气象数据,python优势显著