绘图工具库cedarkit-maps介绍

目录

本文来自笔者 2024 年 3 月 29 日在单位部门内部分享的《cedarkit-maps工具库介绍》。

声明:本文仅代表笔者个人观点,所选材料截止时间为 2024 年 3 月,不代表当前的业务系统现状。 如果想了解最新信息,请咨询 CEMC 相关部门。

本文主要介绍正在开发中的 cedarkit-maps 绘图工具,包括开发背景、设计、实现与当前开发成果。

背景

业务系统现状

绘图任务通常分为准备数据和图形绘制两个步骤。当前业务系统在准备数据方面,简单的绘图任务直接使用模式数据,复杂的绘图任务需要计算得到绘图需要的最终数据。比如 CMA-MESO 使用 NCL 脚本实现降水相态计算,CMA-GEPS 使用内部开发的 Fortarn 程序计算概率。在图形绘制方面主要使用 NCL 和 GrADS 绘图,少部分图片产品使用 Python 绘制。

图 CEMC业务系统绘图任务现状,按功能划分

从工具角度来看,绘图需要使用基础的可视化绘图工具库,例如 NCL、GrADS、Python等。CEMC 使用基础绘图库封装了一系列绘图包,例如 NCL 的 ncllib 和 GrADS 的 map.gs 等脚本,并使用绘图工具库和封装工具库为图片产品开发绘图脚本。同时业务系统的绘图任务需要将业务系统数据与图片产品脚本相关联,通常通过 Shell 脚本调用绘图脚本的方式实现。另外也使用 wgrib2、eccodes 等开源数据处理工具,并为特定的任务开发 Fortran 数据处理程序。绘图脚本中也包含了业务系统的部分信息,比如数据文件命名规范、数据文件路径等。从下图可以看到,业务任务、图片产品、绘图工具库三层之间没有形成明显的边界,导致业务逻辑和绘图实现耦合在一起,影响绘图脚本的适用性。

图 CEMC业务系统绘图任务现状,按工具划分

NCL

NCL 主要用于 CMA-GFS、CMA-MESO,以及 CMA-TYM、CMA-GEPS、CMA-REPS 的部分图形。

CEMC 基于 NCL 开发了 ncllib 库,提供底图、要素查找、文件路径等功能函数,应用在 GFS 和 MESO 的 大部分 NCL 绘图脚本中。ncllib 内置了 NMIC 发布的中国和世界地图底图 shapefile 数据文件。ncllib 库的典型调用如下所示,代码中加载的 ncllib 脚本说明见下表。

load "$GRAPHIC_PRODUCT_LIB_ROOT/lib/common.ncl"

load "$GRAPHIC_PRODUCT_LIB_ROOT/lib/GFS/GFS_GRAPES.ncl"

load "$GRAPHIC_PRODUCT_LIB_ROOT/lib/GFS/$FORECAST_DATA_FORMAT/data_file_name.ncl"

load "$GRAPHIC_PRODUCT_LIB_ROOT/lib/ncl_variable/$FORECAST_DATA_FORMAT/$FORECAST_DATA_CENTER/ncl_variable.ncl"

load "$GRAPHIC_PRODUCT_LIB_ROOT/lib/cn_map_function.ncl"

表 ncllib库常用加载文件说明

文件名说明
cn_map_function_zhyt.ncl绘制地图
GFS/GFS_GRAPES.ncl设置绘图样式参数,输出图片路径
$FORECAST_DATA_FORMAT/$FORECAST_DATA_CENTER/ncl_variable.ncl加载变量
$FORECAST_DATA_FORMAT/data_file_name.ncl数据文件名
common.ncl图片文件名

CMA-GFS、CMA-MESO 的大部分图形使用 NCL 脚本绘制绘制,并使用了 CEMC 封装的 ncllib 库。为了保证各个系统的独立性,每个业务后处理系统都各自包含一份 ncllib 库代码,各个绘图脚本会调用 ncllib 库中的函数实现通用功能。

图 NCL绘图示例,CMA-GFS 东亚区域 500hPa 高度场 + 海平面气压,来源:数值预报_CMA全球天气模式_东亚(10N-60N,60E-140E)_500hpa高度场+海平面气压场

其他业务系统也有大量使用 NCL 脚本的绘图,这类绘图脚本一般由多个相关联的 NCL 脚本组成,通常作为一个独立的绘图包使用。绘图脚本没有使用 ncllib 类似的共享函数库,而是在绘图包的 NCL 脚本中简单实现类似功能。 比如下图展示了 CMA-GEPS 绘制的台风路径图形,包括 CMA-GFS、CMA-TYM、CMA-GEPS 成员与集合平均等多个模式的台风路径预报结果。该绘图包中的 function_geps.ncl 中包含了计算和绘图相关函数,比如计算均值,绘制地图、图例等。

图 NCL绘图示例,CMA-GEPS 台风路径,来源:数值预报_CMA热带气旋预报_CMA模式集合预报_路径

GrADS

GrADS 主要用于 CMA-GEPS、CMA-REPS 和 CMA-TYM。

各个业务系统也基于GrADS封装了一些绘图工具库,比如 CMA-GEPS 的 Gslib 库:

名称说明
basemap.gs底图
map.gs底图
define_colors.gs色表

下图是 CMA-GEPS 使用 GrADS 绘制的概率预报图形。

图 GrADS绘图示例,CMA-REPS 全国 24h 累积降水 概率预报 >= 10.0mm,来源:数值预报_CMA区域集合模式_全国_24h累积降水_概率预报_>=10.0mm

下图是 CMA-REPS 使用 GrADS 绘制的邮票图。

图 GrADS绘图示例,CMA-REPS 邮票图 24 小时累积降水,来源:数值预报_CMA全球集合模式_邮票图_24小时累积降水

下图是 CMA-TYM 使用 GrADS 绘制的台风路径预报图。

图 GrADS绘图示例,CMA-TYM 台风路径,来源:数值预报_CMA台风模式_CMA_TYM_路径 (nmc.cn)

下图是 CMA-MESO 使用 GrADS 绘制的站点时序高度剖面图。

图 GrADS绘图示例,CMA-MESO 站点,相对湿度、温度、风,北京,来源:数值预报_CMA区域模式_单站_相对湿度、温度、风_北京

Python

少量图形使用 Python 绘制。

CMA-GFS、CMA-MESO 和 CMA-TYM 业务系统中使用 Python 绘制风云卫星云图模拟图形,替代之前的 IDL 版本。

注:对于运行在超算平台的数值预报业务系统来说,使用 IDL 等需要授权的商用软件的任务通常不具备良好的迁移性,应当尽量避免业务系统依赖商用软件。

图 Python绘图示例,CMA-GFS 全球 FY4A 卫星云图模拟红外亮温,来源:数值预报_CMA全球天气模式_全球_FY4A卫星云图模拟_红外亮温 (nmc.cn)

上图使用 Nio 读取 GRIB2 数据,使用 Matplotlib 绘图,使用 Basemap 绘制地图。

Nio 打开 GRIB2 文件,并加载变量

filename = 'gmf.gra.' + bdate.strftime('%Y%m%d%H') + f'{ft:03d}.grb2'
fili = Nio.open_file(os.path.join(data_path,filename), "r")
wv = fili.variables["VAR_0_4_227_P0_L1_GLL0"][:]

Basemap 绘制全球地图

map = Basemap(
	projection='cyl',
	llcrnrlon=0,
	llcrnrlat=-90,
	urcrnrlon=360,
	urcrnrlat=90,
	resolution='l',
	ax=ax1
)
map.drawcoastlines()
map.drawparallels(np.arange(-90, 90, 30), labels=[1, 0, 0, 0], fontsize=8)
map.drawmeridians(np.arange(0, 360, 30), labels=[0, 0, 0, 1], fontsize=8)

2023 年,CEMC 面向杭州亚运会开发了区域 3 公里集合预报系统(CMA-HZ-CAEF),并成功实现在国产超算平台实时运行。CEMC 在该系统首次全面使用 Python 绘制业务图片产品,并经过长时间的实时运行已证明业务系统可以使用 Python 绘图来代替原有的 GrADS 和 NCL 的绘图脚本。该绘图工具库将进一步应用在即将业务化升级的中国区域 3KM 集合预报系统中。

