갈루아의 반서재

Numpy로 수치계산하기 (2)


기본 배열 조작


먼저 1부터 10까지의 정수 배열을 만든다.

1
2
3
4
5
6
7
>>> import numpy as np
>>> 
>>> x = np.arange(111)
>>> 
>>> x
 
array([ 1,  2,  3,  4,  5,  6,  7,  8,  910])
cs


조작 테이블을 만들기 위해 먼저 x 를 행과 열로 변환하자. x 는 1D 배열이지만 행과 열 벡터는 2D 배열(매트릭스)이다.

1) reshape( ) 를 이용한 방법

1
2
3
4
5
>>> x_row = x.reshape((1-1))
>>> 
>>> x_row
>>>
array([[ 1,  2,  3,  4,  5,  6,  7,  8,  910]])
cs

reshape( ) 메소드는 파라메터로 새 크기를 받는다. 원소의 총 개수는 그대로 유지해야 한다. 예를 들어, (2, 3) 배열을 (5, 0) 배열로 변환하면 오류가 발생한다. -1 은 해당 축의 크기를 NumPy가 자동으로 인식하라는 것이다.

2개의 대괄호는 x_row 는 하나의 행으로 된 2D 배열을 나타내는 반면 x 는 1D 배열이다.

NumPy 에서 첫 번째 축은 수직행이고, 두 번째 축은 수평행이다.


2) 특수 인덱싱 문법을 사용하는 방법

1
2
3
4
5
6
7
8
9
10
11
12
13
14
>>> x_col = x[:, np.newaxis]
>>>
>>> x_col
 
array([[ 1],
       [ 2],
       [ 3],
       [ 4],
       [ 5],
       [ 6],
       [ 7],
       [ 8],
       [ 9],
       [10]])
cs

여기서 콜론( : )은 첫 번째 축 전체를 선택하고자한다는 의미이고, 새로운 두 번째 축을 만들기 위해 np.newaxis 를 사용했다.


이제 매트릭스를 만들어보자. 두 벡터 간에 매트릭스 곱을 하기 위해 np.dot( )을 사용한다. x_col (10, 1) 배열과 x_row (1, 10) 배열이기 때문에 매트릭스 곱은 (10, 10) 배열이다.

1
2
3
4
5
6
7
8
9
10
11
12
>>> np.dot(x_col, x_row)
>>>
array([[  1,   2,   3,   4,   5,   6,   7,   8,   9,  10],
       [  2,   4,   6,   8,  10,  12,  14,  16,  18,  20],
       [  3,   6,   9,  12,  15,  18,  21,  24,  27,  30],
       [  4,   8,  12,  16,  20,  24,  28,  32,  36,  40],
       [  5,  10,  15,  20,  25,  30,  35,  40,  45,  50],
       [  6,  12,  18,  24,  30,  36,  42,  48,  54,  60],
       [  7,  14,  21,  28,  35,  42,  49,  56,  63,  70],
       [  8,  16,  24,  32,  40,  48,  56,  64,  72,  80],
       [  9,  18,  27,  36,  45,  54,  63,  72,  81,  90],
       [ 10,  20,  30,  40,  50,  60,  70,  80,  90100]])
cs


다른 방법으로 * 기호로 NumPy 곱셈을 하는 방법이 있다. 배열에서 이 연산은 매트릭스 곱셈이 아닌 원소 단위 곱셈이다.

1
2
3
4
5
6
7
8
9
10
11
12
>>> x_row * x_col
>>>
array([[  1,   2,   3,   4,   5,   6,   7,   8,   9,  10],
       [  2,   4,   6,   8,  10,  12,  14,  16,  18,  20],
       [  3,   6,   9,  12,  15,  18,  21,  24,  27,  30],
       [  4,   8,  12,  16,  20,  24,  28,  32,  36,  40],
       [  5,  10,  15,  20,  25,  30,  35,  40,  45,  50],
       [  6,  12,  18,  24,  30,  36,  42,  48,  54,  60],
       [  7,  14,  21,  28,  35,  42,  49,  56,  63,  70],
       [  8,  16,  24,  32,  40,  48,  56,  64,  72,  80],
       [  9,  18,  27,  36,  45,  54,  63,  72,  81,  90],
       [ 10,  20,  30,  40,  50,  60,  70,  80,  90100]])
