本教程首发于极思课堂,极思课堂(极思课堂 - GeoScene Online)是易智瑞信息技术有限公司面向公司客户、合作伙伴、GIS从业者、高校教师与学生以及热衷于GIS技术的极客们打造的一站式GIS前沿技术学习平台。欢迎大家访问、学习与关注。
课程原文链接:极思课堂 - GeoScene Online
https://www.geosceneonline.cn/learn/lesson/a1a37fbf2f614f40a5192f5d53f29e6e
1 引言
台风,这一自然界中极具破坏力的天气现象,是热带气旋的一种强烈形式,主要在热带或亚热带海域生成。当海洋表面的水温超过26.5摄氏度,且大气中的水汽充足、垂直风切变较小时,便为台风的孕育提供了理想的条件。台风的形成起始于一个微弱的热带扰动,随着周围环境能量的输入,这个扰动逐渐增强,旋转加速,中心气压下降,形成一个巨大的旋转风暴系统。
台风登陆时,会带来狂风、暴雨、洪水和风暴潮等灾害,对沿海地区的基础设施、农业生产、居民生活乃至生命安全构成严重威胁。强风可以摧毁建筑物,拔起树木,影响交通;暴雨则易引发洪水,淹没低洼地带;风暴潮则能推高海水,造成海水倒灌,对海岸线造成侵蚀。因此,台风预警和防御工作至关重要,包括提前疏散沿海居民、关闭危险区域、加强堤防和排水系统等措施,以减轻台风带来的损失。
2018年第22号台风“山竹” - 实时路径可视化
以2018年的第22号台风—山竹(英语:Super Typhoon Mangkhut)为例,最开始形成时只是热带扰动,一步一步升级为热带低压、热带风暴,并最终形成为超强台风,并于广东台山海宴镇登陆,造成了巨大的人员财产损失。据 百度百科 资料,截至2018年9月18日17时,台风“山竹”已造成广东、广西、海南、湖南、贵州5省(区)近300万人受灾,5人死亡,1人失踪,160.1万人紧急避险转移和安置。
超强台风—山竹的发展过程
- 2018年9月4日:一个低压区在国际换日线以西海域形成,美国海军研究实验室给予热带扰动编号99W。
- 2018年9月7日:上午8时,日本气象厅将其升格为热带低压,并在同日下午12时05分对其发布烈风警报;晚间9时15分,日本气象厅将其升格为热带风暴,给予国际编号1822,并命名为“山竹”。
- 2018年9月8日:凌晨5时和晚上8时,联合台风警报中心和中国国家气象中心分别将其升格为热带风暴。
- 2018年9月11日:中国国家气象中心和韩国气象厅率先将其升格为超强台风。
山竹 9月11日风云FY-4A卫星云图
- 2018年9月15日:凌晨1时40分,台风"山竹"从菲律宾北部登陆,接近中午时已经离开陆地,以每小时25公里速度吹向南海。
山竹 9月15日风云FY-4A卫星云图
- 2018年9月16日:17时,第22号台风“山竹”(强台风级)在广东台山海宴镇登陆,登陆时中心附近最大风力14级(45米/秒,相当于162公里/小时),中心最低气压955百帕。
山竹 9月16日风云FY-4A卫星云图
- 2018年9月17日:台风“山竹”已减弱为热带低压,到达广西百色市境内,强度继续减弱并远离广东省。
山竹 9月17日风云FY-4A卫星云图
2 课程简介
2018年的第22号台风“山竹”,威力远超在日本登陆的“飞燕”,中心最高风力可能超过17级,是同时期太平洋和大西洋上当之无愧的风王。
对于这种极端气候,本课程详细介绍了如何通过下载和处理风云4号(FY-4A)气象卫星云图数据和台风路径矢量数据,使用GeoScene Pro桌面端软件直观地对台风路径进行可视化,方便我们演示台风的运动轨迹及开展后续的科学研究工作。
台风“山竹”路径动态模拟-2018年09月15日
3 软件及数据
3.1 软件及版本要求
桌面端GIS软件:GeoScene Pro (2.1及以上版本)
3.2 数据来源
3.2.1 风云4号(FY-4A)气象卫星
风云四号(FY-4A)气象卫星是我国第二代静止气象卫星。作为新一代静止轨道定量遥感气象卫星,卫星的辐射成像通道由FY-2G星的5个增加为14个,星上辐射定标精度0.5 K、灵敏度0.2 K、可见光空间分辨率0.5km,与欧美第三代静止轨道气象卫星水平相当,其数据稳定便于获取,同时成像质量好,有利于对于气象的观测。
先进的静止轨道辐射成像仪(Advanced Geostationary Radiation Imager,AGRI)是风云四号静止气象卫星的主要载荷之一,通过精密的双扫描镜机构实现精确和灵活的二维指向,可实现分钟级的区域快速扫描;采用离轴三反主光学系统,高频次获取超过14个波段的地球云图,并利用星上黑体进行高频次红外定标,以确保观测数据的精度。
先进的静止轨道辐射成像仪主要承担获取云图的任务,至少14个通道,是风云二号5通道的近三倍,在风云二号观测云、水汽、植被、地表的基础上,还具备了捕捉气溶胶、雪的能力,并且能清晰区分云的不同相态和高、中层水汽。相比于风云二号单一可见光通道的限制,风云四号首次制作出彩色卫星云图,最快1分钟生成一次区域观测图像。
风云系列卫星的数据可以从 风云卫星遥感数据服务网 下载。
本课程中使用FY-4A卫星AGRI成像仪全圆盘4KML1数据产品,格式为HDF。选择的时间范围为2018年9月10日-9月18日之间的数据。
本课程使用的FY-4A AGRI L1级数据
本课程使用的成像仪全圆盘4KML1数据
3.2.2 台风路径矢量数据
台风路径数据集【1945-2022】 记录了1945年至2022年期间西太平洋范围内的台风路径信息,包括台风编号、名称、登陆地点、强度、路径等多个指标。
4 风云卫星数据预处理
4.1 全圆盘数据坐标转换:标称投影转换为WGS 1984投影
下载的成像仪全圆盘4KML1数据,采用的是标称投影。标称网格是由位于理想位置的地球同步轨道卫星观测视图计算标称网格坐标与地理经纬度间的一种投影变换。国内外大多数静止轨道气象卫星,如风云、ElectroL、H8 AHI、MTG FCI等近几年发射的或预计发射的卫星,均采用CGMS LRIT/HRIT全球规范中定义的标准投影生成L1标称图像。
风云标称影像示例
GeoScene Pro不支持直接导入全圆盘标称影像,所以首先需要对影像做投影转换。可借助 FY4圆盘转等经纬工具 ,将原始的全圆盘标称影像(hdf格式)转换为WGS 1984坐标下的影像(nc格式)。
影像坐标转换过程
最后得到WGS 1984坐标系下的风云影像(nc格式)在GeoScene Pro中显示效果如下:
坐标转换后的风云影像在GeoScene Pro中的真彩色合成效果
4.2 波段计算:影像色彩校正
FY-4A AGRI L1数据是风云四号A星先进的静止轨道辐射成像仪的一级数据产品。该传感器1-3波段接近于可见光的蓝色、绿色和红色通道,但它们并不完全匹配对应的可见光波长,见下表:
AGRI波段波长与可见光波长对照表
所以直接使用卫星的1-3波段进行真彩色合成会导致合成的图像针对植被的部分偏红,影响最终的可视化效果,因此需要矫正。我们新建了三个波段,命名为“ComputedRed”、“ComputedGreen”、“ComputedBlue”,分别表示校正后的红绿蓝波段,这三个波段使用原始影像中的“NOMChannel101”(波段1-蓝)、“NOMChannel102”(波段2-绿)、“NOMChannel103”(波段3-红)波段经过计算得到。
同时由于1景图像中有14个波段,本课程中仅需前3个波段进行真彩色合成,所以在代码中还对不需要的波段进行删除处理,缩减影像大小。
详细代码见下:
import os
import xarray as xr
# 输入文件夹路径
input_folder = r'E:\台风山竹\FY4A_201809_nc_投影变换'
# 输出文件夹路径
output_folder = r'E:\台风山竹\FY4A_201809_nc_色彩校正'
# 确保输出文件夹存在
if not os.path.exists(output_folder):
os.makedirs(output_folder)
# 遍历输入文件夹中的所有文件
for filename in os.listdir(input_folder):
if filename.endswith('.NC'):
# 构建完整的文件路径
nc_path = os.path.join(input_folder, filename)
# 打开 NetCDF 文件
ds = xr.open_dataset(nc_path)
# 复制数据集并应用计算
#增加三个波段“ComputedRed”,“ComputedGreen”,“ComputedBlue”
#使用原有的RGB波段("NOMChannel03","NOMChannel02","NOMChannel01")进行计算,计算公式见下面
#使用计算后的波段能改善最终的显示效果
ds1 = ds.copy()
ds1 = ds1.assign(ComputedRed=(ds1["NOMChannel02"] - 0.14 * ds1["NOMChannel03"]) / 0.86)
ds1 = ds1.assign(ComputedGreen=ds1["NOMChannel01"] * 0.6 + ds1["NOMChannel02"] * 0.33 + ds1["NOMChannel03"] * 0.07)
ds1 = ds1.assign(ComputedBlue=ds1["NOMChannel01"])
#写入新波段的属性信息
ds1["ComputedRed"].attrs = {"long_name":"Computed Red Channel", "units":"K", "coordinates":"latitude longitude", "formula":"('NOMChannel02' - 0.14*'NOMChannel03')/0.86"}
ds1["ComputedGreen"].attrs = {"long_name":"Computed Green Channel", "units":"K", "coordinates":"latitude longitude", "formula":"'NOMChannel01'*0.6 + 'NOMChannel02'*0.33 + 'NOMChannel03'*0.07"}
ds1["ComputedBlue"].attrs = {"long_name":"Computed Blue Channel", "units":"K", "coordinates":"latitude longitude", "formula":"'NOMChannel01'"}
# 删除原始通道变量
ds2 = ds1.drop_vars([
"NOMChannel01", "NOMChannel02", "NOMChannel03", "NOMChannel04",
"NOMChannel05", "NOMChannel06", "NOMChannel07", "NOMChannel08",
"NOMChannel09", "NOMChannel10", "NOMChannel11", "NOMChannel12",
"NOMChannel13", "NOMChannel14"
])
# 构建输出文件路径
output_path = os.path.join(output_folder, filename)
# 保存处理后的数据集到新的 NetCDF 文件
ds2.to_netcdf(output_path)
print(f"Processed and saved {filename} to {output_path}")
print("All files processed successfully.")
色彩校正和波段删减处理代码
处理前后影像的对比结果如下图所示,可以看到经过色彩校正处理后影像植被区域可视化效果较好。
处理前影像真彩色合成图(左图)与处理后影像真彩色合成图效果对比(右图)
5 制作卫星云图时序影像
GeoScene Pro内置时间滑块,可针对带有时间属性的要素或影像图层按照时间的先后顺序进行依次浏览及按时间播放。
按时间顺序浏览卫星云图可清晰的查看台风运动的全过程,为了在GeoScene Pro中实现对时序云图数据的浏览与可视化,我们需要将上一步预处理完成的多维nc数据构建成为多维镶嵌数据集。
5.1 构建镶嵌数据集
在GeoScene Pro中可通过构建镶嵌数据集并加载nc数据的方式对影像进行管理、可视化与分析。
步骤如下:打开GeoScene Pro,在右侧的目录视图中,右键地理数据库,在弹出的菜单中选择“新建”,创建“镶嵌数据集”,选择“GCS_WGS_1984”坐标系进行创建。
创建镶嵌数据集
之后右击生成的镶嵌数据集,点击“添加栅格数据”,将打开“添加栅格至镶嵌数据集”工具。在工具面板中首先选择“输入数据”下的“文件夹”,并选择存放数据的文件夹。接下来,修改栅格类型为“NetCDF”,处理模板选择“Multiband Composite”,选择该处理模板可直接生成RGB合成影像。
添加栅格至镶嵌数据集工具
接下来选择“栅格类型”旁“属性”按钮,在弹出的栅格类型属性框中选择“变量”标签页,在右侧的复选框中选择“ComputedBlue”、“ComputedGreen”、“ComputedRed”这三个波段,点击确认和运行。
使用计算后的RGB波段组成彩色影像
最后使用计算统计数据工具计算栅格数据集的统计数据。
计算统计数据工具
计算统计数据工具运行完成后,可以看到创建的镶嵌数据集能正确加载风云卫星影像(nc格式数据)并将影像中的1-3波段合成显示为真彩色影像。通过查看覆盖区统计表,可以看到全部54景影像已全部加载至镶嵌数据集中。
使用镶嵌数据集对风云卫星影像进行管理与可视化
5.2 构建时间字段
现在我们构建了包含风云卫星影像的镶嵌数据集,接下来需要在镶嵌数据集中增加时间维度,才能进一步构建为包含时间的多维镶嵌数据集,并在Geoscene Pro中使用时间滑块来显示随时间变化的数据。
打开覆盖区的属性表,新建一个名为“BeijingTime”,类型为日期的字段,用于存储每一景影像对应的北京时间。
新建名为“BeijingTime”的时间字段
接下来,使用字段计算器对“BeijingTime”字段进行赋值。在本例中,每景影像的获取时间已提现在影像名称中,所以可以先从“Name”字段中提取时间文本,然后将时间文本转换为日期类型,其中需要注意的是源数据使用的时间为UTC标准时,最后需要转换为北京时间(UTC+8)。
为时间字段赋值
计算后的“BeijingTime”字段如下所示,已正确显示了每一景影像对应的获取时间(北京时间):
计算后的时间字段
5.3 构建多维镶嵌数据集
接下来需要将镶嵌数据集转化为多维镶嵌数据集(在本例中可以理解为增加了时间维度)。步骤如下:
在目录视图中,鼠标放在镶嵌数据集上并点击右键,选择“修改”—“构建多维信息”。
为镶嵌数据集构建多维信息
在打开的“构建多维信息”工具窗口中,变量名称选择“computedred”,维度字段选择上一步创建的北京时间“BeijingTime”字段。点击运行。
构建多维信息工具
查看镶嵌数据集的图层属性,在“源”标签中的“多维信息”栏中,可以看到“StdTime=54”的提示,表示镶嵌数据集已成功转换为多维镶嵌数据集,且具有54个不同的时间。
查看多维镶嵌数据集的多维信息
5.4 启用图层时间属性
包含时间属性的镶嵌数据集构建完成后,如想在GeoScene Pro中按时间查看不同时期的卫星云图,需要为图层配置时间。
同样在图层属性窗口中,点击“时间”标签,在右侧选择“根据属性值过滤图层内容”,时间字段选择“Standard Time”字段,软件会自动计算时间范围,点击确定即可启用图层的时间属性。
启用图层的时间属性
图层启用时间后,将会显示时间滑块,可以按时间播放不同时间的卫星云图。
按时间播放不同时间的卫星云图
6 台风矢量数据处理
在卫星云图时序影像制作完成后,接下来我们需要对台风矢量数据进行处理,具体包括绘制台风点位、台风行进路径,台风风圈、24小时/48小时警戒线等一系列重要的台风关键要素,这些要素制作完成后与台风影像一起显示,能更加清晰的展示台风登陆的全过程。
6.1 绘制台风点位
台风点位数据中包含了台风当前点位的经纬度、台风强度、台风风力、台风速度、台风移动方向、台风移动速度、台风气压、台风7级/10级/12级风圈半径等一系列重要的台风信息,是台风研究中的重要依照。
在Geoscene Pro 中,使用添加数据中的XY点数据,依据数据中的lng、lat字段生成含有台风轨迹数据的点位数据。在显示的XY表转点对话框内,输入表为下载的台风“山竹”的轨迹数据(csv文件),设置输出要素类,X字段为lng,Y字段lat,坐标系为GCS_WGS_1984,如下图所示。
使用XY表转点工具在GeoScene Pro中绘制台风“山竹”的点位
接下来,需要对生成的台风点位,根据台风强度(strong字段)进行唯一值符号化,目的是用不同颜色代表不同台风等级,具体使用的等级样式如下:
不同台风等级的颜色对照表
最后得到台风“山竹”点位的显示效果如下:
在GeoScene Pro中绘制台风点位(颜色表示台风等级)
6.2 绘制台风路径
接下来我们需要根据台风点位生成表示台风路径的线要素。
使用GeoScene Pro中的“点集转线”工具,可依据原有的台风点位数据生成相应的台风路径数据。
点集转线工具
最后对生成的台风路径要素进行符号化处理,样式依据相邻点位的台风等级进行设置,得到效果如下图显示:
在GeoScene Pro中绘制台风路径(颜色表示台风等级)
6.3 绘制台风风圈
台风风圈半径是指台风中心周围不同强度风力所覆盖的区域的半径。具体来说,它是指从台风中心向外,风速达到一定标准的区域的边界距离。台风风圈半径通常用来描述台风的规模和强度,以及其可能带来的影响范围。
台风7级、10级、12级风圈示意图
在气象学中,台风风圈的方向角度范围通常被分为四个象限,分别是东北(NE)、东南(SE)、西南(SW)和西北(NW)。每个象限代表的是围绕台风中心的90度范围。具体来说:
- 东北(NE)象限:指的是从正北方向顺时针旋转到正东方向之间的90度范围,即角度范围0度到90度。
- 东南(SE)象限:指的是从正东方向顺时针旋转到正南方向之间的90度范围,即角度范围90度到180度。
- 西南(SW)象限:指的是从正南方向顺时针旋转到正西方向之间的90度范围,即角度范围180度到270度。
- 西北(NW)象限:指的是从正西方向顺时针旋转到正北方向之间的90度范围,即角度范围270度到360度。
在台风点位数据中,包含了“radius7_qu”、“radius10_qu”、“radius12_qu”三个字段,这三个字段分别表示台风在当前时间点,其7级、10级和12级风圈的半径大小,即在最大风速半径外,近地面各个方向风速衰减至 17.2 米/秒、24.5 米/秒 以及 32.7 米/秒 时离台风中心的距离。
具体来看,例如以 “radius7_qu” 字段中的第一行 {'ne': 250, 'se': 150, 'sw': 150, 'nw': 150} 为例,“ne”表示方向,后面的数字表示该方向的半径大小,即台风在该点位,7级风圈的半径在各个方向依次是:
- 东北(ne):250公里
- 东南(se):150公里
- 西南(sw):150公里
- 西北(nw):150公里
台风“山竹” - 7级、10级、12级风圈大小
在明确了各字段的含义之后,我们可以根据上面的信息分别绘制台风山竹的7级、10级、12级风圈。
在这里,我们以制作12级风圈范围为例。点击台风“山竹”点位图层,打开属性表,为每个方向新建四个字段。例如,在下图中,我们首先为东北方向新建了四个字段:“wind_12u_ne_deg_min”、“wind_12u_ne_deg_max”、“wind_12u_ne_radius_min”、“wind_12u_ne_radius_max”分别代表东北方向扇形区域的起始角、终止角、半径最小距离、半径最大距离。这些字段后续配合“根据要素生成扇形视域工具”能自动生成扇形区域(台风的1/4风圈)。
为东北方向新建4个字段
接下来,使用“字段计算器”,为角度字段“wind_12u_ne_deg_min”、“wind_12u_ne_deg_max”进行赋值,赋值依据如下:东北(NE):分别赋值0和90。
然后对“wind_12u_ne_radius_min”、“wind_12u_ne_radius_max”字段进行赋值,扇形区域半径最小值统一赋值为0,而扇形区域半径最大值依据“radius12_qu”字段,读取对应方向的风圈半径。如下图就是读取ne方向的风圈半径。
读取风圈半径
最后将所有字段都赋值之后的结果如下图:
接下来,可以使用“根据要素生成扇形视域”工具,选择对应的字段类型,生成12级风圈东北方向的风圈范围。
根据要素生成扇形视域工具
最后得到的结果如下,成功绘制了台风“山竹”在各点位东北方向12级风圈范围。
台风“山竹”在各点位东北方向12级风圈的范围
依据前文,重复三次,分别生成其他方向的风圈范围,叠加显示如下:
台风“山竹”在各点位各方向12级风圈的范围
接下来使用“合并”地理处理工具将不同方向的要素进行合并。再使用“融合”地理处理工具,根据特定属性将多个要素合并成一个要素。为了更直观地展示结果,调整要素的符号系统。这包括设置颜色、线型、填充等,以便清晰地区分不同的风圈范围。得到各台风点位12级风圈最终效果如图:
台风“山竹”在各点位各方向12级风圈的范围(最终效果)
最后使用“连接字段”地理处理工具,将原台风点位的属性信息转移到新的台风风圈要素中。
至此,台风“山竹”12级风圈范围的制作完毕,重复上述步骤,可以制作7级和10级风圈范围。
台风“山竹”在各点位各方向7级、10级、12级风圈的范围
接下来使用“合并”地理处理工具,可以将不同等级的风圈要素合并为一个图层。针对时间,新建一个北京时间字段,并赋值。(注意:本例中源数据中的时间字段为北京时间,故不需要额外进行时区转换)
最后,点击GeoScene Pro左侧内容视图中的图层,点击打开属性,开启图层中的时间属性。
开启台风“山竹”风圈图层的时间属性
开启时间属性之后,点击上方的时间面板,可以编辑时间轴进行对于图层的播放,叠加风云云图,最后播放效果如下图所示:
台风“山竹”风圈播放展示
6.4 绘制台风24小时/48小时警戒线
台风24/48小时警戒线是台风追踪过程中的常用指标,含义为跨过警戒线的台风会在24/48小时内登陆。进入24小时警戒线以后,气象部门会每个小时对台风进行监测,密切注视台风动向,并随时通报台风的实际情况。
在GeoScene Pro中绘制24/48小时台风警戒线折点
国家规定的警戒线其转折点的经纬度如下:
- 24小时警戒线折点经纬度:[34, 127], [22, 127], [18, 119], [11, 119], [4.5, 113], [0, 105]
- 48小时警戒线折点经纬度:[34, 132], [15, 132], [0, 120], [0, 105]
使用“添加数据”中的“XY点数据”选项在Geoscene Pro中绘制24/48小时台风警戒线折点。
接下来,调用GeoScene Pro中“点集转线”工具,将点数据转换为线数据,同时进行符号化及标注,最后得到的台风24小时/48小时警戒线效果如下所示。
在GeoScene Pro中绘制24/48小时台风警戒线
7 台风—“山竹”登陆过程模拟
将制作好台风矢量数据叠加到卫星云图时序影像上,并点击GeoScene Pro上方的“时间”窗格,设置播放的时间间隔与时间范围,随后点击时间轴的播放按钮,便可播放台风“山竹”的登陆过程。可清晰地通过风云卫星云图看到台风“山竹”在每个时间点的位置,以及对应的7级、10级、12级风圈的大致范围,整个过程清晰流畅。
下面,我们节选台风“山竹”行进过程中,9月15日、16日、17日的关键时间点进行模拟演示。
- 2018年9月15日
台风"山竹"从菲律宾北部登陆,接近中午时已经离开陆地,以每小时25公里速度吹向南海。
台风"山竹"行进路径 (2018年9月15日 11:00 - 17:00)
- 2018年9月16日
17时,第22号台风“山竹”(强台风级)在广东台山海宴镇登陆,登陆时中心附近最大风力14级(45米/秒,相当于162公里/小时),中心最低气压955百帕。
台风"山竹"行进路径 (2018年9月16日 11:00 - 17:00)
- 2018年9月17日
台风“山竹”已减弱为热带低压,到达广西百色市境内,强度继续减弱并远离广东省。
台风"山竹"行进路径 (2018年9月17日 11:00 - 17:00)
8 课程小结
在本课程中,我们展示了如何利用 GeoScene Pro 软件对气象数据进行可视化处理,特别是以超强台风“山竹”为例,详细讲解了如何利用风云气象卫星云图影像结合台风矢量数据,进行台风登陆过程的实时模拟。我们从数据下载、预处理到使用软件进行管理和可视化,最终实现交互式地图展示。课程重点包括风云卫星云图和台风矢量数据的处理,台风点位、路径、风圈半径的制作,以及台风警戒线的绘制,最终通过模拟台风登陆过程,为决策提供科学依据。
在下一节课程中,我们将讲解如何在Web端使用GeoScene Enterprise 或 GeoScene Online中的 Experience Builder 0代码快速搭建一个台风监测大屏应用。欢迎大家继续学习!
2018年第22号台风 - 山竹 Mangkhut 监测大屏