CMA-HZ-CAEF 绘图库采用多项优秀的设计理念,提高工具库的可维护性:

  • 数据处理与绘图分离data_process.py 包含多个数据计算函数,grib2_process.py 提供处理 GRIB2 数据的类 Grib2Processor。绘图核心算法则保存在 plot_process.py 模块。
  • 配置与算法分离。使用专用模块 global_var.py 保存配置信息,包括数据范围、格点数、分辨率、集合成员数、要素名称、单位等数据信息,以及标题字号、层次、颜色表、颜色条位置等图形样式。应用到其他系统时仅需要修改该配置文件,而不用修改算法文件本身。
  • 提供命令行调用接口。通过 main.py 提供面向用户的命令行接口,并在 config.cfg 提供任务级的配置信息,定义绘制图片的种类。

图 CMA-HZ-CAEF 集合平均图片产品脚本函数调用关系图,仅供参考

图 CMA-HZ-CAEF 海平面气压集合平均图片,来源:系统运行目录

图 CMA-HZ-CAEF 24 小时降水邮票图,其中 m3 和 m10 两个成员没有数据,数据缺失通常因为积分溢出。来源:系统运行目录

其他绘图

另有一些图形使用 Meteoinfo、IDL 等工具绘制。

当前存在的问题

1. 多种绘图工具导致绘图风格不一致,缺少统一绘图工具

每个业务系统都各自维护一套绘图脚本,同样的图形使用不同的脚本绘制,造成代码重复,也增加了维护成本。

例如下图展示的 CMA-GFS 和 CMA-REPS 的 K 指数图片分别使用 NCL 和 GrADS 绘制,可以看到两者从布局到样式均有所不同。

因为 CEMC 长久以来都采用按数值模式切分系统的管理模式,绘图脚本分散在多个部门,再加上 NCEP 已经放弃 NCL 开发,业务系统中从 GrADS 到 NCL 的替代过程实际上已经停滞。 另外,因为 NCL 较高的学习成本,业务系统使用的 ncllib 工具库也没有积极在内部推广使用,CEMC 研发工作仍大量使用 GrADS。 这就造成了研发与业务之间在数据可视化上的不一致,研发人员基本无法自行绘制业务图片产品,从而要求业务系统归档图片产品,而业务工具开发人员无法从更广泛的用户中获得正向反馈,从而失去继续更新完善工具的动力。

需要构建一套统一的绘图工具库,对于某种图片产品,所有系统使用同一个绘图脚本,绘图脚本可以适配多个数值模式数据,以保证各系统绘图风格的统一。随着气象大模型的火爆发展, Python 在气象领域变得越来越流行。 CEMC 业务系统中的少部分图片使用 Python 开源库 Matplotlib 绘制,GetPy 也是基于 Python 开发,CEMC 产品团队也将以 Python 构建产品系统,新版本全球和区域集合预报系统也将采用 Python 构建产品模块和制作流程,研发人员也逐渐开始用 Python 进行绘图。 所以,可以再一次尝试统一绘图工具,面向研发和业务需求开发一款基于 Python 的绘图工具包,从而实现绘图风格的统一,有望能够更广泛地应用到业务和研发工作中。

图 CMA-GFS东亚区域 K 指数,使用 NCL 脚本绘制,来自:数值预报_CMA全球天气模式_东亚(10N-60N,60E-140E)_K指数

图 CMA-REPS 中国区域 K 指数集合平均和离散度,使用 GrADS 绘制,来自:数值预报_CMA区域集合模式_全国_K指数_集合平均和离散度

2. 依赖本地文件,较难支持多种数据源

当前业务系统中的绘图脚本均使用本地文件作为数据源,部分绘图脚本将文件名硬编码到脚本中,导致脚本的适用性较差。

绘图脚本仅支持 HPC 归档目录的 CEMC 短文件名,如果文件名与绘图脚本中规定的文件名不同,需要提前执行重命名操作。例如天擎 (CMADaaS) 使用标准长文件名保存数据,在天擎平台中运行绘图脚本前需要将长文件名链接为短文件名,下表列出了两种文件名的对比。

表 CEMC 数值天气预报模式业务系统产品文件名

系统短文件名,HPC归档长文件名,天擎
CMA-GFSgmf.gra.YYYYMMDDHHFFF.grb2Z_NAFP_C_BABJ_YYYYMMDDHH0000_P_NWPC-GRAPES-GFS-GLB-FFF00.grib2
CMA-MESOrmf.hgra.YYYYMMDDHHFFF.grb2Z_NAFP_C_BABJ_YYYYMMDDHH0000_P_NWPC-GRAPES-3KM-ORIG-FFF00.grb2
CMA-TYMrmf.tcgra.YYYYMMDDHHFFF.grb2Z_NAFP_C_BABJ_YYYYMMDDHH0000_P_NWPC-GRAPES-TYM-ACWP-FFF00.grib2
CMA-GEPSgef.gra.MMM.YYYYMMDDHHFFF.grb2Z_NAFP_C_BABJ_YYYYMMDDHH0000_P_NWPC-GRAPES-GEPS-GLB-FFF00-mMMM.grib2
CMA-REPSmef.gra.MMM.YYYYMMDDHHFFF.grb2Z_NAFP_C_BABJ_YYYYMMDDHH0000_P_NWPC-GRAPES-REPS-CN-FFF00-mMMM.grib2

NCL 和 GrADS 也较难直接使用保存在内存中的数据源,比如天擎 MUSIC 接口检索的要素场,或者从 CEMC 格点数据平台 (GODAS) 检索的要素场。使用这类数据源时,通常需要将检索到的要素场保存到文件系统中,绘图脚本从文件系统中读取数据,相比直接使用内存中的要素场多了 2 次磁盘 IO 操作。

注:NCL 支持 OPeNDAP,但 CEMC 暂未使用 OPeNDAP 数据源。

需要将数据加载与图形绘制解耦,开发单独的数据加载工具库支持从多种数据源加载要素场,绘图工具库使用数据加载工具库生成的数据对象完成图像绘制。

3. 较难与Python生态集成

随着以气象大模型为代表的人工智能在气象领域的应用越来越广泛,Python逐渐成为气象工作者的热门编程语言。Python 领域已经形成完善的科学计算社区,维护了一系列成熟的开源项目,包括 Numpy、Scipy、Matplotlib 等项目,以及地学领域 PANGEO 社区支持的 Xarray、Jupyter、Dask、Iris 等。气象领域开源了大量 Python 项目,比如 ECMWF 的 cfgrib、earthkit 系列工具,英国气象局的 cartopy,新西兰气象局的 cylc,UCAR 的 geocat 系列工具。CEMC 也计划将产品制作系统的编程语言逐步统一为 Python,绘图工具库必须要考虑如何与 Python 科学计算生态的集成。

业务系统绘图当前使用的 NCL 和 GrADS 很难与 Python 集成,尽管有类似 PyNGL 和 PyGrADS 的封装接口,但没有与 Numpy 和 Xarray 等常用科学计算库集成,无法应用社区的开发成果。另外,NCAR 已不再维护 NCL,GrADS 更新也停留在 2019 年发布的 2.2.1 版本,远没有 Python 开源社区项目活跃。无论是开发 NCL 的 NCAR 还是开发 Magics++ 的 ECMWF,都陆续转向基于 Matplotlib 开发 Python 版本的绘图库,很有可能会放弃自有绘图库的开发维护。

为了更好地融入Python科学计算生态中,需要使用 Python 开发一套绘图工具库,利用 Numpy、Scipy、Xarray 提供数据结构和计算功能,基于广泛使用的开源项目 Matplotlib、Cartopy 实现绘图功能。

国内外发展

无论从 NCAR 的 NCL 还是从 ECMWF 的 Magics++ 来看,气象可视化工具正在经历从编译型语言过渡到脚本语言 Python 的趋势。笔者不负责任地推测,可能因为随着老一代程序员的退休,熟悉编译型语言的人才越来越难招,维护老工具的成本逐渐提高,而新一代程序员对 Python 更为熟悉,所以只能被迫转向用 Python 开发新的工具。

从 Python 项目的架构上来看,NCAR 和 ECMWF 都采用多项目组合的方式来替代原有单一工具。

NCAR 开发 GeoCAT 系列软件包,将 NCL 的功能拆分到多个项目中,形成一组既可以搭配使用又可以单独调用的工具包。

