Introduction
在引入了数组概念之后，上周我们进一步了解了如何创建数组和基本的数组计算方法。显然，我们上一期的内容基本都是针对所有数组元素进行的。
那么当我们需要单独计算数组中特定位置或特定区域的值时，又该如何进行操作呢？对应到地球科学中，当我们拥有全球的格网数据，一般而言都需要将我们研究区对应的部分提取出来。
以及另一个问题，我们之前的数组都是在创建时便指定了其维度大小的。但现实是，很多数组可能是分开储存，需要最终聚合在一起才能实现我们所需要的功能的。这又应当如何操作？
❝一个栗子
我们的逐月气温数据分别储存在不同的数组中，我们计算年平均气温时则需要将其聚合成一个数组求均值。
或许你会说，反正只有12个月，那我可以分别读取后从第1项加到第12项，再除以12即可。
倒也不是不行。但是如果有1200个月，足下又该如何应对？
然后，你说，「那我可就要用到循环了」。
不好意思，循环部分我们还没有讲到，我先默认我们还不会。
另外，即使循环可以简单地实现均值，但我们一般会有更多复杂的参数需要计算。
例如，求取每个格网的时间变化趋势，也就是一次函数中的斜率k。对于一堆零散的数组，我不认为使用循环可以简单地解决这个问题。
因此，数组重构尤为重要。当然并不止这一个浅显的例子，它们都有更广阔的应用空间，这些都可以之后慢慢探索。
个人感觉我们将储存了地学数据数组的部分内容提取出来实际上也是对数据结构的一种改变，因此这里的标题也就没有索引的一席之地了。
❝然后标题很奇怪，最后面有个½。那是因为写稿的时候发现，自己呼呼啦啦写了一堆乱七八糟的，第一节长度就已经严重超标了(低估自己胡说八道的能力了，属于是)。
那么，就拆成两期更吧！
索引、切片
索引
前面我们提到过索引的概念：索引(index)，指的是获取指定位置的元素，列表索引方式与数组类似。
并且，由于Python遵循C语言的语法习惯，索引值是从0起算，这与MATLAB、R和Julia等偏数学的语言不一致，需要注意。而当如果数组有多个维度时，则需要使用多个下标来访问特定元素。
值得一提的是，数组的索引分为正序和倒序。当我们按照从头到尾的顺序索引时，则从0起算以确定元素对应的下标；如果我们的元素显然从后往前更方便确定其位置，我们不妨从最后一个元素起算，从后往前依次为-1，-2，-3，以此类推。
下面我们展示案例：
import numpy as np
arr = np.array([1, 2, 3, 4, 5])
print(arr[0])
print(arr[4])
print(arr[-1])
print(arr[-3])
1
5
5
3
当数组有多个维度时，我们需要使用多个下标来访问特定元素。
例如，对于一个二维数组，我们可以用arr[i][j]的方式进行索引，其中i表示行索引，j表示列索引。
arr = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
print(arr[0])
# 以下两种写法最终效果一致
print(arr[0][0])
print(arr[0, 0])
print(arr[-1, -3])
[1 2 3]
1
1
7
切片
显然索引存在一个问题，我们每次只能获取一个位置的值，那么有没有那种更高级的，一次能切下一大块区域组成一个新数组的方式呢？当然，那便是就是更高阶的索引方式——切片。
常规切片
切片，就是将索引值位置修改为使用冒号指定的元素范围。
其语法结构为array[start:stop:step]，意思是从索引值start位置到stop位置以step为步长进行索引，step默认为1，可以不写。
特别注意，切片的索引方式为左闭右开的区间，最后一个索引值的元素是取不到的。
arr = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
print(arr[1:5]) # 正向切片索引值1-4的元素
print(arr[-5:-1]) # 反向切片索引值-5到-2的元素
print(arr[0:10:2]) # 步长为2正向切片
[2 3 4 5]
[6 7 8 9]
[1 3 5 7 9]
arr = np.random.rand(180, 360) # 定义一个1°x1°的格网数据
print(arr[30:61, 60:101]) # 提取经纬度范围为30°N-60°N，60°N-100°E的格网数据
[[0.62427043 0.82133876 0.50815278 ... 0.43005637 0.84763936 0.99592113]
[0.01961363 0.7792136 0.83044477 ... 0.08328295 0.24768328 0.02536089]
[0.87679775 0.81019605 0.31762462 ... 0.12467192 0.74212048 0.59793829]
...
[0.72432394 0.10397449 0.00269136 ... 0.79870904 0.36223087 0.88735507]
[0.29034297 0.05487801 0.1659494 ... 0.4014151 0.41991349 0.85197025]
[0.66423806 0.94671443 0.24224496 ... 0.46171495 0.13174959 0.21355304]]
那么问题来了，如果我们切片时需要取到最后一个元素，stop这里填-1显然是不对的。0？那更不对。那岂不是每次都得用数组的size属性，再将它+1才能确保能取到最后一个元素？当然不会。如果我们切片时选择不填start和stop，那么它就默认能够取到两端的元素。
print(arr[::2]) # 步长为2正向切片
print(arr[::-1]) # 翻转数组
print(arr[1:8:3]) # 步长为3正向切片索引值1-7的元素
[1 3 5 7 9]
[10 9 8 7 6 5 4 3 2 1]
[2 5 8]
当然，如果你有奇数个元素，并且你按照array[::2]的方式切片。显然，最后一个元素是不会出现的。但并不是它没有被取到，而是步长为2跳过了它。
arr = np.array([1, 2, 3, 4, 5, 6, 7, 8])
print(arr[::2])
arr = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9])
print(arr[::2])
[1 3 5 7]
[1 3 5 7 9]
布尔切片
或许你想说，需要知道固定位置的切片方式也过于死板了。有没有那种更高级的，比如我要筛选气温大于30°C的格网所对应区域的降水？
于是，基于逻辑运算的布尔索引就出现了。当我们判断一个数组与一个条件的关系时，数组将会变为值为True和False的掩膜数组。
那么我们就可以用这个掩膜去实现更复杂的提取了。
temperature = np.random.rand(180, 360) * 100 # 定义一个1°x1°的气温格网数据
another_index = np.random.rand(180, 360) # 定义一个1°x1°的指标网格数据
temperature_mask = temperature > 30 # 判断气温是否高于30°C
print(temperature_mask)
print(another_index[temperature_mask]) # 使用掩膜数组索引另一个形状一致的数组
[[False True True ... True False True]
[ True True True ... False True False]
[ True True False ... True False True]
...
[ True True True ... True True False]
[ True True True ... True True True]
[ True False True ... True True True]]
[0.23491348 0.56643524 0.1230861 ... 0.97229519 0.05542499 0.35975104]
特殊切片
此外，最后我们再提一下一种指定不规则位置的特殊索引，以备应对一些可能的特殊情况。
arr = np.array([[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]])
print(arr)
print(arr[[0, 1, 1], [0, 2, 8]])
[[ 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]]
[ 1 13 19]
数组赋值
说了这么多使用索引、切片获取指定位置值，实际上它们还为修改数组元素的值提供无限可能。我们不妨给出一些实际的例子：
不使用索引、切片的情况下：我们只能对整个数组进行操作，如将摄氏度转换为开尔文。
temperature_celsius = np.random.rand(180, 360) * 100
temperature_kelvin = temperature_celsius + 273.15
print(temperature_celsius)
print(temperature_kelvin)
[[46.98955793 48.52445665 87.89904253 ... 11.81341333 16.13229844
80.53932617]
[36.50574919 44.59127172 13.77015774 ... 72.37198475 1.29320209
45.19339505]
[47.77012955 69.52877456 60.50566262 ... 14.42855063 53.12931476
86.16625651]
...
[10.14869845 90.79185507 15.7318585 ... 13.588911 8.37859764
52.9307133 ]
[88.98887715 25.7914286 92.92130093 ... 83.04135583 21.70086893
9.560699 ]
[88.43177587 92.90246587 52.40783702 ... 72.79336751 55.33133393
15.63278645]]
[[320.13955793 321.67445665 361.04904253 ... 284.96341333 289.28229844
353.68932617]
[309.65574919 317.74127172 286.92015774 ... 345.52198475 274.44320209
318.34339505]
[320.92012955 342.67877456 333.65566262 ... 287.57855063 326.27931476
359.31625651]
...
[283.29869845 363.94185507 288.8818585 ... 286.738911 281.52859764
326.0807133 ]
[362.13887715 298.9414286 366.07130093 ... 356.19135583 294.85086893
282.710699 ]
[361.58177587 366.05246587 325.55783702 ... 345.94336751 328.48133393
288.78278645]]
当我们使用索引：已知某个点气温检测到异常偏高，我们可以索引出来后减去异常值。
anomaly = 1 # 定义异常值
print(temperature_celsius[12, 12])
temperature_celsius[12, 12] = temperature_celsius[12, 12] - anomaly
print(temperature_celsius[12, 12])
25.76296894744049
24.76296894744049
当我们使用切片：已知某个区域气温存在再分析数据总体偏高，于是也将其切片出来统一减去异常值。
anomaly = 10
print(temperature_celsius[20:30, 50:60])
temperature_celsius[20:30, 50:60] = temperature_celsius[20:30, 50:60] - anomaly
print(temperature_celsius[20:30, 50:60])
[[25.3579225 83.68788079 96.68081625 90.43210713 76.31520603 78.04790683
54.73767051 7.37831617 54.25792131 89.02225569]
[59.47570906 56.72062105 12.01664326 11.62892288 10.2361407 43.45036545
60.57784347 -0.1711559 35.82533646 86.56153631]
[48.30386448 33.84777879 48.16699698 23.18488938 37.07460779 78.00805894
44.74629321 62.80326902 49.66957877 2.95878244]
[93.10067452 37.61881634 32.04477007 47.14024789 57.61126134 91.2445765
92.85022329 17.74543191 68.76463834 49.70241831]
[ 5.53869259 63.08667929 83.09678338 43.20958447 26.50022741 40.53110859
26.65897495 6.16044341 35.93270875 11.84789225]
[68.42386338 0.27629717 76.34937605 93.75288009 80.91631945 80.7780009
46.0720071 18.1553839 78.5540843 16.10385965]
[78.83995998 76.53035252 33.78898854 74.0211111 78.20355401 20.91878273
97.58933754 1.98253277 68.44972986 53.22561693]
[59.37270815 98.28497602 61.51703955 74.10907248 14.81775848 27.5519606
14.95885246 86.3130667 85.31253824 39.78428638]
[39.1210945 0.29862888 1.23686996 44.10488919 3.94744308 52.64460936
54.84865332 32.83740901 86.03344693 78.71924222]
[15.60075556 73.35991039 19.80827081 72.00264234 50.8875262 65.9904218
36.15142112 58.72433171 77.55760775 93.11467651]]
[[ 15.3579225 73.68788079 86.68081625 80.43210713 66.31520603
68.04790683 44.73767051 -2.62168383 44.25792131 79.02225569]
[ 49.47570906 46.72062105 2.01664326 1.62892288 0.2361407
33.45036545 50.57784347 -10.1711559 25.82533646 76.56153631]
[ 38.30386448 23.84777879 38.16699698 13.18488938 27.07460779
68.00805894 34.74629321 52.80326902 39.66957877 -7.04121756]
[ 83.10067452 27.61881634 22.04477007 37.14024789 47.61126134
81.2445765 82.85022329 7.74543191 58.76463834 39.70241831]
[ -4.46130741 53.08667929 73.09678338 33.20958447 16.50022741
30.53110859 16.65897495 -3.83955659 25.93270875 1.84789225]
[ 58.42386338 -9.72370283 66.34937605 83.75288009 70.91631945
70.7780009 36.0720071 8.1553839 68.5540843 6.10385965]
[ 68.83995998 66.53035252 23.78898854 64.0211111 68.20355401
10.91878273 87.58933754 -8.01746723 58.44972986 43.22561693]
[ 49.37270815 88.28497602 51.51703955 64.10907248 4.81775848
17.5519606 4.95885246 76.3130667 75.31253824 29.78428638]
[ 29.1210945 -9.70137112 -8.76313004 34.10488919 -6.05255692
42.64460936 44.84865332 22.83740901 76.03344693 68.71924222]
[ 5.60075556 63.35991039 9.80827081 62.00264234 40.8875262
55.9904218 26.15142112 48.72433171 67.55760775 83.11467651]]
一些更复杂的应用：假定气温大于50°C时某一指标便会出现异常而超出临界值，我们将它们掩膜出来并使用边界条件进行限制。
another_index = np.random.rand(180, 360) * 10 # 定义一个1°x1°的指标网格数据
print(another_index[temperature_celsius >= 50])
boundary_condition = 0.5
another_index[temperature_celsius >= 50] = boundary_condition # 对气温将大于50°C的值设置为边界条件
print(another_index[temperature_celsius >= 50])
[3.31433926 7.22395752 6.10150125 ... 0.74224767 0.4194778 8.52261236]
[0.5 0.5 0.5 ... 0.5 0.5 0.5]
于是，上一期我们讲到的处理空值的函数nan_to_num便成了一种摆设，我们完全可以自己实现它：
arr = np.array([1, 2, np.nan, 4, 5])
arr[np.isnan(arr)] = -1000
print(arr)
[ 1. 2. -1000. 4. 5.]
Supplementary
另外上次的数组计算部分中，个人认为还有几个比较实用的函数没有提及，这里进行补充。
矩阵乘法.dot：
arr0 = np.array([[1, 2], [3, 4]])
arr1 = np.array([[5, 6], [7, 8]])
print(arr0.dot(arr1))
[[19 22]
[43 50]]
MATLAB中数组会默认进行矩阵计算，不同数组一一对应的运算需要加点进行区分(如，.*，./等)。而在NumPy中则默认使用了更常用的一一对应运算，而将矩阵运算进行了函数封装。
角度制弧度制转换deg2rad和rad2deg：
arr = np.array([0, 30, 45, 60, 90]) # 定义角度数组
print(np.sin(arr * np.pi / 180)) # 角度转弧度，并计算正弦值
print(np.deg2rad(arr)) # 内置函数转换
print(np.sin(np.deg2rad(arr)))
print(np.sin(np.deg2rad(np.rad2deg(np.deg2rad(arr))))) # 俄罗斯套娃一下
[0. 0.5 0.70710678 0.8660254 1. ]
[0. 0.52359878 0.78539816 1.04719755 1.57079633]
[0. 0.5 0.70710678 0.8660254 1. ]
[0. 0.5 0.70710678 0.8660254 1. ]
此处2即英文to的谐音，这与MATLAB中数据类型转换的函数命名风格一致，知道的同学应该已经理解了。
关于弧度制与角度制的转换当然我们可以根据它们的公式定义来计算，但是为了避免有时候就是死活想不起来公式长啥样，还是用上NumPy为我们封装好的“语法糖”吧。
后记
刚刚提到了一个词，叫“语法糖”。所谓语法糖，就是一些封装好的常用功能，广义上来说，所有使用的函数都可以称为语法糖。
但在实际使用过程中，如果我们追求创新，应当是很难完全借助语法糖来实现的。所以，我们要深入理解基本原理与对应的数学关系，才能从根本上构建自己的模型框架。
虽然我们总说不要重复造轮子。这也是Python兴盛的原因之一，你总能找到一些与相关的第三方库去实现一些功能。但如果需要构建一个足够复杂、完善、通用的模型框架，或许还是需要我们每天一铲子，从头开始搭起。
从数学原理上出发，才能实现更多的Innovation和Complexity，而避免沦为过度依赖第三方库的“调包侠”。
Manuscript: RitasCake
Proof: Philero; RitasCake
获取更多资讯，欢迎订阅微信公众号：Westerlies
7万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
