前言
参考文章:Python 解析风云四A卫星L1级别数据以及绘制卫星云图
一、数据下载
- 下载地址: https://satellite.nsmc.org.cn/DataPortal/cn/home/index.html
- 数据格式说明 https://img.nsmc.org.cn/PORTAL/NSMC/DATASERVICE/DataFormat/FY4B/Data/Format/FY4B_AGRI_DISK_L1_4KM_HDF__V1.0_20220610.pdf
由于本人使用的是L1的4km分辨率的数据,因此放的是此数据的格式说明,其他数据的格式官网里也有,可以自行查找并微调代码。
二、解析HDF数据
1. 导入数据
本人使用的是netCDF4导入,由于气象数据大部分是nc格式,所以就统一使用这个库包来处理数据。
from netCDF4 import Dataset
hdf_data_path = "FY4B-_AGRI--_N_DISK_1050E_L1-_FDI-_MULT_NOM_20250306000000_20250306001459_4000M_V0001.HDF"
# netCDF4 解析
nc_obj = Dataset(hdf_data_path)
2. 数据格式与变量名查看
# 直接查看nc_obj变量(下面不是代码,是变量格式)
<class 'netCDF4.Dataset'>
root group (NETCDF4 data model, file format HDF5):
Responser: NSMC
Version Of Software: V1000
Software Revision Date: 2021-06-18
Observing Beginning Date: 2025-03-06
Observing Beginning Time: 00:00:02.345
Observing Ending Date: 2025-03-06
Observing Ending Time: 00:12:54.300
Data Creating Date: 2025-03-06
Data Creating Time: 00:14:23.651
AdditionalAnnotation: shinetek
OBIType: DISK
ProductID: FY4B-_AGRI--_N_DISK_1050E_L1-_FDI-_MULT_NOM_20250306000000_20250306001459_4000M_V0001.HDF
ProductName: FY4B-_AGRI--_N_DISK_1050E_L1-_FDI-_MULT_NOM_20250306000000_20250306001459_4000M_V0001.HDF
Satellite Name: FY-4B
Sensor Name: AGRI
Sensor Identification Code: AGRI
Dataset Name: MULT
File Name: FY4B-_AGRI--_N_DISK_1050E_L1-_FDI-_MULT_NOM_20250306000000_20250306001459_4000M_V0001.HDF
File Alias Name: FY4B-_AGRI--_N_DISK_1050E_L1-_FDI-_MULT_NOM_20250306000000_20250306001459_4000M_V0001.HDF
Version Of Coefficient Index: V1003
Coefficient Index Revision Date: 2023-05-22
Data Quality: 0
Number Of Scans: 695
Incomplete Scans: 65535
QA_Scan_Flag: 0
QA_Pixel_Flag: 0
Begin Line Number: 0
End Line Number: 2747
Begin Pixel Number: 0
End Pixel Number: 2747
NOMCenterLat: 0.0
NOMCenterLon: 105.0
NOMSatHeight: 35786000.0
RegCenterLon: 65535.0
RegCenterLat: 65535.0
RegLength: 2748.0
RegWidth: 2748.0
dEA: 6377830.0
dSamplingAngle: 112.0
dSteppingAngle: 112.0
dObRecFlat: 298.257222101
Earth/Sun Distance Ratio: 0.9920520564
Semimajor axis of ellipsoid: 6378137.0
Semiminor axis of ellipsoid: 6356752.0
Flattening of ellipsoid: 0.003352810681182319
Orbit Point Latitude: [65535. 65535. 65535. 65535.]
Orbit Point Longitude: [65535. 65535. 65535. 65535.]
Circuit A/B Flag: 0
On Board Process Flag: 0
Sun Avoidance Flag: 0
dimensions(sizes):
variables(dimensions):
groups: Calibration, Data, NOMObs, QA, VerSoft
3. 绘制指定通道云图
import matplotlib.pyplot as plt
data = nc_obj['Data/NOMChannel13'][:] # 13通道
plt.imshow(data, cmap='gray')
plt.title("20250306 NOMChannel13")
plt.axis('off')
plt.show()
13通道全球云图
三、数据坐标转换
根据数据格式文档,数据集Data里就是DN值(L1 数据产品,NOMChannel07-15 存放的是定标后辐亮度经过线性变换的 DN 值),但是经纬度需要根据卫星的标称行列号进行转换。FY-4 卫星采用 CGMSLRIT/HRIT 全球规范定义的静止轨道标称投影,地理坐标基于 WGS84 参考椭球计算得到。现在需要根据文档里的公式进行转换,文档如图:
# 卫星参数配置
ea = 6378.137 # 地球半长轴 [km]
eb = 6356.7523 # 地球短半轴 [km]
h = 42164 # 卫星高度 [km]
λD = deg2rad(105.0) # 卫星星下点经度
# 列偏移
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 geolocation_to_linecolumn(lat, lon):
"""地理坐标转行列号(向量化计算)"""
eb2_ea2 = (eb**2)/(ea**2)
φ = np.radians(lat)
λ = np.radians(lon)
φe = arctan(eb2_ea2 * np.tan(φ))
cosφe = np.cos(φe)
re = eb / np.sqrt(1 - (1 - eb2_ea2)*cosφe**2)
λ_diff = λ - λD
r1 = h - re*cosφe*np.cos(λ_diff)
r2 = -re*cosφe*np.sin(λ_diff)
r3 = re*np.sin(φe)
rn = np.sqrt(r1**2 + r2**2 + r3**2)
x = np.degrees(arctan(-r2/r1))
y = np.degrees(arcsin(-r3/rn))
# 4000M分辨率参数
(本人直接使用了4000m分辨率的列偏移和列比例因子,其他分辨率的请按照上面的数组直接修改数据即可)
COFF, CFAC = 1373.5, 10233137
column = COFF + x * CFAC / 65536
line = COFF + y * CFAC / 65536
return np.clip(line, 0, None).astype(int), np.clip(column, 0, None).astype(int)
四、绘制中国DN值图
L1 数据产品,NOMChannel07-15 存放的是定标后辐亮度经过线性变换的 DN 值,直接使用就行。
# 读取数据
nc = Dataset("FY4B-_AGRI--_N_DISK_1050E_L1-_FDI-_MULT_NOM_20250306000000_20250306001459_4000M_V0001.HDF")
dn_data = nc['Data/NOMChannel13'][:] # DN值数据
# 生成地理网格
lon, lat = np.meshgrid(np.linspace(lon_min, lon_max, 1500),
np.linspace(lat_min, lat_max, 1000))
# 转换行列号
lines, cols = geolocation_to_linecolumn(lat, lon)
# 提取有效数据
valid = (lines < dn_data.shape[0]) & (cols < dn_data.shape[1])
result = np.full(lat.shape, np.nan)
result[valid] = dn_data[lines[valid], cols[valid]]
# 指定绘制中国区域范围数据(为了方便出南海小图数据,可自行更改)
lat_min, lat_max = 0.0, 56.0
lon_min, lon_max = 72.0, 136.0
绘图函数完整代码,需调整相关要素可自行更改。
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.io.shapereader as shpreader
import matplotlib.patches as mpatches
from matplotlib.ticker import AutoMinorLocator, FixedLocator
plt.rcParams.update({
'font.family': 'serif', # 主字体类型
'font.serif': 'Times New Roman', # 指定新罗马字体
'mathtext.fontset': 'stix', # 数学公式字体对齐
'axes.unicode_minus': False # 解决负号显示问题
})
def plot_geo_data(data, lon, lat,
main_extent=None,
inset_position=[0.75, 0.1, 0.15, 0.2],
title='Satellite Data Visualization',
shp_path='中国省市县/中国国界与省界-含九段线/bou2_4l.shp',
cmap='jet',
vmin=None,
vmax=None,
cbar_label='Value',
figsize=(15, 10),
dpi=300,
gridline_style={'linewidth':0.8, 'color':'dimgray', 'alpha':0.6}):
"""
绘制地理空间数据可视化图
参数:
data : 2D numpy数组,要可视化的数据
lon : 经度坐标数组
lat : 纬度坐标数组
main_extent : list, optional
主图范围 [lon_min, lon_max, lat_min, lat_max]
默认根据数据范围自动计算
inset_extent : list, optional
缩略图范围 [lon_min, lon_max, lat_min, lat_max]
默认南海区域 [105, 125, 2, 25]
inset_position : list, optional
缩略图位置 [left, bottom, width, height]
title : 图标题 (默认: 'Satellite Data Visualization')
shp_path : shapefile路径
cmap : 颜色映射 (默认: 'jet')
vmin : 颜色条最小值 (默认: 数据最小值)
vmax : 颜色条最大值 (默认: 数据最大值)
cbar_label : 颜色条标签 (默认: 'Value')
figsize : 图像尺寸 (默认: (15,10))
dpi : 分辨率 (默认: 100)
gridline_style : 网格线样式字典 (默认: 灰色虚线)
"""
# 初始化画布和投影
proj = ccrs.PlateCarree()
fig = plt.figure(figsize=figsize, dpi=dpi)
# 主图坐标系
main_ax = fig.add_subplot(111, projection=proj)
# 设置主图范围
if main_extent is not None:
main_ax.set_extent(main_extent, crs=proj)
else: # 自动计算数据范围
buffer = 1 # 边界缓冲度数
main_extent = [lon.min()-buffer, lon.max()+buffer,
lat.min()-buffer, lat.max()+buffer]
main_ax.set_extent(main_extent, crs=proj)
# 设置颜色范围
vmin = np.nanmin(data) if vmin is None else vmin
vmax = np.nanmax(data) if vmax is None else vmax
# 绘制主图数据
img = main_ax.pcolormesh(lon, lat, data,
cmap=cmap,
vmin=vmin,
vmax=vmax,
transform=proj)
# 添加地理要素
try:
main_ax.add_geometries(
shpreader.Reader(shp_path).geometries(),
crs=proj,
edgecolor='black',
facecolor='none',
linewidth=1
)
except Exception as e:
print(f"Error loading shapefile: {e}")
# 配置经纬网
gl = main_ax.gridlines(
crs=proj,
draw_labels=True,
linestyle=(0, (5, 10)),
**gridline_style
)
gl.top_labels = False
gl.right_labels = False
gl.xlabel_style = {'size': 10, 'color': 'black'}
gl.ylabel_style = {'size': 10, 'color': 'black'}
# 添加颜色条
cbar = plt.colorbar(img, ax=main_ax, extend='both', shrink=0.8)
cbar.set_label(cbar_label, fontsize=12, labelpad=12)
cbar.ax.tick_params(labelsize=10)
# 设置标题
plt.title(title, pad=20,
fontdict={'size': 20,
'weight': 'bold'},
loc='center'
)
# ------------------ 南海子图部分 ------------------
# 设置南海显示范围
scs_extent = [105, 125, 2, 25]
# 创建子图(右下角位置)
inset_ax = fig.add_axes(inset_position, projection=proj)
inset_ax.set_extent(scs_extent, crs=proj)
# 绘制相同数据
inset_ax.pcolormesh(lon, lat, data,
cmap=cmap,
vmin=vmin,
vmax=vmax,
transform=proj)
# 添加地理边界
try:
inset_ax.add_geometries(
shpreader.Reader(shp_path).geometries(),
crs=proj,
edgecolor='black',
facecolor='none',
linewidth=0.8 # 更细的边界线
)
except Exception as e:
print(f"Error loading shapefile in inset: {e}")
return fig, main_ax
# dn分布图
fig2, ax2 = plot_geo_data(
data=result,
lon=lon,
lat=lat,
main_extent=[73, 135, 15, 54],
inset_position=[0.61, 0.205, 0.15, 0.2], # [left, bottom, width, height]
title='20250306 FY4B AGRI Channel13 DN Values - China Region',
cmap=plt.cm.get_cmap('jet', 15),
cbar_label='DN Values',
figsize=(12, 8),
gridline_style={'linewidth':1.0, 'color':'dimgray', 'alpha':0.8}
)
plt.show()
五、绘制中国亮温图
根据数据文档,CALChannelXX 为 SR 与亮温(K)的查找表,用于将 DN 值转换成亮温,直接使用索引将DN值转换为亮温。
# 配置参数
CHANNEL = 13
INVALID_EARTH = 65534
INVALID_SPACE = 65535
# 解析校准数据
bt_lut = nc['Calibration/CALChannel13'][:] # 直接获取一维亮温数组
# 初始化亮温数组
bt_result = np.full_like(lat, np.nan, dtype=np.float32)
# 数据提取
valid = (lines < dn_data.shape[0]) & (cols < dn_data.shape[1])
raw_dn = dn_data[lines[valid], cols[valid]]
# 创建有效性掩膜
valid_dn = ~((raw_dn == INVALID_EARTH) | (raw_dn == INVALID_SPACE))
in_range = (raw_dn >= 0) & (raw_dn < len(bt_lut)) # 确保DN在LUT范围内
# 执行转换(直接查表)
final_valid = valid_dn & in_range
bt_values = bt_lut[raw_dn[final_valid].astype(int)]
# 获取所有有效像素的二维坐标
valid_coords = np.argwhere(valid) # shape=(N,2)
# 应用最终有效性筛选
selected_coords = valid_coords[final_valid]
# 分离行列索引
rows = selected_coords[:, 0]
cols = selected_coords[:, 1]
# 直接赋值到二维数组
bt_result[rows, cols] = bt_values
# 亮温分布图
fig1, ax1 = plot_geo_data(
data=bt_result,
lon=lon,
lat=lat,
main_extent=[73, 135, 15, 54],
inset_position=[0.61, 0.205, 0.15, 0.2], # [left, bottom, width, height]
title='20250306 FY4B AGRI Channel13 Brightness Temperature - China Region',
cmap=plt.cm.get_cmap('jet', 15),
cbar_label='Brightness Temperature (K)',
figsize=(12, 8),
gridline_style={'linewidth':1.0, 'color':'dimgray', 'alpha':0.8}
)
plt.show()
完整代码
from netCDF4 import Dataset
import h5py
import matplotlib.pyplot as plt
import numpy as np
import cartopy.crs as ccrs
import cartopy.io.shapereader as shpreader
import cartopy.feature as cfeature
from math import radians as deg2rad
from numpy import arctan, sin, cos, sqrt, arcsin
hdf_data_path = "FY4B-_AGRI--_N_DISK_1050E_L1-_FDI-_MULT_NOM_20250306000000_20250306001459_4000M_V0001.HDF"
# netCDF4 解析
nc_obj = Dataset(hdf_data_path)
data = nc_obj['Data/NOMChannel13'][:]
plt.imshow(data, cmap='gray')
plt.title("20250306 NOMChannel13")
plt.axis('off')
plt.show()
# 卫星参数配置
ea = 6378.137 # 地球半长轴 [km]
eb = 6356.7523 # 地球短半轴 [km]
h = 42164 # 卫星高度 [km]
λD = deg2rad(105.0) # 卫星星下点经度
# 中国区域范围
lat_min, lat_max = 0.0, 56.0
lon_min, lon_max = 72.0, 136.0
# 坐标转换函数
def geolocation_to_linecolumn(lat, lon):
"""地理坐标转行列号(向量化计算)"""
eb2_ea2 = (eb**2)/(ea**2)
φ = np.radians(lat)
λ = np.radians(lon)
φe = arctan(eb2_ea2 * np.tan(φ))
cosφe = np.cos(φe)
re = eb / np.sqrt(1 - (1 - eb2_ea2)*cosφe**2)
λ_diff = λ - λD
r1 = h - re*cosφe*np.cos(λ_diff)
r2 = -re*cosφe*np.sin(λ_diff)
r3 = re*np.sin(φe)
rn = np.sqrt(r1**2 + r2**2 + r3**2)
x = np.degrees(arctan(-r2/r1))
y = np.degrees(arcsin(-r3/rn))
# 4000M分辨率参数
COFF, CFAC = 1373.5, 10233137
column = COFF + x * CFAC / 65536
line = COFF + y * CFAC / 65536
return np.clip(line, 0, None).astype(int), np.clip(column, 0, None).astype(int)
# 读取数据
nc = Dataset("FY4B-_AGRI--_N_DISK_1050E_L1-_FDI-_MULT_NOM_20250306000000_20250306001459_4000M_V0001.HDF")
dn_data = nc['Data/NOMChannel13'][:] # DN值数据
# 生成地理网格
lon, lat = np.meshgrid(np.linspace(lon_min, lon_max, 1500),
np.linspace(lat_min, lat_max, 1000))
# 转换行列号
lines, cols = geolocation_to_linecolumn(lat, lon)
# 提取有效数据
valid = (lines < dn_data.shape[0]) & (cols < dn_data.shape[1])
result = np.full(lat.shape, np.nan)
result[valid] = dn_data[lines[valid], cols[valid]]
# 配置参数
CHANNEL = 13
INVALID_EARTH = 65534
INVALID_SPACE = 65535
# 解析校准数据
bt_lut = nc['Calibration/CALChannel13'][:] # 直接获取一维亮温数组
# 初始化亮温数组
bt_result = np.full_like(lat, np.nan, dtype=np.float32)
# 数据提取
valid = (lines < dn_data.shape[0]) & (cols < dn_data.shape[1])
raw_dn = dn_data[lines[valid], cols[valid]]
# 创建有效性掩膜
valid_dn = ~((raw_dn == INVALID_EARTH) | (raw_dn == INVALID_SPACE))
in_range = (raw_dn >= 0) & (raw_dn < len(bt_lut)) # 确保DN在LUT范围内
# 执行转换(直接查表)
final_valid = valid_dn & in_range
bt_values = bt_lut[raw_dn[final_valid].astype(int)]
# 获取所有有效像素的二维坐标
valid_coords = np.argwhere(valid) # shape=(N,2)
# 应用最终有效性筛选
selected_coords = valid_coords[final_valid]
# 分离行列索引
rows = selected_coords[:, 0]
cols = selected_coords[:, 1]
# 直接赋值到二维数组
bt_result[rows, cols] = bt_values
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.io.shapereader as shpreader
import matplotlib.patches as mpatches
from matplotlib.ticker import AutoMinorLocator, FixedLocator
plt.rcParams.update({
'font.family': 'serif', # 主字体类型
'font.serif': 'Times New Roman', # 指定新罗马字体
'mathtext.fontset': 'stix', # 数学公式字体对齐
'axes.unicode_minus': False # 解决负号显示问题
})
def plot_geo_data(data, lon, lat,
main_extent=None,
inset_position=[0.75, 0.1, 0.15, 0.2],
title='Satellite Data Visualization',
shp_path='中国省市县/中国国界与省界-含九段线/bou2_4l.shp',
cmap='jet',
vmin=None,
vmax=None,
cbar_label='Value',
figsize=(15, 10),
dpi=300,
gridline_style={'linewidth':0.8, 'color':'dimgray', 'alpha':0.6}):
"""
绘制地理空间数据可视化图
参数:
data : 2D numpy数组,要可视化的数据
lon : 经度坐标数组
lat : 纬度坐标数组
main_extent : list, optional
主图范围 [lon_min, lon_max, lat_min, lat_max]
默认根据数据范围自动计算
inset_extent : list, optional
缩略图范围 [lon_min, lon_max, lat_min, lat_max]
默认南海区域 [105, 125, 2, 25]
inset_position : list, optional
缩略图位置 [left, bottom, width, height]
title : 图标题 (默认: 'Satellite Data Visualization')
shp_path : shapefile路径
cmap : 颜色映射 (默认: 'jet')
vmin : 颜色条最小值 (默认: 数据最小值)
vmax : 颜色条最大值 (默认: 数据最大值)
cbar_label : 颜色条标签 (默认: 'Value')
figsize : 图像尺寸 (默认: (15,10))
dpi : 分辨率 (默认: 100)
gridline_style : 网格线样式字典 (默认: 灰色虚线)
"""
# 初始化画布和投影
proj = ccrs.PlateCarree()
fig = plt.figure(figsize=figsize, dpi=dpi)
# 主图坐标系
main_ax = fig.add_subplot(111, projection=proj)
# 设置主图范围
if main_extent is not None:
main_ax.set_extent(main_extent, crs=proj)
else: # 自动计算数据范围
buffer = 1 # 边界缓冲度数
main_extent = [lon.min()-buffer, lon.max()+buffer,
lat.min()-buffer, lat.max()+buffer]
main_ax.set_extent(main_extent, crs=proj)
# 设置颜色范围
vmin = np.nanmin(data) if vmin is None else vmin
vmax = np.nanmax(data) if vmax is None else vmax
# 绘制主图数据
img = main_ax.pcolormesh(lon, lat, data,
cmap=cmap,
vmin=vmin,
vmax=vmax,
transform=proj)
# 添加地理要素
try:
main_ax.add_geometries(
shpreader.Reader(shp_path).geometries(),
crs=proj,
edgecolor='black',
facecolor='none',
linewidth=1
)
except Exception as e:
print(f"Error loading shapefile: {e}")
# 配置经纬网
gl = main_ax.gridlines(
crs=proj,
draw_labels=True,
linestyle=(0, (5, 10)),
**gridline_style
)
gl.top_labels = False
gl.right_labels = False
gl.xlabel_style = {'size': 10, 'color': 'black'}
gl.ylabel_style = {'size': 10, 'color': 'black'}
# 添加颜色条
cbar = plt.colorbar(img, ax=main_ax, extend='both', shrink=0.8)
cbar.set_label(cbar_label, fontsize=12, labelpad=12)
cbar.ax.tick_params(labelsize=10)
# 设置标题
plt.title(title, pad=20,
fontdict={'size': 20,
'weight': 'bold'},
loc='center'
)
# ------------------ 南海子图部分 ------------------
# 设置南海显示范围
scs_extent = [105, 125, 2, 25]
# 创建子图(右下角位置)
inset_ax = fig.add_axes(inset_position, projection=proj)
inset_ax.set_extent(scs_extent, crs=proj)
# 绘制相同数据
inset_ax.pcolormesh(lon, lat, data,
cmap=cmap,
vmin=vmin,
vmax=vmax,
transform=proj)
# 添加地理边界
try:
inset_ax.add_geometries(
shpreader.Reader(shp_path).geometries(),
crs=proj,
edgecolor='black',
facecolor='none',
linewidth=0.8 # 更细的边界线
)
except Exception as e:
print(f"Error loading shapefile in inset: {e}")
return fig, main_ax
# dn分布图
fig2, ax2 = plot_geo_data(
data=result,
lon=lon,
lat=lat,
main_extent=[73, 135, 15, 54],
inset_position=[0.61, 0.205, 0.15, 0.2], # [left, bottom, width, height]
title='20250306 FY4B AGRI Channel13 DN Values - China Region',
cmap=plt.cm.get_cmap('jet', 15),
cbar_label='DN Values',
figsize=(12, 8),
gridline_style={'linewidth':1.0, 'color':'dimgray', 'alpha':0.8}
)
plt.show()
# 亮温分布图
fig1, ax1 = plot_geo_data(
data=bt_result,
lon=lon,
lat=lat,
main_extent=[73, 135, 15, 54],
inset_position=[0.61, 0.205, 0.15, 0.2], # [left, bottom, width, height]
title='20250306 FY4B AGRI Channel13 Brightness Temperature - China Region',
cmap=plt.cm.get_cmap('jet', 15),
cbar_label='Brightness Temperature (K)',
figsize=(12, 8),
gridline_style={'linewidth':1.0, 'color':'dimgray', 'alpha':0.8}
)
plt.show()