ECMWF 正在开发 earthkit 系列软件包,作为 Metview 的继任者,earthkit 将通过简化数据访问、分析和可视化,提供强大和方便的 Python 工具,以加速天气和气候科学工作流(参考《What next for Magics visualisation?》)。earthkit 工具栈同样由多个工具包构成:

同样国家气象中心 nmcdev 小组也开发了 nmc_met 系列软件包:

与分散多软件包对应的另外一种形式是单一集成软件包,将所有功能整合到同一个包中。 比如 unidata 的 MetPy ,dtcenter 的 METPlus ,Met Office 的 iris

nmcdev 也有两个整合面向检验诊断任务的工具包 metdig 和 meteva,其中 CEMC 已将 metdig 应用在评估室开发的天气学分析系统网站中。

  • metdig:天气学诊断分析工具(仅部分 IO 功能依赖 nmc_met_io 库)
  • meteva:全流程检验程序库(全内置)

图 metdig 绘图示例,CMA-GFS 2米温度填充图,来源:天气学分析系统

综上所述,无论从 CEMC 内部科研业务需求,还是从国内外气象领域的相关工作来看,开发基于 Python 的绘图工具已成为一种趋势。为此 cemc-oper 小组正在开展 cedarkit 项目,针对气象数据可视化需求开发一款专用的工具库 cedarkit-maps。

简介

本文介绍 Python 工具栈 cedarkit 中的 cedarkit-maps 绘图工具库,该工具用于对气象数据进行可视化绘图,面向业务和科研需求开发。下图展示了 cedarkit-maps 项目的定位。从功能角度来看,cedarkit-maps 仅负责对气象数据的可视化,使用 cedarkit 工具栈中的其他工具库处理数据,例如数据准备工具库 reki 和数据计算工具库 cedarkit-comp。从系统角度来看,cedarkit-maps 仅提供基础的绘图工具库,具体的图片产品由更上一层的库实现,而将数据与图片相关联的业务系统绘图任务则由最上一层的库实现。

图 cedarkit-maps 绘图工具库的定位

cedarkit项目

cedarkit 项目是 cemc-oper 小组正在开展的一个项目,旨在提供基于 Python 的数值天气预报模式研发支撑工具集。该项目正在积极寻求 CEMC 内部的经费支持,并期望能为 CEMC 基于 Python 开发的各项工具提供底层功能支撑。cedarkit 项目中的大部分工具库均为开源项目,少部分为内部开源,欢迎 CEMC 内外同行一起参与到该项目的开发中。

cedarkit 是 CEMC 地球系统数据分析和可视化工具套件 (CEMC Earth Data Analysis and Rendering Toolkit,cedarkit) 的简称,提供一组 Python 工具集,用于分析地球系统数值预报数据。该套件基于 Python 科学计算社区的开源工具库开发,并开源发布。

cedarkit 项目参考 earthkit 和 geocat 系列工具集的设计理念,按照不同功能将整个工具套件拆分为不同的工具库,目前提供数据准备计算可视化等三个方面的工具库。不同工具库既可以组合使用,也可以单独使用,方便将工具库集成到其他 Python 项目中。多项目组合的形式也方便更多开发者参与项目开发工作,提高开发效率,降低维护成本。

cedarkit 采用分层架构,分为基础层、套件层和应用层。

  • 基础层由 Python 科学计算社区的开源工具库构成,包括
    • 数据结构库 numpy、pandas、xarray
    • 数据编解码库 eccodes
    • 科学计算库 scipy、dask
    • 可视化库 matplotlib、cartopy。
  • 套件层由 cedarkit 套件的多个工具库构成
    • 数据:计划项目 cedarkit-data (暂定),依赖以下项目
    • reki:数据准备
    • cmadaas-client (暂定 nuwe-cmadaas-python):CMADaaS MUSIC 接口客户端,数据源   - 计划项目 godas-client (暂定):GODAS 客户端 (CEMC 格点数据库),数据源,内部开源 (暂定)
    • 计算:正在开发的项目 cedarkit-comp
    • 可视化:正在开发的项目 cedarkit-maps
    • cedarkit-maps-examples:绘图示例,内部开源 (未来计划开源)
  • 应用层是基于 cedarkit 套件开发的应用库,主要应用于业务系统。

图 cedarkit 项目架构图,其中橙色框组件由 Python 科学计算社区开发,白色框是 CEMC 已开发的项目,黄色框是新开发的工具,虚线框是计划开发的工具

数据准备工具库将各类气象数据加载为 Python 常见的科学计算数据对象,形成标准数据结构。

  • 网格点数据使用 xarray.DataArray 表示
  • 表格数据使用 pandas.DataFrame 表示

计算工具库提供各类气象算法,为了方便其他工具库调用,以 numpy.ndarray 为数据结构,并提供适应标准数据结构的包装函数。

可视化工具库为数据准备工具库加载的标准数据结构提供可视化功能。

本文主要介绍 cedarkit 项目的可视化工具库 cedarkit-maps。

目标

cedarkit-maps 工具库的开发目标如下:

  • 快速绘制高质量气象平面
  • 开发面向数值预报模式研发与业务应用的通用绘图工具库,替换现有业务系统的 NCL、GrADS 绘图脚本,支持其他 Python 工具库的可视化绘图需求
  • 提供图片文件输出,并支持 Jupyter Notebook

指导思想

cedarkit-maps 工具库基于以下理念进行设计和开发:

  • 数据。使用统一格式表示数据。提供跨工具一致的数据表示方法。使用开源社区通用库表示数据结构,即数据准备工具库加载的标准数据结构。
  • 计算。数据计算与绘图分离。
  • 定制。绘图功能与图形样式分离。
  • 接口。提供各个组件的通用函数,也提供最终产品的集成函数。绘图脚本能够对接不同业务系统的数据,提高集约化程度。(通过基于 cedarkit-maps 开发的其它库实现)

cedarkit-maps 工具库将继承 CEMC 已有图片产品的绘图风格,吸收 CMA-HZ-CAEF 后处理绘图包的设计理念和功能实现,参考 earthkit、metdig 等开源项目的 API 接口设计,为 CEMC 开发的其他 Python 工具提供可视化功能支撑。

设计

总体设计

cedarkit-maps 封装 Python 科学计算社区的开源工具库的可视化绘图功能,提供高层和低层两类 API 接口,并内置了颜色表和地图资源。

下图是 cedarkit-maps 的架构设计图,采用分层架构设计。 第一层是项目依赖的开源工具库,包括可视化工具库 Matplotlib 和 Cartopy,数据结构库 Numpy、Pandas 和 Xarray。 第二层是低层 API 接口层,包括多个模块,graph 提供等值线、填充图、风场图等图形绘制函数,util 提供绘图框、标题、颜色条、网格线、坐标等底图组件函数,colormap 提供 NCL 色表加载函数,map 提供地图底图加载功能。colormap 和 map 模块使用了项目内置的资源文件 (resources),目前项目内置的资源包括 NCL 的部分色表和来自开源项目 china-shapefile 的中国地图底图 shapefile 文件。另外,map 模块还支持 CEMC 内部项目 cemc-meda-data 中 CEMC 中国地图底图资源。 第三层是高层 API 接口层,提供绘图板 (Panel)、图片 (Chart)、图层 (Layer) 等基础绘图类,并通过 domains 和 style 模块提供绘图定制功能。domains 定制底图布局,内置中国区域、世界区域、集合预报邮票图等多个 CEMC 图片布局。style 定制等值线、填充图、风场等图形的绘图样式。

cedarkit-maps 也具备一定的扩展能力,map 模块支持扩展新的地图底图,domains 模块支持定义新的底图布局,未来将进一步提高工具库的扩展性。

图 cedarkit-maps 架构图,包括 cedarkit-maps 和 cemc-meda-data 两个项目

使用 cedarkit-maps 工具库绘图的通用数据流如下图所示。

  1. 使用数据准备工具库 (reki) 从数据文件中 (GRIB2 file) 加载要素,得到数据要素场 (Data Field);
  2. (可选) 使用计算库 (cedarkit-comp) 对数据要素场进行处理,得到新的数据要素场;
  3. 创建底图布局 (Domain) 和绘图样式 (Style),使用绘图工具库 (cedarkit-maps) 的高层 API 接口对数据要素场进行绘图,得到图片 (plots)。

