炸鸡人博客 基本上无害
Posts with the tag matplotlib:

Cartopy 系列:画中国地图的工具箱 frykit

contourf

前言

最早笔者用 Python 画中国地图时,会准备 bou2_4p.shp 文件,然后封装一个读取 shapefile 并添加到 GeoAxes 上的函数,别的项目要用时就把数据和函数复制粘贴过去。Cartopy 系列:从入门到放弃 里就是这么做的。

后来工作中用到了 Clarmy 开发的 cnmaps 包,只用两行就能快速绘制地图,非常方便。同时萌生了自己实现一个功能类似的包的想法,遂开发出了 frykit

搞颜色系列:单色光光谱

前言

人眼可见色域在色度图中表现为彩色的马蹄形,单色光(monochromatic light)的颜色对应于马蹄的弧形边界。本文想将单色光的颜色按波长线性增大的顺序一字排开,用类似彩虹渐变图的形式展示单色光光谱。用 Python 的 Matplotlib 包来实现的话,很快就能决定画图思路:

  1. 读取 XYZ 颜色匹配函数(CMF)作为 XYZ 三刺激值。
  2. XYZ 变换为 sRGB,接着做 gamma 校正。
  3. 用 RGB 数组构造 ListedColormap 对象,用 plt.colorbar 画出。

RGB 要求范围在 $[0, 1]$,但 CMF 直接计算出的 RGB 既有负数分量,也有大于 1 的分量,所以必须采用一种方法处理范围外的分量。最后的画图效果会因处理方法的不同产生很大差别,例如下图的三条光谱:

three_colorbars.png

就采取了不同的处理方式,因此在发色、颜色过渡,和亮度表现上都大有不同。本文将尝试实现不同的效果并加以分析。完整代码和相关数据见 我的 Github 仓库

搞颜色系列:绘制 CIE 1931 色度图

前言

1920 年代末 Wright 和 Guild 的颜色匹配实验发展出了用红绿蓝三基色(primaries)定量表示所有人眼可见颜色的 CIE RGB 色彩空间,1931 年国际照明委员会(CIE)通过对 CIE RGB 色彩空间做线性变换得到了 CIE XYZ 色彩空间。XYZ 空间里的人眼可见色域(gamut of human vision)是一块从原点出发,向无限远处不断延伸的立体区域。将这块区域投影到 $X + Y + Z = 1$ 的平面上,就能画出方便展示的 CIE 1931 色度图(chromaticity diagram)(图自 维基):

wikipeida-CIE1931xy

Cartopy 系列:裁剪填色图出界问题

前言

裁剪或者说白化,就是让填色图只显示在多边形里面,不显示在多边形外面,例如只显示 GeoAxes.contourf 在中国境内的结果。实现方法为:

from matplotlib.path import Path
from cartopy.mpl.patch import geos_to_path
from cartopy.io.shapereader import Reader

reader = Reader(filepath)
geom = next(reader.geometries())
reader.close()

cf = ax.contourf(X, Y, Z, transform=crs)
geom = ax.projection.project_geometry(geom, crs)
path = Path.make_compound_path(*geos_to_path(geom))
for col in cf.collections:
    col.set_clip_path(path, ax.transData)
  • crs 坐标系上的多边形对象变换到 data 坐标系上。
  • 利用 geos_to_pathmake_compound_path 将变换后的多边形转为 Path 对象。
  • QuadContourSet.collections 里的每个成员调用 set_clip_path 方法,并且指定 data 坐标系。

fig1

CALIPSO L2 VFM 产品的读取和绘制(with Python)

前言

CALIPSO 卫星的 L2 VFM(Vertical Feature Mask)产品根据激光的后向散射和消光信息,将激光通过的各高度层分类为云或气溶胶。该产品在现实中的表现如下图所示:卫星一边在轨道上移动一边向地面发射激光脉冲,相当于在地面上缓缓拉开一幅“画卷”,VFM 描述了“画卷”上云和气溶胶的分布和分类情况。

R-C

处理 VFM 产品的难点在于:

  • VFM 数组呈 (N, 5515) 的形状,N 表示卫星移动时产生了 N 次观测,但 5515 并非表示有 5515 层高度,而是三种水平和垂直分辨率都不同的数据摊平成了长 5515 的数组。因此处理数据时需要参照文档的说明对 5515 进行变形。
  • 文件中的经纬度和时间与 5515 的对应关系。时间数组需要解析成可用的格式。
  • 每个 range bin 的分类结果编码到了 16 位的无符号短整型的每个比特上,需要按位解码。
  • 网上现成的代码偏少。

网上能找到的代码有:

笔者也曾写过两次教程:

本文是对旧教程的翻新,会对 VFM 数据的结构进行更多解释,对代码也进行了更新。本文使用 pyhdf 读取 HDF4 文件,用 Matplotlib 3.6.2 画图。为了方便画图,用了一些自制的函数(frykit)。虽然基于 Python,但希望能给使用其它语言的读者提供一点思路。

完整代码已放入仓库 calipso-vfm-visualization

Matplotlib 系列:手动设置时间序列折线图的刻度

前言

Matplotlib 中画折线图用 ax.plot(x, y),当横坐标 x 是时间数组时,例如 datetimenp.datetime64 构成的列表,xy 的组合即一条时间序列。Matplotlib 能直接画出时间序列,并自动设置刻度。下面以一条长三年的气温时间序列为例:

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

