SC.NumPy 01 | 引入、概念

Introduction

NumPy是 Python 科学计算(Scientific Calculation，SC，也就是本系列标题中缩写字母的含义) 生态系统的基础与最重要的组成部分(窃以为“没有之一”)，同时也是 Python 数据分析领域的“三剑客”之一(NumPy 为数值计算的核心库，Pandas 为强大的数据分析和操作工具，Matplotlib 则提供了可视化能力)。

❝

实际使用过程中，个人感觉 SciPy 在数分领域也具有普遍的适用性，因此本人更习惯将 SciPy 也列为“数分三剑客”之一。众所周知，“御三家”通常有四个，就像“四大天王”有五个，因此“数分三剑客”有四个也没什么奇怪的，论证完毕。

NumPy 的核心在于支持大规模高阶数组(N-dimensional array，NumPy 类型缩写为 ndarray，或者我们直接称之为 array)和矩阵运算，以及与之适配的大量数组操作 API。

从最简单的四则运算到我们使用 Python 构建起庞大的数值模型，NumPy 都将是这座城堡的基石。所以，我们的故事就从 NumPy 开始讲起吧。

列表与数组

列表

在讲解数组之前，有必要先提及列表(List)这一原生 Python 中最常见的容器。

通俗而言，列表就是将一系列元素按顺序装进一个容器中，而这个容器就是我们所说的列表。这些元素可以是任何东西，他可能是一个数字、一串单词组成的字符串、一个 Excel 表格的内容，甚至可以是另一个列表(从组织形式而言，私以为嵌套列表与数组间没有本质区别)。

列表创建时不需要指定其长度，而可以通过append等内置函数便捷地对元素进行增/删/查/改等操作。

❝

Python 中通常使用方括号[]，并用逗号,分隔所需存储的元素来创建列表[1, 'Los Angeles', True, [2, 'Sheffield', 'How are u?']]。

数组

数组是相同数据类型元素的集合，存储在连续的内存位置。每个数组元素都由一个或一组索引标识，可有效地存储和操作大规模数据。

对于数组而言，其创建时须指定维度大小，因而具有固定长度(我们可以分别使用ndim，shape，size和dtype获取数组的维度、形状、大小和数据类型)。数组元素的增/删/改均需赋值为新的数组类型变量才能完成操作。此外，同一个数组中的元素类型应当相同，如果该数组已经被指定为数值型，那么将一个字符型的元素赋值到数组中的操作将不可用。

import numpy as np

# 使用随机数创建一个2000×3000的二维数组(文中未提及的函数均会做简单说明，具体将在后续内容中介绍)
arr = np.random.rand(2000, 3000)

print(arr.ndim)     # 维度
print(arr.shape)    # 形状
print(arr.size)     # 大小
print(arr.dtype)    # 数据类型

arr[100, 200] = 'Python'    # 类型不符的赋值操作将失败

2
(2000, 3000)
6000000
float64

ValueError: could not convert string to float: 'Python'

ls = [1, 'Los Angeles', True, [2, 'Sheffield', 'How are u?']]
ls[0] = 'Zürich'    # 列表可接受不同类型赋值

print(ls)

['Zürich', 'Los Angeles', True, [2, 'Sheffield', 'How are u?']]

那么相较于列表，数组的优势在哪儿？

数组支持广播运算和向量化操作使得它的运行速度和代码量均大幅优化，同时 NumPy 提供了大量用于执行高级数学运算的函数，可以便捷地进行如矩阵乘法、傅里叶变换和微积分等运算。

NumPy 相对于列表的运算速率，我们可以进行一个简单的试验：

import time

t1 = time.time()    # 获取当前时间戳
x = np.arange(1000000) + 1    # 将1000000以内的所有自然数+1
t2 = time.time()

print('Runtime: %.2f ms' % ((t2 - t1) * 1000))  # 计算运行时间

Runtime: 2.00 ms

t1 = time.time()

ls = []
for i in range(1000000):
    ls.append(i + 1)

t2 = time.time()
print('Runtime: %.2f ms' % ((t2 - t1) * 1000))

Runtime: 96.00 ms

显然，数组的运行速率显著优于列表。

数组轴

前面我们提到「从组织形式而言，嵌套列表与数组间没有本质区别」，这里就涉及到了维度dim的概念。

NumPy 中每个维度均对应了一条轴axis(其中第一个维度为axis=0，第二个维度为 1，以此类推)，这使得更丰富的操作变得可能，如分别对每一行/每一列分别计算均值/最值等。

arr = np.random.rand(2000, 3000)