图 cedarkit-maps 数据流图。蓝色框表示工具库,白色框表示文件,黄色圆框表示内部数据结构

API设计

cedarkit-maps 提供多层 API 接口。

高层 API 提供简单易用的面向对象接口,可以使用内置的底图布局绘制图形。 低层 API 提供函数式接口,封装底层绘图库 Matplotlib 和 Cartopy 的函数,可以更精确控制图形组件。

高层API功能通过调用低层API实现。

高层API

参考 earthkit-maps 设计面向用户的高层 API,提供简单易用的面向对象式接口,隐藏内部绘图实现。

高层 API 主要由以下 3 种类构成:

  • Panel:绘图板,如下图所示
    • Chart:一个独立的图片,通常是一个地图,可能包含多个子图
    • Layer:图片中的一个子图图层,比如南海区域子图,在子图中绘制图形
  • XYDomain:底图布局
    • MapDomain:地图布局
  • Style:图形样式
    • ContourStyle:等值线/填充图样式
    • BarbStyle:风向杆图样式

下图展示了 PanelChartLayer 的层次关系。

图 cedarkit-maps 高层 API 接口绘图板示例,绘图板 Panel 包含多个图片 Chart,Chart 包含多个图层 Layer,图形在 Layer 中绘制

下面代码使用内置的东亚区域 (EastAsiaMapDomain) 绘制 2 米温度填充图(省略部分代码),展示高层 API 的一般用法。通常步骤如下:

  • 准备数据 t_2m_field
  • 创建绘图样式 t_2m_style
    • 定义色表 t_2m_color_map
    • 定义层次 t_2m_level
  • 创建布局 domain
  • 创建绘图板 panel
  • 绘制图形 panel.plot(...),多次调用绘制叠加图形
  • 定制绘图板
    • 设置标题 domain.set_title(...)
    • 添加颜色条 domain.add_colorbar(...)

代码 cedarkit-maps 高层 API 使用示例

from cedarkit.maps.style import ContourStyle  
from cedarkit.maps.chart import Panel  
from cedarkit.maps.domains import EastAsiaMapDomain

# 省略 t_2m_field,t_2m_color_map 和 t_2m_color_map 的定义代码

# 样式,等值线样式
t_2m_style = ContourStyle(  
    colors=t_2m_color_map,  
    levels=t_2m_level,  
    fill=True,  
)

# 布局,东亚区域
domain = EastAsiaMapDomain()

# 画板
panel = Panel(domain=domain)

# 绘图
panel.plot(t_2m_field, style=t_2m_style)  

# 定制画板
domain.set_title(  
    panel=panel,  
    graph_name="2m Temperature (C)",  
    system_name=system_name,  
    start_time=start_time,  
    forecast_time=forecast_time,  
)  
domain.add_colorbar(panel=panel, style=t_2m_style)  

低层API

低层 API 提供函数式接口,封装了 Matplotlib 和 Cartopy 的绘图功能,提供方便使用的绘图函数。

注:cedarkit-maps 正在开发过程中,未来低层 API 接口可能会有调整,建议使用高层 API 接口。

按照不同类型可以分为以下几类:

  • graph 模块:图形绘制函数
  • util 模块:底图组件函数
  • colormap 模块:颜色表函数
  • map 模块:地图底图模块
graph模块

graph 模块为每类图形提供绘图函数:

函数说明
add_contourf填充图
add_contour等值线
add_contour_label等值线/填充图标签
add_barb风向杆图
下面代码展示如何使用 add_contour 函数绘制海平面气压填充图,省略了部分代码。

代码 cedarkit-maps 低层 API 接口使用示例

import cartopy.crs as ccrs
import numpy as np
from cedarkit.maps.colormap import get_ncl_colormap
from cedarkit.maps.graph import add_contourf

# 省略:ax, prmsl

projection = ccrs.PlateCarree()
prmsl_contour_lev = np.arange(1000, 1035, 5)
prmsl_color_map = get_ncl_colormap(  
    "rainbow+white+gray",  
    index=np.array([-1, 80, 100, 120, 140, 170, 190, 205])  
)

cs = add_contourf(  
    ax,  
    field=prmsl,  
    projection=projection,  
    levels=prmsl_contour_lev,  
    cmap=prmsl_color_map,  
    zorder=0,  
)
util模块

util 模块提供底图组件相关函数,部分函数参数通过特定的类设置,例如 set_map_box_title 使用GraphTitle 中提供的标题文本和位置配置来添加标题。下表列出了 util 模块的主要函数

注:util 模块包含待弃用函数,未在表格中列出

表 cedarkit.maps.util 模块主要函数说明

类别函数/类说明备注
绘图框draw_map_box为图形添加矩形边框
标题GraphTitle标题配置,边框四角 + 主标题Class
fill_graph_title生成标题文本
set_map_box_title为图形边框设置标题
颜色条GraphColorbar颜色条配置Class
add_map_box_colorbar为图形添加颜色条
区域set_map_box_area设置绘图区域范围和比例
坐标轴set_map_box_axis设置坐标轴刻度和标签
draw_map_box_gridlines绘制坐标轴网格线
绘图实用函数clear_xarray_plot_components清除 Xarray 自动绘图生成的图片组件
clear_axes隐藏坐标轴、边框线、坐标短线、坐标标签
colormap模块

colormap 模块提供对 cedarkit-maps 内置的 NCL 颜色和颜色表的支持。

表 cedarkit.maps.colormap 模块主要函数说明

函数说明
get_ncl_colormap从内置的 NCL 颜色表返回 Matplotlib 色表对象
generate_colormap_using_ncl_colors根据 NCL 颜色名称列表返回 Matplotlib 色表对象
下面代码使用 generate_colormap_using_ncl_colors 函数将 NCL 颜色名列表转为 Matpltolib 色表对象 matplotlib.colors.ListedColormap
from cedarkit.maps.colormap import generate_colormap_using_ncl_colors

rain_1h_colormap = generate_colormap_using_ncl_colors([  
    "white",  
    "DarkOliveGreen3",  
    "ForestGreen",  
    "DeepSkyBlue",  
    "blue",  
    "magenta",  
    "DeepPink4",  
    "brown",
], "rain_1h_colormap")
map模块

map 模块为图片提供与中国边界相关的地图底图,默认使用 Cartopy 内置的底图,同时也支持 cedarkit-maps 内置的 china-shapfiles 项目底图和 cemc-meda-data 项目中的 CEMC 内部地图底图。

基础底图类 MapBase 包含一系列函数,返回特定的地图对象。通过继承该类可以支持新的地图底图。MapBase 的成员函数如下:

cedarkit.maps.map.MapBase 成员函数说明

成员函数说明
__init__创建函数,需要设置绘图种类 MapType
coastline海岸线
land陆地
rivers河流
lakes湖泊
china_coastline中国海岸线
china_borders中国国界
china_provinces中国省界
china_rivers中国河流
china_nine_lines九段线
global_borders全球国界

默认使用内置的 china-shapefile 项目底图,可以使用 set_default_map_package 修改默认底图。例如下面代码使用 CEMC 内部地图底图。

注:需要单独安装 CEMC 内部软件包 cemc-meda-data。

from cedarkit.maps.map import set_default_map_package

set_default_map_package("cedarkit.maps.map.cemc")

下图展示了 MapBase 各个主要函数对应的地图组件,底图来自 cemc-meda-data 项目。

图 cedarkit.maps.map.MapBase 成员函数示例,底图来自 cemc-meda-data 项目

示例

CMA-GFS

CMA-MESO

CMA-REPS

实现

本节介绍 cedarkit-maps 核心部分的具体实现。

数据

使用单独的数据加载工具库将数据加载为标准数据结构。 对于格点数据,将其加载为 xarray.DataArray 对象,例如 reki 库支持 GRIB、GrADS 和 NetCDF 等格点数据格式。一个典型的格点数据结构包含以下坐标信息:

  • longitude:经度
  • latitude:纬度
  • 层次,不同类型层次使用不同名称
    • pl:等压面,hPa
    • ml:模式面
    • eccodes 支持的层次,例如:
      • heightAboveGround:距地面高度,米
  • step:预报时效
  • valid_time:预报时间
  • start_time:起报时间
  • member:集合成员编号,可选,仅集合预报数据

当加载单个要素场时,reki 返回以 latitudelongitude 为坐标的二维数组。 批量加载多层次要素场时,返回以层次、latitudelongitude 为坐标的三位数组。