cs



NumPy 로 배열계산


NumPy 로 배열을 다뤄보자.NYC 택시 데이터를 로드하여 판다 DataFrame 의 .values 로 ndarray 로 택시 승하차 위치를 구한다.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
>>> import numpy as np
>>> import pandas as pd
>>> import matplotlib.pyplot as plt
>>> import seaborn as sns
>>> %matplotlib inline
>>> 
>>> data = pd.read_csv('data/nyc_data.csv', parse_dates = ['pickup_datetime''dropoff_datetime'])
>>> 
>>> pickup = data[['pickup_longitude''pickup_latitude']].values
>>> dropoff = data[['dropoff_longitude''dropoff_latitude']].values
>>> pickup
 
array([[-73.955925,  40.781887],
       [-74.005501,  40.745735],
       [-73.969955,  40.79977 ],
       ..., 
       [-73.993492,  40.729347],
       [-73.978477,  40.772945],
       [-73.987206,  40.750568]])
cs

NumPy 에서는 pickup[ i, j ] 와 같은 형태로 주어진 원소를 검색할 수 있다. pickup[ i, j ] 는 0 부터 시작하는 i 번째 행이며, 0 부터 시작하는 j 번째 열이다.

그리고 라인 5에서와 같이 시작 위치, 긑 위치, 간격을 선택하여 슬라이싱도 가능하다. 여기서는 [1, 1], [3, 1], [5, 1]에 있는 원소를 선택했다. 파이썬 슬라이싱 문법에서는 end 에 해당하는 값은 포함하지 않는다. start 나 end 가 생략되면 start 는 시작과 end 는 끝으로 설정된다. step 은 기본 1이다. 예를 들어, 1:  은 1:n:1 과 같고 n 은 축의 길이이다.

1
2
3
4
5
6
7
8
9
>>> print(pickup[31])
 
40.755081
 
>>> pickup[1:7:2,1:]
 
array([[ 40.745735],
       [ 40.755081],
       [ 40.768978]])
cs


실제 CSV 데이터에서 보면 다음과 같다.

모든 승차 위치의 경도와 위도를 다음과 같이 구할 수 있다.

1
2
3
4
5
6
7
8
9
10
11
>>> lon = pickup[:, 0]
>>> lon
 
array([-73.955925-74.005501-73.969955, ..., -73.993492-73.978477,
       -73.987206])
 
>>> lat = pickup[:, 1]
>>> lat
 
array([ 40.781887,  40.745735,  40.79977 , ...,  40.729347,  40.772945,
        40.750568])
cs



NumPy 필터링 연산


주어진 위치에서 출발한 모든 운행을 선택해보자. lon_min 과 lon_max 사이에 있는 경도로, 모든 운행에서 선택해보자. 결과는 lon 벡터에 있는 원소에 대한 Boolean 벡터다.

1
2
3
4
5
6
7
>>> lon_min, lon_max = (-73.98330-73.98025)
>>> lat_min, lat_max = (40.7672440.76871)
>>> 
>>> in_lon = (lon_min <= lon) & (lon <= lon_max)
>>> in_lon
 
array([False, False, False, ..., False, False, False], dtype=bool)
cs

sum( ) 메소드를 이용하여 모든 원소의 합을 반환할 수 있다. False 원소는 0 으로, 그리고 True 원소는 1 로 변환한다. 그러므로, 합은 True 인 원소의 개수가 된다.

1
2
3
>>> in_lon.sum()
 
69163
cs

마찬가지로 위도도 처리해보자.