print(arr.mean())    # 对所有值求均值
print(arr.sum(axis=0))      # 按行求和
print(arr.mean(axis=1))     # 按列求均值

0.4999523616278683
[999.26900252  994.98832343  994.0841875  ... 1025.45738118  994.97250643 992.98966004]
[0.49987789 0.50101981 0.49383497 ... 0.50501385 0.49313256 0.50345248]

值得注意的是，数组轴的顺序并非固定不变，这或许与我们的常理相悖，但其中依然有其自洽的逻辑。

NumPy 中，最基础的数组单位为行向量(一维数组)。当我们创建一个 0-9 的自然数向量arr = array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])，由于所有元素均在第 0 行，此时并不存在行的概念，仅存在一个axis=0的轴(即我们常说的 x 轴)，索引值对应元素所在列(如，arr[2]的值即为2)。

通过在行的方向上堆叠多个行向量形成二维数组，行的概念便产生了，于是axis 推广为 2 个(0 和 1，或者说 y 和 x 轴)。但由于行索引相较列索引高一级(因为行是由堆叠行向量才产生的，而 NumPy 的设计中索引方式是由外而内)，y 对应的轴占据了axis=0，x 对应的轴则相应往后推。

多维数组的产生以此类推。

arr1D = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
arr2D = np.array([[ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9],
                  [10, 11, 12, 13, 14, 15, 16, 17, 18, 19],
                  [20, 21, 22, 23, 24, 25, 26, 27, 28, 29],
                  [30, 31, 32, 33, 34, 35, 36, 37, 38, 39]])

print(arr1D[2])
print(arr2D[0, 2])

2
2

空间坐标系中，通常使用(x, y, z)来表示一个元素的位置。而在 NumPy 中，当我们已知二维数组中行和列分别是 y 和 x 轴(或者说，其axis 分别为 0 和 1)，通过堆叠多个二维数组组成三维数组时，三维数组中元素的索引方式通常是array[z, y, x](此时，z 轴的 axis为 0，而行、列则变为 1 和 2)。

这是前文中一维升二维情况的三维推广，当我们通过堆叠多个二维数组组成三维数组，索引元素时首先需要确定该元素在哪一个二维数组中，在确定二维数组所在层级后再对其行列位置索引。

❝

索引(index)，指的是获取指定位置的元素，如二维数组A中第 2 行第 3 列位置的元素可以通过A[2, 3]或A[2][3]索引，列表索引方式与数组类似。并非本节核心内容，后续将进一步讲解，此处不赘述。

需要注意的是，Python 中的索引值与 C 等传统程序设计语言一致，均是从 0 开始。上面的第 2 行第 3 列代表的是数学意义上第 3 行第 4 列位置的值。

实际上，A[2, 3]和A[2][3]的索引方式存在一定的差异，阅读完本文后可以思考一下它们的差异在何处？

arr0 = np.random.rand(1000, 800)
arr1 = np.random.rand(1000, 800)
arr3d = np.stack((arr0, arr1))     # 沿新轴堆叠，默认为axis=0

print(arr3d.shape)

(2, 1000, 800)

再回到本节最开始的那句话，假定我们将多个普通非嵌套元素组成的列表List1D嵌套组合成一个新列表，此时列表中不存在列表、字典等可嵌套索引的容器，新列表List2D的数据结构应当与二维数组相似。

而我们进一步将多个嵌套后的列表二阶嵌套，类似三维数组的列表List3D便也形成。此时，如果我们需要索引(x0，y0，z0)处的元素，我们可以通过类似数组索引的方式实现List3D[z0][y0][x0]。

list1d = [1, 2, 3, 4]
list2d = [list1d, list1d, list1d, list1d]
list3d = [list2d, list2d]

print(list3d)
print(list2d[1][2])
print(list3d[0][1][2])

[[[1, 2, 3, 4], [1, 2, 3, 4], [1, 2, 3, 4], [1, 2, 3, 4]], [[1, 2, 3, 4], [1, 2, 3, 4], [1, 2, 3, 4], [1, 2, 3, 4]]]
3
3

应用案例

如果单纯从数学角度理解或许以上内容显得晦涩而抽象，我们不妨给出一个贴近地球科学的案例。

假定我们拥有一个庞大的气温数据集，其维度大小为 100×8×180×360(此处使用 NumPy 生成的随机数替代)，分别对应时间/气压层高度/纬度/经度(time/level/lat/lon)。

temp_pres_level = np.random.rand(100, 8, 180, 360) * 10
print(temp_pres_level.shape)

