import numpy as np3 NumPy and Xarray
3.1 Topics Covered
- NumPy data types
- Data type, size, shape, and number of dimensions of an array
- Creating, reshaping, and slicing NumPy arrays
- Numeric and comparison operations on arrays
- How Xarray differs from NumPy
- Creating, subsetting, and summarizing labeled arrays with Xarray
3.2 Introduction
The NumPy and Pandas libraries are central to data science in Python. NumPy allows for the efficient analysis and processing of data arrays with varying sizes, shapes, and numbers of dimensions while Pandas allows for reading in and working with data tables. In this chapter, we focus on NumPy. We also introduce Xarray, which builds on NumPy by adding labeled dimensions and coordinates to arrays. Pandas is discussed in Chapter 4.
3.3 Creating NumPy Arrays
The NumPy library allows for creating and working with arrays. It is of specific value when you want to perform mathematical operations on values stored in arrays. It is also very fast and memory efficient due to its reliance on the C language. As mentioned in Chapter 1, arrays are similar to lists in that they store a series of values or elements. However, arrays can be expanded to include many dimensions. For example, an image could be represented as an array with three dimensions: height, width, and channels. If you work with deep learning, tensors are the primary data model used to read and manipulate data and are essentially multidimensional arrays that can be stored in RAM or within GPU memory for faster computation. In short, array-based calculations and manipulations are essential to data science, so NumPy is an important library to learn if you work in the Python environment.
The complete documentation for NumPy can be found here.
Before you can use NumPy, you must make sure that it is installed into your conda environment, as described in the Preface. Once NumPy is installed, you will need to import it in order to use it in your code. It is common to assign NumPy an alias of np to simplify your code.
Lists can be converted to NumPy arrays using the array() function. Once the list object is converted to an array, the type is defined as numpy.ndarray, which indicates that it is a NumPy array specifically. Since this array only has one dimension, it is specifically called a vector.
lst1 = [3, 6, 7, 8, 9]
arr1 = np.array(lst1)
print(type(lst1))<class 'list'>print(type(arr1))<class 'numpy.ndarray'>print(arr1)[3 6 7 8 9]A two-dimensional array is known as a matrix. In the example below, we generate a matrix from a list of lists.
lst2 = [[3, 6, 7, 8, 9], [3, 6, 7, 8, 9], [3, 6, 7, 8, 9]]
arr2 = np.array(lst2)
print(arr2)[[3 6 7 8 9]
[3 6 7 8 9]
[3 6 7 8 9]]Again, one of the powerful advantages of NumPy arrays is the ability to store data in arrays with many dimensions. In the example below, we create a three-dimensional array from a list of lists of lists. This would be similar to an image with dimensions of image height, image width, and image channels (for example, red, green, and blue). A four-dimensional array could represent a time series (height, width, channels, and time) or a video containing multiple frames (frame height, frame width, channels, and frame number).
lst3 = [[[3, 6, 7, 8, 9], [3, 6, 7, 8, 9], [3, 6, 7, 8, 9]],
[[3, 6, 7, 8, 9], [3, 6, 7, 8, 9], [3, 6, 7, 8, 9]],
[[3, 6, 7, 8, 9], [3, 6, 7, 8, 9], [3, 6, 7, 8, 9]]]
arr3 = np.array(lst3)
print(arr3)[[[3 6 7 8 9]
[3 6 7 8 9]
[3 6 7 8 9]]
[[3 6 7 8 9]
[3 6 7 8 9]
[3 6 7 8 9]]
[[3 6 7 8 9]
[3 6 7 8 9]
[3 6 7 8 9]]]The cell below provides some examples of NumPy functions for generating arrays.
The arange() function returns an array of evenly spaced values and accepts start, stop, step, and data type parameters. In the example, we create an array of evenly spaced values from 0 to 100 with a step size of 5. We use 101 as opposed to 100 since the last value in the provided range is not included. We specifically define the data type as integer, but NumPy can infer a data type if it is not provided.
The linspace() function is similar to arange(); however, a number of samples is specified as opposed to a step size. In the example, since 5 samples are requested, 5 evenly spaced values between 0 and 100 are returned.
The ones() function is used to return an array of 1s. In the example, we generate a three-dimensional array where the first dimension has a length of 3, the second a length of 4, and the third a length of 4. The shape and dimensions of the array are specified using a tuple.
Similar to ones(), zeros() generates an array of zeros.
It is also possible to generate random values between 0 and 1 (random.rand()) and a specified number of random integer values between two values (random.randint()).
arr4 = np.arange(0, 101, 5, dtype="int")
print(arr4)[ 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85
90 95 100]arr5 = np.linspace(0, 100, 5, dtype="int")
print(arr5)[ 0 25 50 75 100]arr6 = np.ones((3, 4, 4))
print(arr6)[[[1. 1. 1. 1.]
[1. 1. 1. 1.]
[1. 1. 1. 1.]
[1. 1. 1. 1.]]
[[1. 1. 1. 1.]
[1. 1. 1. 1.]
[1. 1. 1. 1.]
[1. 1. 1. 1.]]
[[1. 1. 1. 1.]
[1. 1. 1. 1.]
[1. 1. 1. 1.]
[1. 1. 1. 1.]]]arr7 = np.zeros((3, 4, 2))
print(arr7)[[[0. 0.]
[0. 0.]
[0. 0.]
[0. 0.]]
[[0. 0.]
[0. 0.]
[0. 0.]
[0. 0.]]
[[0. 0.]
[0. 0.]
[0. 0.]
[0. 0.]]]arr8 = np.random.rand(3, 4, 5)
print(arr8)[[[0.88609588 0.83451541 0.51107762 0.36350961 0.7258199 ]
[0.40608208 0.24466723 0.16234896 0.23221965 0.09050055]
[0.13016631 0.58720394 0.15044583 0.21939141 0.78378492]
[0.14127942 0.46236845 0.84406729 0.84216381 0.63389861]]
[[0.25337333 0.76563718 0.65080751 0.23369803 0.09891386]
[0.20027669 0.81924743 0.74921102 0.10899274 0.90617973]
[0.39937862 0.1401036 0.58231025 0.74309508 0.60386222]
[0.83866278 0.20860347 0.83738251 0.75463004 0.43994647]]
[[0.83773455 0.13222103 0.37260129 0.20028061 0.60999762]
[0.18556061 0.12699235 0.88211062 0.8376589 0.42964443]
[0.33175242 0.3369349 0.18649587 0.03069557 0.53719916]
[0.08406271 0.95630543 0.51588475 0.10103673 0.72980306]]]arr9 = np.random.randint(1, 200, 7)
print(arr9)[ 30 126 183 173 40 113 10]3.4 NumPy Data Types
NumPy provides additional and more specific data types in comparison to base Python. Here, we provide a brief explanation of commonly used data types.
-
bool_: BooleanTrueorFalse -
int8: 8-bit signed integer (-128 to 127) -
int16: 16-bit signed integer (-32,768 to 32,767) -
int32: 32-bit signed integer (-2,147,483,648 to 2,147,483,647) -
int64: 64-bit signed integer (-9,223,372,036,854,775,808 to 9,223,372,036,854,775,807) -
uint8: 8-bit unsigned integer (0 to 255) -
uint16: 16-bit unsigned integer (0 to 65,535) -
uint32: 32-bit unsigned integer (0 to 4,294,967,295) -
uint64: 64-bit unsigned integer (0 to 18,446,744,073,709,551,615) -
float16: half precision float -
float32: single precision float -
float64: double precision float
Signed integers can differentiate positive and negative values while unsigned integers cannot. Float data can store decimal values while integer data cannot. There are also data types for complex numbers, which we do not discuss here.
Below, we demonstrate how to define the data type with the dtype parameter. In all cases, we use ones() to create an array with three elements. For both int8 and int16, 1 as an integer value is returned. When the data are defined as float16, 1 as a float value is returned (1.). Lastly, when the type is set to bool_, Boolean True is returned since 1 indicates True and 0 indicates False.
The data type impacts the amount of memory needed. For example, an int8 requires less memory than an int16. We generally try to use the data type that can represent the data with the least amount of memory unless a specific data type is needed in an analysis.
arr1 = np.ones((3), dtype="int8")
print(arr1)[1 1 1]print(arr1.dtype)int8arr1 = np.ones((3), dtype="int16")
print(arr1)[1 1 1]print(arr1.dtype)int16arr1 = np.ones((3), dtype="float16")
print(arr1)[1. 1. 1.]print(arr1.dtype)float16arr1 = np.ones((3), dtype="bool_")
print(arr1)[ True True True]print(arr1.dtype)bool3.5 Understanding and Manipulating Array Shape and Dimensions
Let’s spend some time discussing the dimensions and shape of an array. The shape of an array relates to the length of each dimension. The len() function returns the length of the first dimension (in this case 3). To obtain a tuple of the lengths for all dimensions, you must use the shape property. So, the array generated has three dimensions with lengths of 3, 4, and 4, respectively. The number of dimensions is returned with the ndim property. The size property returns the number of elements in the array. There are 48 elements in the example array: 3 X 4 X 4 = 48. The dtype property provides the data type.
arr6 = np.ones((3, 4, 4))
print("Length of first dimension: " + str(len(arr6)))Length of first dimension: 3print("Shape of array: " + str(arr6.shape))Shape of array: (3, 4, 4)print("Number of dimensions: " + str(arr6.ndim))Number of dimensions: 3print("Size of array: " + str(arr6.size))Size of array: 48print("Data type of array: " + str(arr6.dtype))Data type of array: float64NumPy has a built-in method for changing the shape of an array: reshape(). Note that the number of elements, or size, of the array must perfectly fill the new shape. In the first example, we maintain the number of dimensions but change the shape or length of each dimension. In the next two examples, we convert the three-dimensional array to two-dimensional arrays. Lastly, we convert the array to a one-dimensional array, or vector, with a length of 48.
arr6b = arr6.reshape(4, 4, 3)
arr6c = arr6.reshape(4, 12)
arr6d = arr6.reshape(12, 4)
arr6e = arr6.reshape(48)
print(arr6b)[[[1. 1. 1.]
[1. 1. 1.]
[1. 1. 1.]
[1. 1. 1.]]
[[1. 1. 1.]
[1. 1. 1.]
[1. 1. 1.]
[1. 1. 1.]]
[[1. 1. 1.]
[1. 1. 1.]
[1. 1. 1.]
[1. 1. 1.]]
[[1. 1. 1.]
[1. 1. 1.]
[1. 1. 1.]
[1. 1. 1.]]]print(arr6c)[[1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
[1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
[1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
[1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]]print(arr6d)[[1. 1. 1. 1.]
[1. 1. 1. 1.]
[1. 1. 1. 1.]
[1. 1. 1. 1.]
[1. 1. 1. 1.]
[1. 1. 1. 1.]
[1. 1. 1. 1.]
[1. 1. 1. 1.]
[1. 1. 1. 1.]
[1. 1. 1. 1.]
[1. 1. 1. 1.]
[1. 1. 1. 1.]]print(arr6e)[1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]As a second example, here we convert a vector into a matrix.
arr10 = np.random.randint(1, 1200, 100)
arr10b = arr10.reshape(10, 10)
print(arr10b)[[ 481 685 889 158 275 638 408 60 742 1055]
[ 94 83 458 451 185 459 989 903 527 782]
[ 707 624 313 256 575 7 108 301 421 829]
[ 331 525 689 695 822 129 71 696 912 397]
[ 594 393 313 246 842 905 606 214 567 222]
[ 79 28 258 287 885 568 385 771 485 867]
[ 920 63 96 1016 611 54 1162 415 23 32]
[1006 226 1183 1163 258 664 509 913 319 570]
[ 206 1134 1101 734 1104 161 1181 243 1149 772]
[ 126 1071 428 598 820 1007 474 1119 593 203]]When reshaping, it is possible to have NumPy determine the appropriate size of a single dimension to fill an array with the available elements. This is accomplished using -1 in the array dimension location when applying the reshape() method.
arr10 = np.random.randint(1, 1200, 1000)
arr10b = arr10.reshape(-1, 10, 10)
arr10c = arr10.reshape(10, -1, 10)
arr10d = arr10.reshape(10, 10, -1)
print(arr10.shape)(1000,)print(arr10b.shape)(10, 10, 10)print(arr10c.shape)(10, 10, 10)print(arr10d.shape)(10, 10, 10)3.6 NumPy Array Indexing
Similar to lists, NumPy arrays are indexed. So, values from the array can be extracted or referenced using their associated index. Since arrays often have multiple dimensions, indexes must also extend into multiple dimensions. See the comments below for general array indexing rules. Remember that indexing starts at 0, the first index provided is included, and the last index provided is not included. Extracting portions of an array is known as slicing.
arr11 = np.linspace(0, 50, 50, dtype="int")
arr12 = arr11.reshape(2,5,5)
print("Original array")Original arrayprint(arr12)[[[ 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]
[40 41 42 43 44]
[45 46 47 48 50]]]print("All values in first index of first dimension")All values in first index of first dimensionprint(arr12[0]) #This will extract just the values from the first index in the first dimension[[ 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]]print("All values in second index of first dimension")All values in second index of first dimensionprint(arr12[1]) #This will extract just the values from the second index in the first dimension[[25 26 27 28 29]
[30 31 32 33 34]
[35 36 37 38 39]
[40 41 42 43 44]
[45 46 47 48 50]]print("All values in first index of first dimension and first index of second dimension")All values in first index of first dimension and first index of second dimensionprint(arr12[0][0]) #This will extract all values occurring in the first index of both the first and second dimensions[0 1 2 3 4]print("A single value specified with three indexes, one for each dimension")A single value specified with three indexes, one for each dimensionprint(arr12[1, 3, 3]) #This will extract a specific value based on an index in all three dimensions43print("Incorporating ranges")Incorporating rangesprint(arr12[1, 0:2, 0:2]) #All values in second index of first dimension that are also included in the first to second index of the second and third dimensions[[25 26]
[30 31]]print("Using colons")Using colonsprint(arr12[:, 0:2, 0:2]) #Only a colon means select all values in a dimension[[[ 0 1]
[ 5 6]]
[[25 26]
[30 31]]]print(arr12[:,2:,0:2]) #Can also use colons to select all values before or after an index[[[10 11]
[15 16]
[20 21]]
[[35 36]
[40 41]
[45 46]]]Once values have been selected using index notation, they can be changed. In the example below, we convert all values in the first index of the first dimension and the first index of the second dimension to 0.
arr12[0][0] = 0
print(arr12)[[[ 0 0 0 0 0]
[ 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]
[40 41 42 43 44]
[45 46 47 48 50]]]3.7 Boolean Arrays
It is also possible to create arrays of Boolean values, as demonstrated below.
arr13 = np.array([True, False, True, False, True, False, True, False, False])
arr13b = arr13.reshape(3, 3)
print(arr13b)[[ True False True]
[False True False]
[ True False False]]Comparison operators can be used to compare each value in an array to a value and return the Boolean result to the associated position in a new array.
arr10 = np.random.randint(1, 1200, 100)
arr10b = arr10.reshape(10, 10)
print(arr10b)[[ 372 1135 607 180 515 197 188 811 962 546]
[ 543 809 933 871 1136 317 488 939 553 1070]
[1179 371 309 829 932 319 1004 930 1062 769]
[ 145 497 616 440 733 622 1039 96 344 728]
[ 222 1092 346 912 140 923 904 807 1102 497]
[ 280 333 1190 326 153 14 1163 65 506 957]
[1080 915 1185 313 666 458 212 1049 32 910]
[ 7 801 595 1003 359 1164 556 83 62 1083]
[ 174 416 923 746 144 933 716 1156 919 1124]
[1159 651 369 1047 754 263 47 650 750 96]]arr10bool = arr10b > 150
print(arr10bool)[[ True True True True True True True True True True]
[ True True True True True True True True True True]
[ True True True True True True True True True True]
[False True True True True True True False True True]
[ True True True True False True True True True True]
[ True True True True True False True False True True]
[ True True True True True True True True False True]
[False True True True True True True False False True]
[ True True True True False True True True True True]
[ True True True True True True False True True False]]3.8 Copy vs. View
In Chapter 1, we explained how, for mutable data types, setting a variable equal to another variable results in a reference to the original object or data in memory. So, changes to the original object or the new object will change both since they reference the same data.
The behavior for NumPy arrays is similar. Two methods are available for creating a new variable relative to an existing variable: copy() and view().
When using view(), the variable will reference the same data or object in memory. So, changes to the original variable or the new variable created using view() will result in changes to both. Using view(), you can also reference portions of an array. This allows you to work with a subset of the data values without copying or replicating the data in memory.
In contrast to view(), copy() will copy the data or object in memory, so changes made to either the original or the copied object will not impact the other.
The three examples below demonstrate this behavior. In the first example, arr2 is created as a view of arr1 while arr3 is created as a copy of arr1. A subsequent change to arr1 changes arr1 and arr2 but not arr3. In the second example, a change to arr2, a view of arr1, impacts both arr1 and arr2 but not arr3. Lastly, changes to arr3 impact only arr3 and not arr1 or arr2 since it is a copy of arr1 as opposed to a view.
In summary, if you want to make a copy of an array as opposed to referencing the original data, you should use the copy() method.
import copy
arr1 = np.array(25)
arr2 = np.array(25)
arr3 = np.array(25)
arr1 = np.random.randint(1, 100, 25)
arr1 = arr1.reshape(5, 5)
print(arr1)[[79 32 32 34 49]
[70 79 93 37 8]
[19 18 71 75 87]
[34 61 61 75 91]
[85 19 23 10 18]]arr2 = arr1.view()
arr3 = arr1.copy()
arr1[:, :] = 0
print(arr1)[[0 0 0 0 0]
[0 0 0 0 0]
[0 0 0 0 0]
[0 0 0 0 0]
[0 0 0 0 0]]print(arr2)[[0 0 0 0 0]
[0 0 0 0 0]
[0 0 0 0 0]
[0 0 0 0 0]
[0 0 0 0 0]]print(arr3)[[79 32 32 34 49]
[70 79 93 37 8]
[19 18 71 75 87]
[34 61 61 75 91]
[85 19 23 10 18]]import copy
arr1 = np.array(25)
arr2 = np.array(25)
arr3 = np.array(25)
arr1 = np.random.randint(1, 100, 25)
arr1 = arr1.reshape(5, 5)
print(arr1)[[86 2 78 30 99]
[93 32 61 91 23]
[84 81 60 9 89]
[10 95 11 21 84]
[26 19 23 4 80]]arr2 = arr1.view()
arr3 = arr1.copy()
arr2[:, :] = 0
print(arr1)[[0 0 0 0 0]
[0 0 0 0 0]
[0 0 0 0 0]
[0 0 0 0 0]
[0 0 0 0 0]]print(arr2)[[0 0 0 0 0]
[0 0 0 0 0]
[0 0 0 0 0]
[0 0 0 0 0]
[0 0 0 0 0]]print(arr3)[[86 2 78 30 99]
[93 32 61 91 23]
[84 81 60 9 89]
[10 95 11 21 84]
[26 19 23 4 80]]import copy
arr1 = np.array(25)
arr2 = np.array(25)
arr3 = np.array(25)
arr1 = np.random.randint(1, 100, 25)
arr1 = arr1.reshape(5, 5)
print(arr1)[[51 15 46 46 67]
[17 11 14 2 52]
[86 55 50 74 5]
[48 34 17 52 41]
[83 26 16 11 32]]arr2 = arr1.view()
arr3 = arr1.copy()
arr3[:, :] = 0
print(arr1)[[51 15 46 46 67]
[17 11 14 2 52]
[86 55 50 74 5]
[48 34 17 52 41]
[83 26 16 11 32]]print(arr2)[[51 15 46 46 67]
[17 11 14 2 52]
[86 55 50 74 5]
[48 34 17 52 41]
[83 26 16 11 32]]print(arr3)[[0 0 0 0 0]
[0 0 0 0 0]
[0 0 0 0 0]
[0 0 0 0 0]
[0 0 0 0 0]]3.9 Array Arithmetic and Operations
It is generally easy to perform mathematical operations on arrays, as demonstrated below. In all cases, the same operation is applied to all elements in the array.
arr14 = np.random.randint(1, 1200, 25)
arr14b = arr14.reshape(5, 5)
print(arr14b)[[ 593 526 654 1197 1037]
[ 353 1010 1160 823 871]
[ 301 720 1163 640 40]
[1099 407 450 862 260]
[1195 403 429 829 926]]print(arr14b+21)[[ 614 547 675 1218 1058]
[ 374 1031 1181 844 892]
[ 322 741 1184 661 61]
[1120 428 471 883 281]
[1216 424 450 850 947]]print(arr14b-52)[[ 541 474 602 1145 985]
[ 301 958 1108 771 819]
[ 249 668 1111 588 -12]
[1047 355 398 810 208]
[1143 351 377 777 874]]print(arr14b*2)[[1186 1052 1308 2394 2074]
[ 706 2020 2320 1646 1742]
[ 602 1440 2326 1280 80]
[2198 814 900 1724 520]
[2390 806 858 1658 1852]]print(arr14b/3)[[197.66666667 175.33333333 218. 399. 345.66666667]
[117.66666667 336.66666667 386.66666667 274.33333333 290.33333333]
[100.33333333 240. 387.66666667 213.33333333 13.33333333]
[366.33333333 135.66666667 150. 287.33333333 86.66666667]
[398.33333333 134.33333333 143. 276.33333333 308.66666667]]print(arr14b**2)[[ 351649 276676 427716 1432809 1075369]
[ 124609 1020100 1345600 677329 758641]
[ 90601 518400 1352569 409600 1600]
[1207801 165649 202500 743044 67600]
[1428025 162409 184041 687241 857476]]It is also possible to perform mathematical operations on sets of arrays as long as they have the same shape. In such cases, elements are matched based on having the same position within the array.
arr14 = np.random.randint(1, 1200, 25)
arr14b = arr14.reshape(5, 5)
print(arr14b)[[ 398 1019 732 162 1115]
[ 261 838 200 916 224]
[ 747 583 794 442 335]
[ 61 781 127 482 195]
[ 47 288 499 116 1162]]print(arr14b+arr14b)[[ 796 2038 1464 324 2230]
[ 522 1676 400 1832 448]
[1494 1166 1588 884 670]
[ 122 1562 254 964 390]
[ 94 576 998 232 2324]]print(arr14b-arr14b)[[0 0 0 0 0]
[0 0 0 0 0]
[0 0 0 0 0]
[0 0 0 0 0]
[0 0 0 0 0]]To summarize the results from above, it is possible to:
- perform mathematical operations between an array with any shape and a scalar (i.e., a single value).
- perform mathematical operations between arrays that have the same shape.
There are some other cases in which it is possible to perform mathematical operations using a technique known as broadcasting. The following rules summarize when broadcasting can be used and how.
- If two arrays have a different number of dimensions, the shape of the array with fewer dimensions is padded with ones on its leading side (for example, to multiply an array of shape (3) by an array of shape (3,3), the first array is treated as having shape (1,3)).
- If the shape of the arrays does not match in a dimension, the array with shape equal to 1 in that dimension is stretched to match the other shape.
- If in any dimension the sizes disagree and neither is equal to 1, an error is raised.
In the example below, an array of shape (6, 6) is multiplied by an array of shape (6). This requires that the second array be broadcast to a shape of (1, 6).
arr1 = np.random.randint(1, 100, 36)
arr1b = arr1.reshape(6, 6)
arr2 = np.ones((6))
arr2[:] = 2
print(arr1b)[[75 31 9 93 27 63]
[86 23 34 41 44 33]
[97 58 60 86 74 49]
[87 39 18 8 59 49]
[56 63 85 95 45 35]
[62 79 70 84 80 24]]print(arr2)[2. 2. 2. 2. 2. 2.]print(arr1b*arr2)[[150. 62. 18. 186. 54. 126.]
[172. 46. 68. 82. 88. 66.]
[194. 116. 120. 172. 148. 98.]
[174. 78. 36. 16. 118. 98.]
[112. 126. 170. 190. 90. 70.]
[124. 158. 140. 168. 160. 48.]]NumPy provides mathematical functions and methods for performing common tasks. The block of code below provides some examples.
arr14 = np.random.randint(1, 1200, 25)
arr14b = arr14.reshape(5, 5)
print(np.max(arr14b))1148print(np.min(arr14b))14print(np.sqrt(arr14b))[[20.4450483 22.97825059 23.04343724 24.85960579 23.60084744]
[10.86278049 30.33150178 30.46309242 28.82707061 3.74165739]
[16.0623784 18.35755975 33.88214869 25.84569597 32.24903099]
[ 8.48528137 9.05538514 33.86738844 17.29161647 14.03566885]
[15.45962483 30.46309242 15.42724862 25.11971337 13.7113092 ]]3.10 Xarray
3.10.1 How Xarray Differs from NumPy
NumPy arrays are fast and flexible, but they only keep track of the position of each value. A NumPy array does not know what its dimensions represent. When working with geospatial and environmental data, however, dimensions almost always have meaning. A multispectral image has band, row (y), and column (x) dimensions, and a stack of monthly climate grids has time, latitude, and longitude dimensions. With NumPy, it is up to you to remember that axis 0 is time, axis 1 is y, and axis 2 is x. It is also up to you to separately store the dates, coordinates, and units associated with the data. This makes it easy to make mistakes, such as calculating a mean over the wrong axis.
The Xarray library addresses these issues by adding labels to NumPy arrays. Here are the key differences between Xarray and NumPy:
-
Named dimensions: Each dimension is given a name, such as “time”, “y”, or “x”. Operations can then reference dimensions by name (for example,
mean(dim="time")) as opposed to by axis number (mean(axis=0)). - Coordinates: Each position along a dimension can be assigned a label, such as a date, a map coordinate, or a band wavelength. Data can then be selected using these labels as opposed to only using indexes.
- Attributes: Arbitrary metadata, such as units or a description, can be stored with the data.
- Broadcasting and alignment by name: When performing operations on multiple arrays, Xarray matches dimensions based on their names as opposed to their order and shape.
- Datasets: Multiple related variables that share dimensions and coordinates can be stored together in a single object, similar to how a NetCDF file stores multiple variables.
- File input/output: Xarray can read and write common scientific data formats, such as NetCDF and Zarr, and, with the rioxarray extension, geospatial raster formats such as GeoTIFF. It can also work with datasets that are too large to fit in memory using Dask.
Xarray does not replace NumPy. The data in an Xarray object are stored as a NumPy array, and most NumPy functions can be applied to Xarray objects directly. In general, we use NumPy for lower-level array manipulation and when preparing data for other libraries that expect NumPy arrays, such as scikit-learn. We use Xarray when working with labeled, multidimensional data, such as time series of raster grids or climate model output.
Xarray provides two primary data structures:
- DataArray: a single labeled, multidimensional array (the labeled equivalent of a NumPy array)
- Dataset: a collection of DataArrays that share dimensions and coordinates
It is common to assign Xarray an alias of xr.
import xarray as xr3.10.2 Creating a DataArray
To compare NumPy and Xarray, we recreate the three-dimensional array that was used above to demonstrate indexing. Imagine that this array represents two monthly grids of mean temperature in degrees Celsius, one for June and one for July, each with 5 rows and 5 columns of cells.
A DataArray is created using the DataArray() function. In the example, we provide the NumPy array, a tuple of dimension names, a dictionary of coordinates for each dimension, a name for the variable, and a dictionary of attributes. The y coordinates decrease from the top to the bottom of the grid, which is the convention for raster data. Note that the printed DataArray includes the dimension names and lengths, the data values, the coordinates, and the attributes.
arr12 = np.linspace(0, 50, 50, dtype="int").reshape(2, 5, 5)
da12 = xr.DataArray(
arr12,
dims=("time", "y", "x"),
coords={
"time": np.array(["2024-06-01", "2024-07-01"], dtype="datetime64[ns]"),
"y": [500, 400, 300, 200, 100],
"x": [10, 20, 30, 40, 50]
},
name="temp",
attrs={"units": "degrees C", "description": "Monthly mean temperature"}
)
print(da12)<xarray.DataArray 'temp' (time: 2, y: 5, x: 5)> Size: 400B
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],
[40, 41, 42, 43, 44],
[45, 46, 47, 48, 50]]])
Coordinates:
* time (time) datetime64[ns] 16B 2024-06-01 2024-07-01
* y (y) int64 40B 500 400 300 200 100
* x (x) int64 40B 10 20 30 40 50
Attributes:
units: degrees C
description: Monthly mean temperatureMany of the properties available for NumPy arrays are also available for DataArrays. In addition, the dims property returns the dimension names, and the sizes property returns the length of each dimension associated with its name. The values property returns the underlying data as a NumPy array.
print("Dimension names: " + str(da12.dims))Dimension names: ('time', 'y', 'x')print("Shape of array: " + str(da12.shape))Shape of array: (2, 5, 5)print("Dimension sizes: " + str(dict(da12.sizes)))Dimension sizes: {'time': 2, 'y': 5, 'x': 5}print("Number of dimensions: " + str(da12.ndim))Number of dimensions: 3print("Size of array: " + str(da12.size))Size of array: 50print("Data type of array: " + str(da12.dtype))Data type of array: int64print(type(da12.values))<class 'numpy.ndarray'>Coordinates and attributes can be accessed using the coords and attrs properties or by using the coordinate name.
print(da12.coords)Coordinates:
* time (time) datetime64[ns] 16B 2024-06-01 2024-07-01
* y (y) int64 40B 500 400 300 200 100
* x (x) int64 40B 10 20 30 40 50print(da12.attrs){'units': 'degrees C', 'description': 'Monthly mean temperature'}print(da12["time"].values)['2024-06-01T00:00:00.000000000' '2024-07-01T00:00:00.000000000']3.10.3 Indexing and Selecting
Xarray provides two methods for selecting data. The isel() method selects using integer indexes, similar to NumPy, while the sel() method selects using coordinate labels. In both cases, the dimension is referenced by name, so the order in which the dimensions are listed does not matter.
In the example below, the same value is extracted in three ways: (1) using NumPy indexing, (2) using isel() and indexes, and (3) using sel() and coordinate labels. Note that the last two dimensions are provided in reverse order in the sel() example, which is not a problem since the dimensions are referenced by name. The item() method returns the single value as a standard Python number.
print(arr12[1, 3, 3])43print(da12.isel(time=1, y=3, x=3).item())43print(da12.sel(time="2024-07-01", x=40, y=200).item())43Ranges are defined using slice(). When using isel(), slices behave like NumPy: the first index is included and the last index is not. In contrast, when using sel(), both the start and end labels are included. In the example, the first selection is equivalent to arr12[:, 0:2, 0:2], and the second selection returns all cells with x coordinates between 10 and 20, inclusive. Note that the coordinates are subset along with the data.
print(da12.isel(y=slice(0, 2), x=slice(0, 2)))<xarray.DataArray 'temp' (time: 2, y: 2, x: 2)> Size: 64B
array([[[ 0, 1],
[ 5, 6]],
[[25, 26],
[30, 31]]])
Coordinates:
* time (time) datetime64[ns] 16B 2024-06-01 2024-07-01
* y (y) int64 16B 500 400
* x (x) int64 16B 10 20
Attributes:
units: degrees C
description: Monthly mean temperatureprint(da12.sel(x=slice(10, 20)))<xarray.DataArray 'temp' (time: 2, y: 5, x: 2)> Size: 160B
array([[[ 0, 1],
[ 5, 6],
[10, 11],
[15, 16],
[20, 21]],
[[25, 26],
[30, 31],
[35, 36],
[40, 41],
[45, 46]]])
Coordinates:
* time (time) datetime64[ns] 16B 2024-06-01 2024-07-01
* y (y) int64 40B 500 400 300 200 100
* x (x) int64 16B 10 20
Attributes:
units: degrees C
description: Monthly mean temperatureSince time is stored as a date, it is also possible to select using partial dates. For example, "2024-07" will select all data from July 2024.
print(da12.sel(time="2024-07"))<xarray.DataArray 'temp' (time: 1, y: 5, x: 5)> Size: 200B
array([[[25, 26, 27, 28, 29],
[30, 31, 32, 33, 34],
[35, 36, 37, 38, 39],
[40, 41, 42, 43, 44],
[45, 46, 47, 48, 50]]])
Coordinates:
* time (time) datetime64[ns] 8B 2024-07-01
* y (y) int64 40B 500 400 300 200 100
* x (x) int64 40B 10 20 30 40 50
Attributes:
units: degrees C
description: Monthly mean temperature3.10.4 Arithmetic and Summary Statistics
Mathematical operations work the same as with NumPy arrays. However, the result is a DataArray that maintains the dimension names and coordinates. NumPy functions, such as np.sqrt(), can also be applied directly to DataArrays. In the example, we convert the temperature values from degrees Celsius to degrees Fahrenheit. Since the units have changed, we update the attributes using assign_attrs().
da12_f = (da12 * 9/5 + 32).assign_attrs(units="degrees F")
print(da12_f)<xarray.DataArray 'temp' (time: 2, y: 5, x: 5)> Size: 400B
array([[[ 32. , 33.8, 35.6, 37.4, 39.2],
[ 41. , 42.8, 44.6, 46.4, 48.2],
[ 50. , 51.8, 53.6, 55.4, 57.2],
[ 59. , 60.8, 62.6, 64.4, 66.2],
[ 68. , 69.8, 71.6, 73.4, 75.2]],
[[ 77. , 78.8, 80.6, 82.4, 84.2],
[ 86. , 87.8, 89.6, 91.4, 93.2],
[ 95. , 96.8, 98.6, 100.4, 102.2],
[104. , 105.8, 107.6, 109.4, 111.2],
[113. , 114.8, 116.6, 118.4, 122. ]]])
Coordinates:
* time (time) datetime64[ns] 16B 2024-06-01 2024-07-01
* y (y) int64 40B 500 400 300 200 100
* x (x) int64 40B 10 20 30 40 50
Attributes:
units: degrees F
description: Monthly mean temperatureprint(np.sqrt(da12))<xarray.DataArray 'temp' (time: 2, y: 5, x: 5)> Size: 400B
array([[[0. , 1. , 1.41421356, 1.73205081, 2. ],
[2.23606798, 2.44948974, 2.64575131, 2.82842712, 3. ],
[3.16227766, 3.31662479, 3.46410162, 3.60555128, 3.74165739],
[3.87298335, 4. , 4.12310563, 4.24264069, 4.35889894],
[4.47213595, 4.58257569, 4.69041576, 4.79583152, 4.89897949]],
[[5. , 5.09901951, 5.19615242, 5.29150262, 5.38516481],
[5.47722558, 5.56776436, 5.65685425, 5.74456265, 5.83095189],
[5.91607978, 6. , 6.08276253, 6.164414 , 6.244998 ],
[6.32455532, 6.40312424, 6.4807407 , 6.55743852, 6.63324958],
[6.70820393, 6.78232998, 6.8556546 , 6.92820323, 7.07106781]]])
Coordinates:
* time (time) datetime64[ns] 16B 2024-06-01 2024-07-01
* y (y) int64 40B 500 400 300 200 100
* x (x) int64 40B 10 20 30 40 50
Attributes:
units: degrees C
description: Monthly mean temperatureWhere Xarray really differs from NumPy is in summarizing data along dimensions. With NumPy, you must provide the axis number. With Xarray, you provide the dimension name. In the first example below, both lines calculate the mean temperature in each cell across the two months. The Xarray version is easier to interpret, and the coordinates for the remaining dimensions are maintained. In the second example, we calculate the maximum temperature in each month by summarizing over both spatial dimensions.
print(arr12.mean(axis=0))[[12.5 13.5 14.5 15.5 16.5]
[17.5 18.5 19.5 20.5 21.5]
[22.5 23.5 24.5 25.5 26.5]
[27.5 28.5 29.5 30.5 31.5]
[32.5 33.5 34.5 35.5 37. ]]print(da12.mean(dim="time"))<xarray.DataArray 'temp' (y: 5, x: 5)> Size: 200B
array([[12.5, 13.5, 14.5, 15.5, 16.5],
[17.5, 18.5, 19.5, 20.5, 21.5],
[22.5, 23.5, 24.5, 25.5, 26.5],
[27.5, 28.5, 29.5, 30.5, 31.5],
[32.5, 33.5, 34.5, 35.5, 37. ]])
Coordinates:
* y (y) int64 40B 500 400 300 200 100
* x (x) int64 40B 10 20 30 40 50
Attributes:
units: degrees C
description: Monthly mean temperatureprint(da12.max(dim=["y", "x"]))<xarray.DataArray 'temp' (time: 2)> Size: 16B
array([24, 50])
Coordinates:
* time (time) datetime64[ns] 16B 2024-06-01 2024-07-01
Attributes:
units: degrees C
description: Monthly mean temperature3.10.5 Comparisons and Masking
Comparison operators work the same as in NumPy and return a Boolean DataArray. To demonstrate this, we reuse the random array of integers from the Boolean Arrays section above. The where() method can then be used to mask values: values that meet the condition are maintained while all other values are replaced with NaN. Note that the result is converted to the float data type since integers cannot store NaN. A second argument can be provided to where() to replace values that do not meet the condition with a specific value as opposed to NaN.
arr10 = np.random.randint(1, 1200, 100)
arr10b = arr10.reshape(10, 10)
da10 = xr.DataArray(arr10b, dims=("y", "x"))
print(da10 > 150)<xarray.DataArray (y: 10, x: 10)> Size: 100B
array([[ True, True, True, True, True, True, False, True, True,
True],
[ True, True, True, True, True, True, True, False, True,
True],
[ True, True, True, True, True, True, True, True, True,
True],
[ True, True, True, False, True, True, False, True, True,
True],
[ True, True, True, True, True, False, True, True, True,
True],
[ True, True, True, True, False, True, True, True, True,
True],
[ True, True, True, True, True, True, True, True, True,
True],
[ True, True, True, True, True, True, True, True, True,
True],
[ True, True, True, True, True, True, True, True, True,
True],
[ True, True, True, True, True, True, True, True, True,
True]])
Dimensions without coordinates: y, xprint(da10.where(da10 > 150))<xarray.DataArray (y: 10, x: 10)> Size: 800B
array([[ 490., 438., 1191., 278., 228., 300., nan, 390., 711.,
1035.],
[ 883., 519., 849., 927., 901., 164., 462., nan, 515.,
484.],
[ 435., 382., 1070., 267., 677., 256., 313., 821., 811.,
463.],
[ 969., 1159., 813., nan, 344., 763., nan, 1193., 637.,
947.],
[ 913., 1034., 555., 954., 1138., nan, 246., 429., 818.,
698.],
[ 592., 350., 876., 607., nan, 720., 1120., 664., 583.,
361.],
[ 235., 920., 1141., 1104., 1072., 409., 1083., 835., 341.,
576.],
[ 711., 720., 361., 895., 1097., 903., 642., 362., 359.,
614.],
[ 324., 267., 195., 1077., 219., 1168., 924., 973., 407.,
199.],
[1178., 282., 185., 665., 730., 1186., 240., 1180., 295.,
848.]])
Dimensions without coordinates: y, xprint(da10.where(da10 > 150, 0))<xarray.DataArray (y: 10, x: 10)> Size: 400B
array([[ 490, 438, 1191, 278, 228, 300, 0, 390, 711, 1035],
[ 883, 519, 849, 927, 901, 164, 462, 0, 515, 484],
[ 435, 382, 1070, 267, 677, 256, 313, 821, 811, 463],
[ 969, 1159, 813, 0, 344, 763, 0, 1193, 637, 947],
[ 913, 1034, 555, 954, 1138, 0, 246, 429, 818, 698],
[ 592, 350, 876, 607, 0, 720, 1120, 664, 583, 361],
[ 235, 920, 1141, 1104, 1072, 409, 1083, 835, 341, 576],
[ 711, 720, 361, 895, 1097, 903, 642, 362, 359, 614],
[ 324, 267, 195, 1077, 219, 1168, 924, 973, 407, 199],
[1178, 282, 185, 665, 730, 1186, 240, 1180, 295, 848]],
dtype=int32)
Dimensions without coordinates: y, x3.10.6 Broadcasting by Dimension Name
Above, we multiplied an array of shape (6, 6) by an array of shape (6). Following the NumPy broadcasting rules, the second array was treated as a single row and applied to every row (in other words, each column was multiplied by a different value). What if you instead wanted to multiply each row by a different value? With NumPy, you would need to reshape the second array to (6, 1). With Xarray, you simply assign the second array to the appropriate dimension name.
In the example, we use the values 1 through 6 as opposed to 2 so that the result is easier to interpret. When the multiplier is assigned to the “x” dimension, each column is multiplied by a different value. When it is assigned to the “y” dimension, each row is multiplied by a different value. The NumPy equivalent is provided for comparison.
arr1 = np.random.randint(1, 100, 36)
arr1b = arr1.reshape(6, 6)
mult = np.array([1, 2, 3, 4, 5, 6])
da1 = xr.DataArray(arr1b, dims=("y", "x"))
mult_x = xr.DataArray(mult, dims="x")
mult_y = xr.DataArray(mult, dims="y")
print(da1.values)[[49 20 97 18 20 1]
[ 4 70 58 50 91 53]
[15 98 19 33 54 68]
[89 13 37 5 22 58]
[98 95 32 38 43 47]
[29 87 94 12 82 22]]print((da1 * mult_x).values)[[ 49 40 291 72 100 6]
[ 4 140 174 200 455 318]
[ 15 196 57 132 270 408]
[ 89 26 111 20 110 348]
[ 98 190 96 152 215 282]
[ 29 174 282 48 410 132]]print((da1 * mult_y).values)[[ 49 20 97 18 20 1]
[ 8 140 116 100 182 106]
[ 45 294 57 99 162 204]
[356 52 148 20 88 232]
[490 475 160 190 215 235]
[174 522 564 72 492 132]]print(np.array_equal(da1 * mult_y, arr1b * mult.reshape(6, 1)))TrueIf two arrays have different dimension names, Xarray combines them into an array with all of the dimensions. For example, multiplying a DataArray with a “y” dimension by a DataArray with an “x” dimension results in a two-dimensional array.
print((mult_y * mult_x).dims)('y', 'x')print((mult_y * mult_x).values)[[ 1 2 3 4 5 6]
[ 2 4 6 8 10 12]
[ 3 6 9 12 15 18]
[ 4 8 12 16 20 24]
[ 5 10 15 20 25 30]
[ 6 12 18 24 30 36]]3.10.7 Datasets
A Dataset stores multiple DataArrays that share dimensions and coordinates. In the example, we combine the temperature grids in degrees Celsius and degrees Fahrenheit into a single Dataset using a dictionary in which the keys are the variable names. Variables can then be accessed using bracket or dot notation, and selections and summaries can be applied to all variables at once.
ds = xr.Dataset({"temp_c": da12, "temp_f": da12_f})
print(ds)<xarray.Dataset> Size: 896B
Dimensions: (time: 2, y: 5, x: 5)
Coordinates:
* time (time) datetime64[ns] 16B 2024-06-01 2024-07-01
* y (y) int64 40B 500 400 300 200 100
* x (x) int64 40B 10 20 30 40 50
Data variables:
temp_c (time, y, x) int64 400B 0 1 2 3 4 5 6 7 ... 42 43 44 45 46 47 48 50
temp_f (time, y, x) float64 400B 32.0 33.8 35.6 37.4 ... 116.6 118.4 122.0print(ds["temp_f"].sel(time="2024-07-01"))<xarray.DataArray 'temp_f' (y: 5, x: 5)> Size: 200B
array([[ 77. , 78.8, 80.6, 82.4, 84.2],
[ 86. , 87.8, 89.6, 91.4, 93.2],
[ 95. , 96.8, 98.6, 100.4, 102.2],
[104. , 105.8, 107.6, 109.4, 111.2],
[113. , 114.8, 116.6, 118.4, 122. ]])
Coordinates:
* y (y) int64 40B 500 400 300 200 100
* x (x) int64 40B 10 20 30 40 50
time datetime64[ns] 8B 2024-07-01
Attributes:
units: degrees F
description: Monthly mean temperatureprint(ds.mean(dim=["y", "x"]))<xarray.Dataset> Size: 48B
Dimensions: (time: 2)
Coordinates:
* time (time) datetime64[ns] 16B 2024-06-01 2024-07-01
Data variables:
temp_c (time) float64 16B 12.0 37.04
temp_f (time) float64 16B 53.6 98.673.10.8 Converting Between NumPy, Xarray, and Pandas
Since many libraries expect NumPy arrays, it is common to convert between the two. The values property or to_numpy() method returns the data as a NumPy array, and a NumPy array can be converted to a DataArray using DataArray(), as demonstrated above. A DataArray can also be converted to a Pandas DataFrame, which is discussed in Chapter 4, using to_dataframe(). Each cell becomes a row, and the coordinates become the index.
Similar to NumPy, setting a variable equal to a DataArray will reference the same object in memory. Use the copy() method to create an independent copy.
arr12_back = da12.to_numpy()
print(type(arr12_back))<class 'numpy.ndarray'>print(arr12_back.shape)(2, 5, 5)print(da12.to_dataframe().head(10)) temp
time y x
2024-06-01 500 10 0
20 1
30 2
40 3
50 4
400 10 5
20 6
30 7
40 8
50 9
da12_copy = da12.copy()3.10.9 Xarray and Geospatial Data
Xarray is commonly used to work with geospatial raster data and Earth science data. The open_dataset() function reads NetCDF and other supported files as a Dataset. The rioxarray extension, which is not included in the pybook environment, adds a rio accessor for reading and writing geospatial raster formats, such as GeoTIFF, and for working with coordinate reference systems. A multiband image read using rioxarray is returned as a DataArray with “band”, “y”, and “x” dimensions and the map coordinates of each cell.
The code below is provided as an example only and is not executed.
import rioxarray
img = rioxarray.open_rasterio("datasets/my_image.tif")
print(img.dims)
print(img.rio.crs)
red = img.sel(band=3)
nir = img.sel(band=4)
ndvi = (nir - red)/(nir + red)3.11 Concluding Remarks
As mentioned above, NumPy is central to using Python for analyzing data. So, an understanding of NumPy is important for data scientists and geospatial data scientists. Additional libraries make use of NumPy to expand Python’s data science functionality. Xarray, which we explored in this chapter, is one example: it adds dimension names, coordinates, and attributes to NumPy arrays to simplify working with labeled, multidimensional data. In Chapter 4, we explore another library that builds on NumPy: Pandas.
3.12 Questions
- Explain two ways in which a NumPy array differs from a Python list.
- Explain the difference between signed and unsigned integer data types. Why might an 8-bit unsigned integer (
uint8) data type be appropriate for some image data? - For a three-dimensional array representing a multispectral image, explain what the
shape,ndim, andsizeproperties describe. - Explain the purpose of the
reshape()method. What does providing -1 for a dimension accomplish? - Explain the difference between a copy and a view of an array. Why does this difference matter?
- What is broadcasting? Give an example of an operation between two arrays with different shapes that is made possible by broadcasting.
- How does Xarray represent the dimensions of an array differently than NumPy? What is one benefit of this approach for geospatial data?
- Explain the difference between the
isel()andsel()methods in Xarray, including how they differ when slicing. - Explain the difference between a DataArray and a Dataset.