df = pd.read_csv('test.csv', index_col=0, parse_dates=True)
series = df.loc['2012':'2014', 'T']

fig, ax = plt.subplots(figsize=(10, 4))
ax.plot(series.index, series)
ax.set_ylabel('Temperature (℃)')

print(ax.xaxis.get_major_locator())
print(ax.xaxis.get_major_formatter())
<matplotlib.dates.AutoDateLocator object at 0x000001AC6BF89A00>
<matplotlib.dates.AutoDateFormatter object at 0x000001AC6BF89B20>

fig_1

编写科研作图代码有更好的方法

这是物理海洋学家 Ken Hughes 在 2021 年发表的博客文章,原文标题为 A better way to code up scientific figures。以 Matplotlib 和 Matlab 为例,强调了模块化思想对于科研作图代码的帮助。我很少看到关于作图方法论的文章,所以翻译出来交流学习。

我画一张出版级别的科研配图一般需要写 100 - 200 行代码,这个长度有点点危险,因为很容易写出能正常运行但又一团糟的东西。如果代码片段都很短还可以从头重写,但如果代码有上千行,提前做好规划会更明智一些。不过在这两种极端情况之间潜藏着另一种吸引人的做法:写出一段当时感觉无比连贯,但以后会让你吃苦头的脚本。

假设你想画一张中等复杂度的图片,类似下面这张:

data_overview-1

Cartopy 系列:探索 shapefile

前言

Cartopy 可以通过 feature 模块向地图添加国界 BORDER 和省界 STATES,因其底层采用的 Natural Earth 地图数据并不符合我国的政治主张,所以我们经常需要自备 shapefile 文件来画中国省界,以下面的代码为例

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.io.shapereader as shpreader

extents = [70, 140, 0, 60]
crs = ccrs.PlateCarree()
fig = plt.figure()
ax = fig.add_subplot(111, projection=crs)
ax.set_extent(extents, crs)

filepath = './data/bou2_4/bou2_4p.shp'
reader = shpreader.Reader(filepath)
geoms = reader.geometries()
ax.add_geometries(geoms, crs, lw=0.5, fc='none')
reader.close()

plt.show()

图就不放了,这段代码足以应付大部分需要画省界的情况。然而我在无脑粘贴代码的过程中逐渐产生了疑惑:为什么 shapefile 会由三个文件组成?省界是以何种形式存储在文件中?Cartopy 和 Matplotlib 又是怎样将省界画出来的?调查一番源码后总结出了这段代码底层实现的流程:

  • 利用 PyShp 包读取 shapefile 文件中的每个形状。
  • 利用 Shapely 包将形状转换为几何对象。
  • 利用 Cartopy 包将几何对象投影到地图所在的坐标系上。
  • 用投影后的坐标构造 Matplotlib 的 Path 对象,最后画在地图上。

本文的目的即是从头到尾解说一下这段流程,希望加深对 shapefile 格式,Matplotlib 和 Cartopy 包的理解。令人意外的是,随着探索的不断深入,我发现自己自然而然地学会了如何实现省份填色、省份合并,地图白化等,以前看起来十分困难的操作。本文也会一并介绍这些应用。

Matplotlib 系列:网格数据与 pcolor

前言

Matplotlib 的 pcolor 函数能够绘制由一个个四边形(quadrilateral)单元构成的网格数据的彩色图像,相比绘制等值填色图的 contourf 函数,不会产生过度的平滑效果,能忠实反映像元的数值大小,因而在科学可视化中也很常用。本文并不打算介绍该函数的种种,只想着重讨论网格数据的显示效果、shading 参数发挥的作用,以及 pcolorpcolormesh 这对双胞胎间的差异。本文基于 Matplotlib 3.3.4。

图解网格数据

pcolor 全名 pseudo color,即伪彩色。函数签名为

pcolor([X, Y], C, **kw)

其中 XY 分别是网格的横纵坐标,C 是网格单元内变量的数值。之所以称之为“伪”,是因为 pcolor 并不像 imshow 那样直接用 RGB(A) 数组表示颜色,而是将 C 的数值归一化之后,在一个颜色查找表中查找对应的颜色,进而用颜色差异表现数值大小(原理详见 Matplotlib 系列:colormap 的设置)。C 数组的形状为 (ny, nx)XY 的形状要比 C 大上一圈,即 (ny + 1, nx + 1)ny 在前表示纵坐标会随数组的行号变动,nx 在后表示横坐标会随数组的列号变动。pcolor 对网格数据的显示效果如下图所示

regular_and_irregular_grids

Matplotlib 系列:图解 quiver

前言

Matplotlib 中用箭头表示风场或电磁场等矢量场时需要用到 quiver 方法,据字典,quiver 一词的意思是颤动、颤抖或箭袋,貌似也就最后一个意思跟箭头搭得上边。相比于其它画图方法,quiver 的参数又多又容易混淆,所以本文将以图解的方式逐一介绍。这些参数按功能可分为三种:控制箭头位置和数值的、控制箭头长度和角度的,以及控制箭头尺寸和形状的。下面会按照这个分组顺序来解说。本文代码基于 Matplotlib 3.3.4。