用户可以根据自己的需求按照特定的维度合并多个要素场。比如按照 step 坐标合并得到时间序列数据。

下面代码从 CMA-MESO 的 GRIB2 产品中加载 2 米温度要素场,并将单位转换为摄氏度。

from reki.data_finder import find_local_file
from reki.format.grib import load_feild_from_file

file_path = find_local_file(  
    "cma_meso_3km/grib2/orig",  
    start_time=start_time,  
    forecast_time=forecast_time,  
    storage_base="/mnt/m/archive"  
)  
  
t_2m_field = load_field_from_file(  
    file_path,  
    parameter=t_info.eccodes_short_name,  
) - 273.15
print(t_2m_field)

输出 2 米温度的概要信息如下:

<xarray.DataArray '2t' (latitude: 1671, longitude: 2501)> Size: 33MB
array([[-10.73962891, -10.77962891, -10.80962891, ..., -26.80962891,
        -27.34962891, -27.47962891],
       [ -3.86962891,  -3.90962891,  -2.67962891, ..., -21.41962891,
        -23.03962891, -27.60962891],
       [ -3.71962891,  -3.80962891,  -2.66962891, ..., -20.38962891,
        -22.30962891, -28.05962891],
       ...,
       [ 26.84037109,  26.86037109,  26.84037109, ...,  25.87037109,
         25.88037109,  26.13037109],
       [ 26.84037109,  26.85037109,  26.86037109, ...,  25.87037109,
         25.88037109,  26.16037109],
       [ 26.89037109,  26.87037109,  26.87037109, ...,  25.88037109,
         25.90037109,  26.16037109]])
Coordinates:
    time               datetime64[ns] 8B 2024-02-05T12:00:00
    step               timedelta64[ns] 8B 1 days
    valid_time         datetime64[ns] 8B 2024-02-06T12:00:00
    heightAboveGround  int64 8B 2
  * latitude           (latitude) float64 13kB 60.1 60.07 60.04 ... 10.03 10.0
  * longitude          (longitude) float64 20kB 70.0 70.03 70.06 ... 145.0 145.0

下面代码加载 CMA-REPS 所有成员 048 时效 24 小时降水并返回 3 维数组 (numberlatitudelongitude) :

代码 加载 CMA-REPS 所有成员的 24 小时降水要素

import pandas as pd
import xarray as xr
from reki.data_finder import find_local_file
from reki.format.grib.eccodes import load_field_from_file

forecast_time = pd.to_timedelta("48h")
interval = pd.to_timedelta("24h")
previous_forecast_time = forecast_time - interval

fields = []
for number in np.arange(0, 15):
    logger.info(f"loading apcp for {number}...")
    file_path = find_local_file(
        data_type,
        start_time=start_time,
        forecast_time=forecast_time,
        number=number,
    )
    apcp_field = load_field_from_file(
        file_path,  
        parameter=parameter,
    )
    file_path_previous = find_local_file(
        data_type,
        start_time=start_time,
        forecast_time=previous_forecast_time,
        number=number,
    )
    apcp_field_previous = load_field_from_file(
        file_path_previous,
        parameter=parameter
    )  
    apcp_24h_field = apcp_field - apcp_field_previous
    fields.append(apcp_24h_field)
  
apcp_24h_field = xr.concat(fields, dim="number")

当然,可以使用其他工具库加载格点数据,只要数据结构与上述保持一致即可,即二维数组的维度是 (latitudelongitude)。

计算

使用单独的气象计算工具库完成数据计算。

例如正在开发的 cedarkit-comp 工具库包含了 NCL 的 smth9 函数,实现 9 点平滑。cedarkit-comp 工具库提供函数式接口,并使用基础数据结构 numpy.ndarray 以方便算法开发和使用。如果需要处理 cedarkit 工具栈使用的 xarary.DataArray 对象,需要使用辅助函数 apply_to_xarray_values。下面代码从 CMA-MESO 数据中加载雷达组合反射率 (CR),并执行 2 次 smth9 平滑。

代码 cedarkit-comp 库调用示例

from reki.format.grib.eccodes import load_field_from_file
from cedarkit.comp.smooth import smth9  
from cedarkit.comp.util import apply_to_xarray_values

cr_field = load_field_from_file(  
    file_path,  
    parameter=cr_info.wgrib2_name,  
)

cr_field = apply_to_xarray_values(cr_field, lambda x: smth9(x, 0.5, -0.25, False))
cr_field = apply_to_xarray_values(cr_field, lambda x: smth9(x, 0.5, -0.25, False))

绘图组件

底图布局

底图布局 (Domain) 提供预定义的图片布局,自动绘制地图、坐标轴、标题等绘图组件,作为高层 API 接口 Panel 的参数,方便用户快速创建图像。

下图以 cedakit-maps 内置的东亚/中国区域 (EastAsiaMapDomain) 说明底图布局预先定义的常用绘图组件。该布局包括两个图层,中国区域主图层 (main layer) 和 南海区域子图层 (sub layer),每个图层都绘制地图 (map),并在右下角标注了地图信息 (map info)。两个图层都绘制了经纬度网格线 (axis gridline),主图层还标注了经纬度标签 (axis tick/label)。主图层被绘图框 (map box) 包围。仅需将 EastAsiaMapDomain 对象作为参数传递给 Panel 对象,上述绘图组件均默认自动生成。同时,EastAsiaMapDomain 还提供了一些定制函数,用于添加附加绘图组件,包括绘图框四角的标题 (title)和绘图框右侧的颜色条 (colorbar)。

图 cedarkit-maps 内置底图布局 EastAsiaMapDomain 示意图

基础布局类

当前只有一个基础布局类,即二维平面图布局抽象类 XYDomain。该抽象类包含一个成员函数 render_panel ,完成布局的绘制任务。Panel__init__ 函数需要 XYDomain 作为参数,并会调用 render_panel 函数绘制底图布局,所以对于 XYDomain 的定制需要在传递给 Panel 前完成,Panel 创建后对 XYDomain 配置的修改将不会影响当前 Panel 的底图布局。

下图展示 cedarkit-maps 底图布局在高层 API 中的使用流程。首先创建布局对象 XYDomain,可以对布局对象进行定制,将对象传递给 Panel.__init__ 函数用于创建绘图板对象 Panel__init__ 函数会调用 XYDomain.render_panel 函数绘制地图布局,包括创建子图 Chart 和图层 Layer,创建绘图区域 (Axes),绘制坐标轴 (Axis)、文本 (text),地图 (Map) 和绘图框 (Map box)。

图 cedarkit-maps 底图布局在高层 API 中的使用流程示例

内置布局

内置布局均继承自 XYDomain。子类 MapDomain 是带地图布局的基类。TimeStempAndLevelXYDomain 表示时序高度剖面图布局。

图 cedarkit-maps 内置底图布局类的关系图,所有类均继承自 XYDomain

当前内置的底图布局包括以下几种。

表 cedarkit-maps内置的底图布局

名称说明
东亚EastAsiaMapDomain东亚,或中国范围,带南海子图
中国区域CnAreaMapDomain中国区域,例如华北、华南等
北半球NorthPolarMapDomain北半球范围
欧亚大陆EuropeAsiaMapDomain欧亚大陆范围
全球GlobalMapDomain全球范围
全球区域GlobalAreaMapDomain全球任意区域
中国区域集合预报邮票图EnsCnMapDomain中国区域,集合预报邮票图

东亚/中国

东亚区域 (cedarkit.maps.domains.EastAsiaMapDomain) 用于绘制中国范围图形,包括中国区域和南海子图两个图层。默认区域范围北纬 15 - 55,东经 70 - 140,使用简易圆柱投影 (cartopy.crs.PlateCarree)。

图 cedarkit-maps 内置底图布局,东亚区域 (EastAsiaMapDomain)

import numpy as np
import pandas as pd
import matplotlib.colors as mcolors
  
from cedarkit.maps.chart import Panel
from cedarkit.maps.domains import EastAsiaMapDomain
from cedarkit.maps.map import set_default_map_package
from cedarkit.maps.style import ContourStyle

set_default_map_package("cedarkit.maps.map.cemc")
  
domain = EastAsiaMapDomain()
panel = Panel(domain=domain)

