深度文章 | 我国FY-4A气象卫星产品的处理与分析
本帖最后由 好奇心Log 于 2024-5-1 22:17 编辑风云系列气象卫星
我国早在二十世纪七十年代就开始发展我国的气象卫星，分别实现了极轨卫星和静止卫星的业务化运行，是继美国、俄罗斯之后第三个同时拥有极轨气象卫星和静止气象卫星的国家。在世界气象组织（WMO）协调和管理下，世界各国已经实现了气象观测资料的交换和获取全球气象观测资料的目的，全球业务气象卫星探测系统就是在WMO协调下，通过各气象卫星发射和运行的国家共同努力建成的天基探测系统。中国的业务气象卫星风云三号和风云二号已经成为了全球业务气象卫星探测系统的重要成员。他们在覆盖范围和分辨率上相辅相成，成为中国获取全球资料和满足区域灾害性天气和环境监测、气象服务和地球系统科学研究的重要工具。
风云一号和三号为极轨气象卫星，风云二号和四号为静止气象卫星。
风云系列卫星的发射时间轴
风云极轨气象卫星基本信息
风云静止气象卫星基本信息
风云卫星系列仪器
以上信息来自国家卫星气象中心：
【1】http://gsics.nsmc.org.cn/nsmc/cn/satellite/index.html
【2】http://gsics.nsmc.org.cn/nsmc/cn/instrument/index.html
风云四号气象卫星
风云四号气象卫星是我国第二代静止气象卫星，充分考虑海洋、农业、林业、水利以及环境、空间科学等领域的需求，以实现综合利用。风云四号的主要发展目标是：卫星X射线太阳观测，加强空间天气监测预警。风云四号卫星计划发展光学和微波两种类型的卫星。
2016年12月11日，风云四号科研试验卫星在西昌卫星发射中心成功发射。
2016年12月17日定点于东经99.5度赤道上空静止轨道位置，并正式命名为风云四号A星（又称FY-4A，简称“A星”）。
2017年2月27日，随着我国新一代静止气象卫星风云四号A星获取首批图像和数据，世界第一幅静止轨道地球大气高光谱图正式亮相，与此同时，我国首次获取彩色卫星云图和闪电分布图。
2017年5月25日起，FY-4A卫星漂移到风云静止卫星的主业务位置（东经105度）附近的东经104.7度赤道上空，提早做好业务值班准备。
2017年9月25日，风云四号A星正式交付用户投入使用，目前定点在105E，FY-4A的成功发射在技术体制上实现了更新换代，该星考核寿命为5年。
风云四号B星为风云四号系列气象卫星首颗业务卫星，于2021年6月3日0时17分，在西昌卫星发射中心搭乘长征三号乙运载火箭发射。6月10日17时7分，风云四号B星（以下简称“B星”）成功定点东经123.5°赤道上空静止轨道位置。B星是我国新一代静止轨道气象卫星风云四号系列卫星的首发业务星，B星考核寿命由A星的5年提升至7年。它的成功发射，标志着我国新一代静止轨道卫星观测系统正式进入业务化发展阶段，对确保我国静止气象卫星升级换代和连续、可靠、稳定业务运行意义重大。
B星发射后与A星一起双星组网进一步满足我国及“一带一路”沿线国家和地区气象监测预报、应急防灾减灾等服务需求。
FY-4A卫星装载了多通道扫描成像辐射计、干涉式大气垂直探测仪、闪电成像仪和空间环境监测仪器。
FY-4A卫星扫描成像辐射计主要承担获取云图的任务，共14通道，是风云二号5通道的近三倍，在风云二号观测云、水汽、植被、地表的基础上，还具备了捕捉气溶胶、雪的能力，并且能清晰区分云的不同相态和高、中层水汽。相比于风云二号单一可见光通道的限制，风云四号首次制作出彩色卫星云图，最快1分钟生成一次区域观测图像。
干涉式大气垂直探测仪是国际上第一台在静止轨道上以红外高光谱干涉分光方式探测大气垂直结构的精密遥感仪器，在静止轨道上从二维观测进入三维综合观测。
闪电成像仪为亚太地区首次研制发射的同类载荷，测试数据印证了其可对我国及周边区域闪电进行探测，进而实现强对流天气的监测和跟踪，提供闪电灾害预警。此外，空间环境监测仪器具备监测太阳活动和空间环境的能力，探测通道数量、探测精度显著提高。
http://www.nsmc.org.cn/NSMC/Channels/FY4A_en.html
FY4A官方应用平台
网址：http://rsapp.nsmc.org.cn/geofy/
FY-4A数据下载
注册账号：http://data.cma.cn/
实时数据(30天内)
下载地址：https://fy4.nsmc.org.cn/data/en/data/realtime.html
历史数据 (2018-03-12 -- now)
网址：http://satellite.nsmc.org.cn/PortalSite/Data/Satellite.aspx
FY4A AGRI L1数据
多通道扫描成像辐射计（AGRI） 是风云四号静止气象卫星的主要载荷之一，通过精密的双扫描镜机构实现精确和灵活的二维指向，可实现分钟级的区域快速扫描；采用离轴三反主光学系统，高频次获取14波段的地球云图，并利用星上黑体进行高频次红外定标，以确保观测数据的精度。
1.全圆盘 FY4A-_AGRI--_N_DISK _1047E_L1-_FDI-_MULT_NOM_20210726030000_20210726031459_4000M_V0001.HDF
2.中国区 FY4A-_AGRI--_N_REGC _1047E_L1-_FDI-_MULT_NOM_20210726033000_20210726033417_1000M_V0001.HDF
全圆盘标识：DISK, 中国区标识：REGC.
咸迪,方翔,贾煦,等. 风云四号气象卫星天气应用平台及其应用. 卫星应用,2020(2):20-24. DOI:10.3969/j.issn.1674-9030.2020.02.008.
导入模块
!pip install satpy==0.29.0
import os, glob
from satpy.scene import Scene
from pyresample import get_area_def
import matplotlib.pyplot as plt
import gdal
from numpy import deg2rad, rad2deg, arctan, arcsin, tan, sqrt, cos, sin
import warnings
warnings.filterwarnings('ignore')
解析FY-4A数据
Satpy支持的卫星数据（ https://satpy.readthedocs.io/en/latest/index.html#reader-table ）
# 加载FY4A文件
filenames = glob.glob('/home/mw/input/FY4A9106/FY4A-_AGRI*4000M_V0001.HDF')
# 创建scene对象
scn = Scene(filenames, reader='agri_l1')# 查看可用的通道
scn.available_dataset_names()['C01',
'C02',
'C03',
'C04',
'C05',
'C06',
'C07',
'C08',
'C09',
'C10',
'C11',
'C12',
'C13',
'C14',
'satellite_azimuth_angle',
'satellite_zenith_angle',
'solar_azimuth_angle',
'solar_glint_angle',
'solar_zenith_angle']
遥感数据可视化
绘制单波段影像
# 以红外通道为例
ir_channel = 'C12'
scn.load()
scn.save_dataset(ir_channel, filename='C12.png')# 在Notebook中显示
scn.show(ir_channel)
# 查看可用的合成选项
scn.available_composite_names()['ash',
'cloud_phase',
'cloud_phase_distinction',
'cloud_phase_distinction_raw',
'cloud_phase_raw',
'color_infrared',
'day_microphysics',
'day_microphysics_agri',
'day_microphysics_eum',
'dust',
'fire_temperature_awips',
'fog',
'green',
'green_snow',
'ir108_3d',
'ir_cloud_day',
'land_cloud',
'land_cloud_fire',
'natural_color',
'natural_color_raw',
'natural_enh',
'overview',
'overview_raw',
'snow',
'snow_fog',
'true_color']
绘制真彩色影像
composite = 'true_color'
scn.load()
scn.save_dataset(composite, filename='true_color.png')Required file type 'agri_l1_4000m_geo' not found or loaded for 'satellite_zenith_angle'
Required file type 'agri_l1_4000m_geo' not found or loaded for 'solar_azimuth_angle'
Required file type 'agri_l1_4000m_geo' not found or loaded for 'solar_zenith_angle'
Required file type 'agri_l1_4000m_geo' not found or loaded for 'satellite_azimuth_angle'
scn.show(composite)
绘制自定义区域
area_id = '01'
x_size = 332
y_size = 276
area_extent = (-551748.414108, -658226.923988, 551748.414108, 669774.963851)
projection = '+proj=laea +lat_0=31.0 +lon_0=120.0 +ellps=WGS84'
description = "myarea"
proj_id = 'yanhua_120._31.0'
areadef = get_area_def(area_id, description, proj_id, projection,x_size, y_size, area_extent)
yanhua_scene = scn.resample(areadef)
yanhua_scene.show(composite)
7月24日9时，国家卫星气象中心迅速进入台风二级应急响应状态，加强监测服务，利用极轨卫星和静止卫星协同开展对台风“烟花”的实时、连续、动态监测。7月25日12时30分前后，台风“烟花”在浙江省舟山普陀沿海登陆，登陆时中心附近最大风力13级（38米/秒），中心最低气压为965百帕。7月26日9时50分前后，台风“烟花”在浙江嘉兴平湖市沿海再次登陆，登陆时中心附近最大风力10级（强热带风暴，28米/秒）。
风云四号气象卫星先进静止轨道辐射成像仪（AGRI）可以对台风实施定位和定强，也可以反映台风的云顶高度信息。
（来源：国家卫星气象中心遥感应用服务中心）
网址：http://www.cma.gov.cn/2011xzt/2021zt/2021tf/20200902/2018070902/202107/t20210725_581629.html
遥感数据重投影
圆盘标称投影
上述绘制特定区域的脚本其实对投影做了处理，这里介绍一下原理和方法，有兴趣的可以自己尝试建立转换关系。
圆盘标称投影数据时静止气象卫星常见的数据产品，也就是Geostationary投影，主要的投影参数有中央经度和卫星高度。
王静,刘成保,杨磊,商建,张志清.静止气象卫星标称网格的计算方法及其在风云四号中的应用.光学学报,2018,38(12):135-143.
行列号与经纬度互换
网址：https://github.com/Mo-Dabao/BiteFengyun
ea = 6378.137# 地球的半长轴
eb = 6356.7523# 地球的短半轴
h = 42164# 地心到卫星质心的距离
λ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 + x * 2**-16 * CFAC
line = LOFF + y * 2**-16 * LFAC
return line, columndef linecolumn2latlon(line, column, resolution):
"""
(line, column) → (lat, lon)
resolution：文件名中的分辨率{'0500M', '1000M', '2000M', '4000M'}
"""
# Step1.求x,y
x = deg2rad((column - COFF) / (2**-16 * CFAC))
y = deg2rad((line - LOFF) / (2**-16 * LFAC))
# Step2.求sd,sn,s1,s2,s3,sxy
cosx = cos(x)
cosy = cos(y)
siny = sin(y)
cos2y = cosy**2
hcosxcosy = h * cosx * cosy
cos2y_ea_eb_siny_2 = cos2y + (ea / eb * siny)**2
sd = sqrt(hcosxcosy**2 - cos2y_ea_eb_siny_2 * (h**2 - ea**2))
sn = (hcosxcosy - sd) / cos2y_ea_eb_siny_2
s1 = h - sn * cosx * cosy
s2 = sn * sin(x) * cosy
s3 = -sn * siny
sxy = sqrt(s1**2 + s2**2)
# Step3.求lon,lat
lon = rad2deg(arctan(s2 / s1) + λD)
lat = rad2deg(arctan(ea**2 / eb**2 * s3 / sxy))
return lat, loninterp_steps = {"0500M": 0.005,
"1000M": 0.01,
"2000M": 0.02,
"4000M": 0.04}
lat_S, lat_N = 0, 50
lon_W, lon_E = 80, 130
# 先乘1000取整是为了防止浮点数的精度误差累积
lat_S = int(1000 * lat_S)
lat_N = int(1000 * lat_N)
lon_W = int(1000 * lon_W)
lon_E = int(1000 * lon_E)
interp_steps = {x: int(y * 1000) for x, y in interp_steps.items()}
# 开始经纬度转行列号（内存充足情况下）
lc = {}# 保存各分辨率经纬度对应的行列号为字典
for resolution, interp_step in interp_steps.items():
lat = np.arange(lat_N, lat_S-1, -interp_step) / 1000
lon = np.arange(lon_W, lon_E+1, interp_step) / 1000
lon, lat = np.meshgrid(lon, lat)# 构造经纬度网格
line, column = latlon2linecolumn(lat, lon, resolution)
#l = l[:, :, newaxis]
#c = c[:, :, newaxis]
#lc = concatenate((l, c), axis=2)
FY-4A经纬度查找表
rawfile = '/home/mw/input/FY4A9106/FullMask_Grid_4000.raw'
dim = 2748
data = np.fromfile(rawfile,dtype=float,count=dim*dim*2)
latlon = np.reshape(data,(dim, dim, 2))
lat = latlon[:,:,0]
lon = latlon[:,:,1]
因为.raw是全圆盘，周围四个角上都是填充值9999，中间的数据是正确的
风云数据服务网可以下载到二进制格式的数据，为了方便读取，这里提供了转为TIFF格式的查找表
百度云链接：：https://pan.baidu.com/s/1ttdraWVabmClSkbK2n90sg 密码：6sro
来源：https://www.ixxin.cn/2018/08/11/fy4ageoluttiff/
<span style="background-color: rgb(255, 255, 255);">dataset = gdal.Open('/home/mw/input/FY4A9106/FullMask_Grid_4000.tif')</span>上述两种方法都可以读取FY-4A经纬度查找表，这里比较一下两种方法所得的经纬度数据，可以发现是完全一致的。
(lat==lat1).all()True
(lon==lon1).all()True
FY-4A更多应用
定量降水估计
风云四号A星降水估计2021年7月24日08:00时
三维立体扫描
台风内部大气温度三维分布（2021年7月24日11:00时）
综合灾害风险
台风登陆前后，浙江北部、上海、江苏南部等地具有风雨灾害风险，特别是浙江东北部沿海局地有中等及以上风险。
（来源：国家卫星气象中心遥感应用服务中心） 网址：http://www.zgqxb.com.cn/zx/yw/202107/t20210724_581591.html
脚本获取
在好奇心Log公众号后台回复卫星处理，欢迎大家点击阅读原文去线上直接练习。
文章来源于微信公众号：好奇心Log
页:
[1]