(100, 8, 180, 360)

其中，最基础的维度为经度(lon，长度为 360)和纬度(lat，长度为 180)，数组中任意一个空间分布表征了某一时间、某一气压场高度的气温场状态。

# 时段10气压层3的气温空间分布(均以0起算，不再赘述)
print(temp_pres_level[10, 3, :, :])
print(temp_pres_level[10, 3, :, :].shape)

[[5.15967091 8.32741165 8.38750895 ... 5.13925573 6.31967397 3.6473191 ]
 [2.43346555 1.47431224 3.64297239 ... 1.16419275 2.16354937 0.56699932]
 [8.49879345 9.36592812 4.35082365 ... 7.25290327 7.51761726 1.49693121]
 ...
 [0.06623796 4.7796034  9.29730585 ... 8.81317218 3.70472341 0.65723977]
 [8.36701033 5.70102812 8.31637686 ... 0.65968372 1.86492099 2.29193739]
 [3.33517635 6.73889334 4.89346057 ... 2.65087397 7.88711282 0.84554989]]

(180, 360)

任一时间切片中，均存在由该时间段 8 个气压层高度上的气温场堆叠而成的气温-高度三维空间分布。

# 时段10的三维气温场
print(temp_pres_level[10, :, :, :])
print(temp_pres_level[10, :, :, :].shape)

[[[7.2727462  5.69704542 4.36228912 ... 0.01613352 4.09021786 1.49154325]
  [4.72277386 6.4618104  8.19231906 ... 6.10729153 8.25230224 2.4833938 ]
  [6.24862484 0.52868575 3.56514689 ... 4.93829487 1.80267463 1.38778563]
  ...
  [3.44555925 2.83745865 2.86904983 ... 0.97567335 9.02052484 9.77716271]
  [0.07316817 8.18822885 4.06921978 ... 5.89807453 2.29668037 7.08444082]
  [0.01961915 9.70148522 5.26756195 ... 7.89959187 3.69218049 2.27333866]]

 ...

 [[3.72730666 3.83916958 2.72706794 ... 4.78328168 0.25756904 0.17664715]
  [6.1722617  8.44116822 6.61624912 ... 2.26677992 4.59884788 7.84832724]
  [7.7379857  3.7521109  7.10844922 ... 5.38552642 3.82394269 2.73134341]
  ...
  [5.02510603 5.72130781 1.14054717 ... 5.76591574 9.38298676 3.84753024]
  [0.90644445 9.5964141  6.82524456 ... 2.30565578 7.56616989 1.36951395]
  [2.43210955 1.84263774 4.02661893 ... 7.22185607 4.18110154 1.39632165]]]

(8, 180, 360)

当我们指定各个维度的索引值，我们就能得到特定时段/气压层高度/经度/纬度的气温。

# 时段10气压层0纬度90经度180位置气温值
print(temp_pres_level[10, 0, 90, 180])

4.262264606486996

因此，我们不难获得大量气候要素的分布。

print(temp_pres_level.mean(axis=0))   # 时间平均气温(输出过长，该部分省略)
print(temp_pres_level.mean(axis=1))   # 气压层平均气温
print(temp_pres_level.min(axis=3))   # 纬圈最低气温
print(temp_pres_level.max(axis=2))   # 经圈最高气温
print(temp_pres_level.mean(axis=(2, 3)))  # 时间-气压层的空间平均气温
print(temp_pres_level[10].mean(axis=1))   # 时间段10的气压层-经度平均气温分布
print(temp_pres_level[10].mean(axis=2))   # 时间段10的气压层-纬度平均气温分布

结语

鸽了不知道多久的公众号，终于写出了第一篇(所有推送内容均同步至和鲸社区、腾讯云与CSDN，欢迎关注)。属实费老大劲儿了，后面熟悉了应该能高效一点。从设计风格到主题样式，从主体内容到行文结构，很多东西比预想的复杂，当然也有些部件实际上难度被高估了。

「万事开头难」，虽说中间和结尾也难说会简单到哪里去，但总算是开始填坑了。当然，未来也可能会挖很多最终都没填完的坑。更多关于创建公众号的过程和奇思妙想欢迎移步『序章』听我絮絮叨叨。

另外，欢迎大家留言评论哈。根据大家提出的比较主流的问题适当地出一些番外篇，个人感觉也是不错的突破创作瓶颈的方式。

那么，咱们下周见(也许)。

Manuscript: RitasCake

Proof: Philero; RitasCake

获取更多资讯，欢迎订阅微信公众号：Westerlies