map_colors = np.array([
    (255, 255, 255),
    (0, 0, 0),
    (255, 255, 255),
    (0, 200, 200),
    (0, 210, 140),
    (0, 220, 0),
    (160, 230, 50),
    (230, 220, 50),
    (230, 175, 45),
    (240, 130, 40),
    (250, 60, 60),
    (240, 0, 130),
    (0, 0, 255),
    (255, 140, 0),
    (238, 18, 137)
], dtype=float) / 255
colormap = mcolors.ListedColormap(map_colors)
  
wind_speed_colormap = mcolors.ListedColormap(colormap(np.array([2, 3, 4, 5, 6, 7, 8, 9, 10, 11])))
wind_speed_contour_lev = np.array([3.4, 5.5, 8, 10.8, 13.9, 17.2, 20.8, 24.5, 28.5])
wind_speed_style = ContourStyle(
    colors=wind_speed_colormap,
    levels=wind_speed_contour_lev,
    fill=True,
)  

domain.set_title(
    panel=panel,
    graph_name="10m Wind, 10m Wind Speed(m/s, shadow)",
    system_name="CMA-MESO",
    start_time=pd.to_datetime("2024-03-21 00:00"),
    forecast_time=pd.to_timedelta("24h"),
)
domain.add_colorbar(panel=panel, style=wind_speed_style)

panel.show()

中国区域

中国区域 (cedarkit.maps.domains.CnAreaMapDomain) 继承自东亚区域 (cedarkit.maps.domains.EastAsiaMapDomain),默认不绘制南海子图。

图 cedarkit-maps 内置底图布局,中国区域 (CnAreaDomain)

import pandas as pd  
import matplotlib.colors as mcolors  
  
from cedarkit.maps.chart import Panel  
from cedarkit.maps.domains import CnAreaMapDomain  
from cedarkit.maps.map import set_default_map_package  
from cedarkit.maps.style import ContourStyle  
  
from cedarkit_maps_examples.systems.meso.area.plot_area import get_plot_area  
from cedarkit.maps.colormap import get_ncl_colormap

set_default_map_package("cedarkit.maps.map.cemc")

plot_area = get_plot_area("NorthChina")

domain = CnAreaMapDomain(area=plot_area.area)
panel = Panel(domain=domain)

color_map = get_ncl_colormap("BlAqGrYeOrReVi200")
t_2m_level = [-24, -20, -16, -12, -8, -4, 0, 4, 8, 12, 16, 20, 24, 28, 32]
color_index = [2, 12, 22, 32, 42, 52, 62, 72, 82, 92, 102, 112, 122, 132, 142, 152]
			   
t_2m_color_map = mcolors.ListedColormap(color_map(color_index))
t_2m_style = ContourStyle(
    colors=t_2m_color_map,
    levels=t_2m_level,
    fill=True,
)
domain.add_colorbar(panel=panel, style=t_2m_style)

domain.set_title(
    panel=panel,
    graph_name="2m Temperature (C)",
    system_name="CMA-MESO",
    start_time=pd.to_datetime("2024-03-20 00:00"),
    forecast_time=pd.to_timedelta("24h"),
)

panel.show()

北半球

北半球底图布局 (cedarkit.maps.domains.NorthPolarMapDomain) 绘制北半球图形,使用北半球立体投影 (cartopy.crs.NorthPolarStereo),默认中心经度东经 110 度。

图 cedarkit-maps 内置底图布局,北半球 (NorthPolarMapDomain)

import numpy as np  
import pandas as pd  
  
from cedarkit.maps.chart import Panel  
from cedarkit.maps.domains import NorthPolarMapDomain  
from cedarkit.maps.map import set_default_map_package  
from cedarkit.maps.style import ContourStyle  
from cedarkit.maps.colormap import get_ncl_colormap

set_default_map_package("cedarkit.maps.map.cemc")  
domain = NorthPolarMapDomain()  
panel = Panel(domain=domain)  
  
contour_lev = np.arange(-16, 25, 2)  
color_map = get_ncl_colormap(  
    "rainbow+white+gray",  
    index=np.arange(38, 207, 8),  
)  
t_style = ContourStyle(  
    colors=color_map,  
    levels=contour_lev,  
    fill=True,  
)  
  
domain.set_title(  
    panel=panel,  
    graph_name="850 hPa TEMP ($^\circ$C) shadow",  
    system_name="CMA-GFS",  
    start_time=pd.to_datetime("2024-03-21 00:00"),  
    forecast_time=pd.to_timedelta("24h"),  
)  
domain.add_colorbar(panel=panel, style=t_style)  
  
panel.show()

欧亚大陆

欧亚大陆底图布局 (cedarkit.maps.domains.EuropeAsiaMapDomain) 绘制欧亚大陆区域图形,使用兰勃特等角圆锥投影 (cartopy.crs.LambertConformal),中心经度东经 95 度。

图 cedarkit-maps 内置底图样式,欧亚大陆区域 (EuropeAsiaMapDomain)

import numpy as np
import pandas as pd

from cedarkit.maps.chart import Panel
from cedarkit.maps.domains import EuropeAsiaMapDomain
from cedarkit.maps.map import set_default_map_package
from cedarkit.maps.colormap import get_ncl_colormap
from cedarkit.maps.style import ContourStyle

set_default_map_package("cedarkit.maps.map.cemc")
domain = EuropeAsiaMapDomain()
panel = Panel(domain=domain)

prmsl_color_map = get_ncl_colormap(
    "rainbow+white+gray",
    index=np.array([-1, 80, 100, 120, 140, 160, 180, 200, 205])
)
prmsl_contour_lev = np.arange(1000, 1036, 5)
prmsl_style = ContourStyle(
    colors=prmsl_color_map,
    levels=prmsl_contour_lev,
    fill=True,
)

domain.set_title(
    panel=panel,
    graph_name="500 hPa HGT (10gpm) line + MSLP (hPa) shadow",
    system_name="CMA-GFS",
    start_time=pd.to_datetime("2024-03-21 00:00"),
    forecast_time=pd.to_timedelta("24h"),
)
domain.add_colorbar(panel=panel, style=prmsl_style)

panel.show()

全球

全球底图布局 (cedarkit.maps.domains.GlobalMapDomain) 绘制全球范围图形,使用简易圆柱投影 (cartopy.crs.PlateCarree),默认中心经度东经 80 度。默认使用 Cartopy 的底图,并填充陆地颜色。

图 cedarkit-maps 内置底图布局,全球 (GlobalMapDomain)

import pandas as pd
import numpy as np

from cedarkit.maps.chart import Panel
from cedarkit.maps.domains import GlobalMapDomain
from cedarkit.maps.map import set_default_map_package
from cedarkit.maps.style import ContourStyle, ColorbarStyle
from cedarkit.maps.colormap import get_ncl_colormap

set_default_map_package("cedarkit.maps.map.cemc")
domain = GlobalMapDomain()
panel = Panel(domain=domain)

t_contour_lev = np.arange(-80, 60 + 4, 4)
colormap = get_ncl_colormap("temp_19lev", count=len(t_contour_lev) + 1, spread_start=0, spread_end=17)
t_style = ContourStyle(
    colors=colormap,
    levels=t_contour_lev,
    fill=True,
    colorbar_style=ColorbarStyle(
        label_levels=t_contour_lev[::2]
    )
)

domain.set_title(
    panel=panel,
    graph_name="850 hPa HGT (10gpm)",
    system_name="CMA-GFS",
    start_time=pd.to_datetime("2024-02-28 00:00"),
    forecast_time=pd.Timedelta(hours=24),
)
domain.add_colorbar(panel, style=t_style)

panel.show()

全球区域

全球区域底图布局 (cedarkit.maps.domains.GlobalAreaMapDomain) 绘制全球范围任意区域的图形,继承自全球底图布局 (GlobalMapDomain),并增加了国境线。

图 cedarkit-maps 内置底图布局,全球区域 (GlobalAreaMapDomain)

import numpy as np
import pandas as pd

from cedarkit.maps.chart import Panel
from cedarkit.maps.domains import GlobalAreaMapDomain
from cedarkit.maps.map import set_default_map_package
from cedarkit.maps.style import ContourStyle
from cedarkit.maps.colormap import get_ncl_colormap

set_default_map_package("cedarkit.maps.map.cemc")
area = [45, 110, -10, 45]
domain = GlobalAreaMapDomain(area=area)
panel = Panel(domain=domain)
  