1
2
3
4
>>> in_lat = (lat_min <= lat) & (lat <= lat_max)
>>> in_lat.sum()
 
16110
cs

위도와 경도에 해당하는 사각에 속한 모든 운행을 구해보자.

1
2
3
4
>>> in_lonlat = in_lon & in_lat
>>> in_lonlat.sum()
 
3998
cs

np.nonzero( ) 함수는 다음과 같이 불린 배열에서 True에 해당하는 인덱스를 반환한다.

1
2
3
>>> np.nonzero(in_lonlat)[0]
 
array([   901,   1011,   1066, ..., 845749845903846080])
cs

마지막으로 하차 좌표가 필요하다. T 속성은 매트릭스의 전치이다. 여기서 dropoff.T 는 (2, N)의 배열로, 첫 행은 경도이며, 두 번째 행은 위도이다.

1
2
3
4
5
6
7
>>> lon1, lat1 = dropoff.T
>>> lon1, lat1
 
(array([-73.963181-73.964943-73.954567, ..., -74.013725-73.963814,
        -73.970909]),
 array([ 40.777832,  40.755722,  40.787392, ...,  40.702332,  40.773922,
         40.795815]))
cs



배열에서의 수학연산


모든 승하차 위치를 NumPy 배열로 만들었다. 모든 운행에서 두 위치간의 거리를 구해보자. 두 위치간 거리를 구하는 공식은 아래 링크를 참고한다.

※ Formula and code for calculating distance based on two lat/lon locations
http://jan.ucc.nau.edu/~cvm/latlon_formula.html


모든 택시 운행의 직선거리를 계산하여 히스토그램으로 나타내보면 아래와 같다.

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
>>> EARTH_R = 6372.8
>>> 
>>> def geo_distance(lon0, lat0, lon1, lat1):
>>>     lat0 = np.radians(lat0)
>>>     lon0 = np.radians(lon0)
>>>     lat1 = np.radians(lat1)
>>>     lon1 = np.radians(lon1)
>>>     dlon = lon0 - lon1
>>>     y = np.sqrt(
>>>         (np.cos(lat1) * np.sin(dlon)) ** 2
>>>         + (np.cos(lat0) * np.sin(lat1)
>>>         - np.sin(lat0) * np.cos(lat1) * np.cos(dlon)) ** 2)
>>>     x = np.sin(lat0) * np.sin(lat1) + np.cos(lat0) * np.cos(lat1) * np.cos(dlon)
>>>     c = np.arctan2(y, x)
>>>     return EARTH_R * c
>>> 
>>> distances = geo_distance(lon, lat, lon1, lat1)
>>> 
>>> lon[:10], lat[:10], lon1[:10], lat1[:10], distances[:10]
 
(array([-73.955925-74.005501-73.969955-73.991432-73.966225,
        -73.955238-73.98558 , -73.999413-73.99218 , -74.006554]),
 array([ 40.781887,  40.745735,  40.79977 ,  40.755081,  40.773716,
         40.768978,  40.752197,  40.730137,  40.738644,  40.74054 ]),
 array([-73.963181-73.964943-73.954567-73.991417-73.955399,
        -73.994064-73.988777-73.983047-73.973099-73.980469]),
 array([ 40.777832,  40.755722,  40.787392,  40.755085,  40.782597,
         40.720299,  40.737244,  40.741844,  40.749355,  40.730408]),
 array([  7.59535877e-01,   3.59342825e+00,   1.89062515e+00,
          1.33984499e-03,   1.34431257e+00,   6.32615091e+00,
          1.68484565e+00,   1.89684673e+00,   2.00118747e+00,
          2.47044801e+00]))
 
>>> plt.hist(distances[in_lonlat], np.linspace(0.10.50))
>>> plt.xlabel('Trip distance (km)')
>>> plt.ylabel('Number of trips')
>>> 
<matplotlib.text.Text at 0x7f7c6ff29c10>
cs