t_contour_lev = np.array([0, 4, 8, 12, 16, 20, 24, 28, 32, 36])
colormap = get_ncl_colormap("WhBlGrYeRe", index=np.array([2, 8, 14, 20, 45, 55, 65, 72, 80, 90, 100]))
t_style = ContourStyle(
    colors=colormap,
    levels=t_contour_lev,
    fill=True,
)

domain.set_title(
    panel=panel,
    graph_name=f"2m Temperature ($^\circ$C,shadow)",
    system_name="CMA-GFS",
    start_time=pd.to_datetime("2024-02-28 00:00"),
    forecast_time=pd.Timedelta(hours=24),
)
domain.add_colorbar(panel, style=t_style)

panel.show()

中国区域集合邮票图

中国区域集合邮票图底图布局 (cedarkit.maps.domains.EnsCnMapDomain) 用于绘制集合预报邮票图,使用简易圆柱投影 (cartopy.crs.PlateCarree)。

注:集合预报邮票图底图布局正在开发中,未来可能会有变化。

图 cedarkit-maps 内置底图布局,中国区域集合预报邮票图 (EnsCnMapDomain)

import pandas as pd
import matplotlib.colors as mcolors

from cedarkit.maps.chart import Panel, Schema
from cedarkit.maps.domains.ens_cn import EnsCNMapDomain
from cedarkit.maps.style import ContourStyle

color_list = ['#fefefe', '#c9fbc2', '#02db01', '#03c5c4', '#203dfc', '#f40984', '#a9a9a9']
color_map = mcolors.ListedColormap(color_list)
contour_style = ContourStyle(
    levels=[0.1, 10, 25, 50, 100, 250],
    colors=color_map,
    fill=True,
)

domain = EnsCNMapDomain(enable_max=True)
panel = Panel(domain=domain, schema=Schema(figsize=(10, 8)))
domain.set_title(
    panel=panel,
    graph_name="rain24",
    system_name="CMA-REPS",
    start_time=pd.to_datetime("2023-01-01 00:00"),
    forecast_time=pd.to_timedelta("48h")
)
domain.add_colorbar(panel=panel, style=contour_style)

panel.show()

地图

cedarkit-maps 的 cedarkit.maps.map 模块负责地图加载,主要解决 Cartopy 的中国地图问题。目前支持两种地图底图,通过类的形式封装,所有地图包类继承自 MapBase

  • 默认地图 (cedarkit.maps.default.DefaultMap):使用内置的开源项目 china-shapefiles 资源文件绘制中国地图。
  • CEMC 地图 (cedarkit.maps.map.cemc.CemcMap):使用内部项目 cemc-meda-data 中的 CEMC 地图绘制地图。

图 cedarkit-maps 内置地图包

cedarkit-maps 默认使用 DefaultMap 底图包,可以使用 set_default_map_package 修改默认使用的地图包,底图布局类调用 get_map_class 函数获取当前默认的地图包类。

下面代码将默认地图包设为项目内置的 CEMC 地图包。

注:需要单独安装内部项目库 cemc-meda-data 才能使用 CEMC 地图包

from cedarkit.maps.map import set_default_map_package

set_default_map_package("cedarkit.maps.map.cemc")

增加新的地图包需要创建一个新的模块,继承 MapBase 类并重载各个成员函数,在模块中将子类名称赋值给新模块的 map_class 变量。通过类似 module_a.new_map_module 的模块导入路径使用新增的地图包。

图 cedakrit-maps 新增地图包示例,新地图包 new_map_module 包含继承自 MapBase 的子类 NewMap,包中 map_class 变量赋值为 NewMap。使用 set_default_map_package 函数设置新地图包后,Domain 对象可以正确创建 NewMap 对象。

内置地图

默认地图

默认地图 (cedarkit.maps.map.default.DefaultMap) 使用 Cartopy 内置的 NaturalEarthFeature 绘制海岸线、陆地、河流和湖泊,使用开源项目 china-shapefiles 中国地图绘制中国边界、省界和南海九段线。china-shapefiles 项目已内置在 cedarkit-maps 的资源目录 resources 中,无需单独安装。

图 cedarkit-maps 内置地图包,默认底图 (DefaultMap),使用 china-shapefiles 绘制中国地图

CEMC地图

CEMC 地图 (cedarkit.maps.map.cemc.CemcMap) 继承自默认底图,使用内部项目 cemc-meda-data 中的 CEMC 底图绘制中国海岸线、中国边界、中国省界、中国河流、南海九段线,并在 Global 模式中绘制全球国界。

图 cedarkit-maps 内置地图包,CEMC 地图包 (CemcMap),使用 cemc-meda-data 绘制中国地图

颜色

cedarkit-maps 使用 Matplotlib 的颜色功能,也支持 NCL 的颜色,并内置了部分 NCL 颜色表。

get_ncl_colormap 函数加载内置的 NCL 颜色表,内置的色表保存在 resources/colormap/ncl 目录中。下面代码加载 precip_diff_12lev 为雨雪生成颜色表。

from cedarkit.maps.colormap import get_ncl_colormap

rain_snow_color_map = get_ncl_colormap(
	"precip_diff_12lev", 
	index=np.array([6, 5, 4, 3, 2, 1])
)

generate_colormap_using_ncl_colors 函数将 NCL 颜色名转为 Matplotlib 颜色表,NCL 颜色保存在 resources/colormap/ncl/ncl_colors.csv 文件中。下面代码使用 NCL 颜色名为降水生成颜色表。

from cedarkit.maps.colormap import generate_colormap_using_ncl_colors

rain_color_map = generate_colormap_using_ncl_colors(  
    [  
        "transparent",  
        "PaleGreen2",  
        "ForestGreen",  
        "DeepSkyBlue",  
        "blue1",  
        "magenta1",  
        "DeepPink3",  
        "DarkOrchid4"  
    ],  
    name="rain"  
)

颜色条的位置与底图布局密切相关,因此向图形中添加颜色条的函数通常放在底图布局类中。例如下面代码使用东亚底图布局,调用 add_colorbar 在绘图板中为海平面气压填充图添加颜色条。

from cedarkit.maps.domains import EastAsiaMapDomain

# ...skip...
domain = EastAsiaMapDomain()
panel = Panel(domain=domain)
panel.plot(mslp_field, style=mslp_style)
# ...skip...
domain.add_colorbar(panel=panel, style=mslp_style)

文本

cedarkit-maps 中的文本主要是如下图所示的主标题 (main_label) 和绘图框四角标题 (top_left_label, top_right_label, bottom_left_label 和 bottom_right_label),图中还列出了定位上述标题需要的五个位置参数 (top, bottom, left, right 和 main_post)。以上信息均包含在图形标题设置类 GraphTitle 中。

图 cedarkit-maps 标题设置类 GraphTitle 示意图

与颜色条类似,标题的位置与底图布局密切相关,因此向图形中添加标题的函数通常放在底图布局类中。例如下面代码使用东亚底图布局,调用 set_title 在绘图板中为海平面气压填充图添加四角标题。

import pandas as pd
from cedarkit.maps.chart import Panel
from cedarkit.maps.domains import EastAsiaMapDomain

domain = EastAsiaMapDomain()
panel = Panel(domain=domain)
# ...skip...
domain.set_title(
    panel=panel,
    graph_name="500 hPa HGT (10gpm) line + MSLP (hPa) shadow",
    system_name="CMA-GFS",
    start_time=pd.to_datetime("2024-03-20 00:00"),
    forecast_time=pd.to_timedelta("24h"),
)

图形与样式

为标准数据结构封装 Matplotlib 和 Xarray 的绘图功能,提供多种图形的绘图函数。 高层 API 使用样式类 (Style) 作为参数。Style 类定义绘制图形需要的定制参数,部分参数名称与 Matplotlib 绘图函数的参数相同,每类图形有各自的样式类。 低层 API 则直接支持 Matplotlib 或 Xarray 绘图函数的参数。

下表给出当前已实现图形的高层 API、样式、低层 API 及调用的底层函数。

图形高层API样式低层API底层函数
等值线Layer.contourfContourStyleadd_contourfxarryDataArray.plot.contourf
填充图Layer.contourContourStyleadd_contourxarryDataArray.plot.contour
- 标签Layer.contour_labelContourLabelStyleadd_contour_labelmatplotlibAxes.clabel
- 颜色条Domain.add_colorbarColorbarStyleadd_map_box_colorbarmatplotlibFigure.colorbar
风场图Layer.barbBarbStyleadd_barbmatplotlibAxes.barbs

调用接口

本节以 CMA-GFS 世界气象中心网站全球 850 hPa 温度图片为例说明绘图样式在不同 API 层次中如何传递。

下面是温度图片示例和绘制代码,省略了变量设置和数据加载部分的代码。

图 cedarkit-maps 绘图示例,CMA-GFS 世界气象中心网站全球 850 hPa 温度场

import numpy as np  
import pandas as pd  
import xarray as xr  
  
from cedarkit.maps.style import ContourStyle, ColorbarStyle  
from cedarkit.maps.colormap import get_ncl_colormap  
from cedarkit.maps.chart import Panel  
from cedarkit.maps.domains import GlobalMapDomain  

# ...skip...

t_contour_lev = np.arange(-80, 60 + 4, 4)  
colormap = get_ncl_colormap("temp_19lev", count=len(t_contour_lev) + 1, spread_start=0, spread_end=17)  
t_style = ContourStyle(  
    colors=colormap,  
    levels=t_contour_lev,  
    fill=True,  
    colorbar_style=ColorbarStyle(  
        label_levels=t_contour_lev[::2]  
    ),  
    label=True,  
    label_style=ContourLabelStyle(  
        levels=t_contour_lev[::2],  
        colors="black",  
        background_color="white",  
        fontsize=6,  
    )  
)  
  
# plot  
domain = GlobalMapDomain()  
panel = Panel(domain=domain)  
panel.plot(t_field, style=t_style)  
  
domain.set_title(  
    panel=panel,  
    graph_name=f"{level} hPa TEMP ($^\circ$C)",  
    system_name=system_name,  
    start_time=start_time,  
    forecast_time=forecast_time,  
)  
domain.add_colorbar(panel, style=t_style)

图 cedarkit-maps 绘图样式在多层 API 接口中的调用链路,包括高层 API、低层 API 和社区开源库 API

高层 API

Panel.plot 方法用于在绘图板中绘制图形,接收绘图数据要素场 (t_field) 和绘图样式 (style)两个参数。

注意:本例中绘图数据是单个要素场,在 Panel.plot 内部会被包装成具有单个元素的列表对象。

panel.plot(t_field, style=t_style)

t_styleContourStyle 对象,定义了填充图的颜色表 (colors)、层次 (levels) 和填充开关 (fill)。对象中的 colorbar_styleColorbarStyle 对象,定义颜色条需要标注的层次 (label_levels)。对象中的 label_styleContourLabelStyle 对象,定义等值线/填充图标注的层次 (levels)、颜色 (colorsbackground_color) 和文本大小 (fontsize)。

Panel.plot 会调用绘图板中每个子图的绘图函数,即 Chart.plot。每个子图使用数据列表对象中的一个元素,并直接使用绘图样式对象。

graph = self.charts[i].plot(data=d, style=style, layer=layer)

Chart.plot 会在子图中的每个图层进行绘图,根据 style 参数的类型选择合适的绘图函数,例如本例中绘图样式是 ContourStyle 对象,并设置了 fill 属性,所以选择填充图绘制函数 layer.contourf ,并传递数据和绘图样式参数。

if isinstance(style, ContourStyle):
    if style.fill:
        result = layer.contourf(data=data, style=style)
    else:
        result = layer.contour(data=data, style=style)

Layer.contourf 函数调用低层 API 接口 cedarkit.maps.graph.add_contourf 函数实现填充图绘制,从 Layer 对象中获取绘图坐标轴 (ax) 和投影 (projection) 属性,并将绘图样式转换为函数参数 (levels, cmap)。

contour = add_contourf(  
    self.ax,
    field=data,
    levels=style.levels,
    projection=self.projection,
    cmap=style.colors,
    **kwargs
)

低层 API

add_contourf 函数调用开源项目 xarray 的绘图接口 xarray.DataArray.plot.contouf 实现填充图绘制,从数据对象 field 中计算得到取值范围 (vmin, vmax),传递样式相关参数 (levels, cmap(包含在 kwargs 中)),并设置了额外的默认参数 (add_colorbar, add_labels, extend)。

c = field.plot.contourf(
    ax=ax,
    transform=projection,
    vmin=min_level,
    vmax=max_level,
    levels=levels,
    add_colorbar=False,
    add_labels=False,
    extend="both",
    ylim=ylim,
    **kwargs
)

样式简介

等值线/填充图

ContourStyle 类定义等值线/填充图样式。

属性说明默认值
colors颜色None
levels层次None
linewidths线宽,用于等值线图None
linestyles线型,用于等值线图None
fill是否填充False
label是否添加标签False
label_style标签样式,见 ContourLabelStyleNone
colorbar_style颜色条样式,见 ColorbarStyleNone

等值线标签

ContourLabelStyle 类定义等值线/填充图标签样式。

属性说明默认值
fontsize字号None
inlineTrue
inline_spacing5
fmtNone
colors颜色None
manualFalse
zorderNone

风场图

BarbStyle 类定义风场图样式。

属性说明默认值
length4
linewidth0.5
pivot“middle”
barbcolor“red”
flagcolor“red”
barb_incrementsNone
colorbar_style颜色表样式,见 ColorbarStyleNone

颜色表

ColorbarStyle 类定义颜色表样式。

属性说明默认值
loc色表位置None
label名称None
label_levels标注层次None

总结

cedarkit-maps 项目提供气象数据二维绘图功能,基于 Python 科学计算社区的开源库开发,吸收国内外同行及 CEMC 绘图工具包的优秀设计理念,为科研和业务提供统一的气象绘图底层工具库。cedarkit-maps 提供面向对象的高层 API 接口,也提供函数式的低层 API 接口。目前已支持绘制等值线、填充图、风场图等多种图形,内置中国、北半球、欧亚大陆、世界等多个底图布局,已能实现绘制 CEMC 部分业务风格的图片,并通过 cemc-meda-data 项目提供 CEMC 地图底图支持。

cedarkit-maps 项目是 cedarkit 工具套件的成员,与 cedarkit 工具套件中的数据准备工具库 reki 和气象算法工具库 cedarkit-comp 一起提供从数据到可视化的完整工具链。

当前 cedarkit 项目仍处于早期开发阶段,后续将进一步完善现有功能,提供更强的定制和扩展能力,并继续开发新的特性。

下一步计划

开发阶段

cedarkit-maps 工具库的最终目标是为科研和业务提供气象数据的通用可视化绘图功能。开发阶段如下:

阶段名称说明
STAGE-1业务产品级固定样式图形,业务系统绘图 (NMC、WMC-BJ网站)
STAGE-2绘图组件级基础绘图功能,作为底层组件应用到 CEMC 的其他 Python 工具库中
STAGE-3数据计算级研发绘图,提供自动绘图功能
STAGE-4通用工具级从数据、分析到可视化的完整工具链 (cedarkit 工具套件)

当前正处于第一个开发阶段,开发业务绘图的脚本。

  • CMA-MESO 1KM 绘图脚本,待实现(2024)
  • CMA-REPS 3KM 绘图脚本,已有,考虑与 cedarkit-maps 整合(2024)
  • CMA-GFS 绘图脚本,待实现(2025)
  • CMA-TYM 绘图脚本,待实现(2025)

计划 2024 年借助 MESO 1KM 业务化的机会完成前两项任务,2025 年继续完成后两项任务,随后进入到第二个开发阶段。后续会根据工作任务需求适时调整开发计划。

开发组织方式

cedarkit 当前仅是个人向的探索项目,尚未成为 CEMC 的工作项目,后续将继续推动 Python 底层工具软件库进入到 CEMC 的正式工作项目中,争取尽早解决 cedarkit 项目的定位问题。

cedarkit-maps 工具库和 cedarkit 工具套件将持续探索团队合作的开发方式,通过多种途径寻求项目支持。

附录

本文所涉及的项目库:

cedarkit-maps:气象绘图库,开源

cedarkit-comp:气象算法库,开源

reki:数据准备工具库,开源

cemc-meda-data:CEMC地图底图资源库,CEMC 内部开源(需要CEMC组权限)

cedarkit-maps-examples:cedarkit-maps 项目绘图示例库,CMA 内部开源