갈루아의 반서재

Numpy로 수치계산하기 (3)


NumPy 로 밀도지도 그리기


하루 중 특정 시간대의 승하차 위치에 대한 2D 밀도 지도를 계산하고 그려보자. 먼저, 택시 운행을 선택해보자. 데이터셋에서 저녁 운행수는 242,818건이다.

1
2
3
4
5
>>> evening = (data.pickup_datetime.dt.hour >= 19).values
>>> n = np.sum(evening)
>>> n
 
242818
cs


밀도 지도를 그리는 방법은 다음과 같다. 먼저 저녁 운행수 n 에 대한 모든 승하차 위치를 생각해본다. 이러한 위치에 대한 2n 을 만든다. 모든 위치에서 승차는 -1 의 가중치를, 하차는 +1 의 가중치를 부여한다. 주어진 위치에서 산술 밀도는 이 위치에서 출발하거나 도착한 사람을 반영한다.

2n 에 대한 weights 벡터를 생성하기 위해 0 으로 이루어진 벡터를 먼저 만들고 절반은 -1 (승차), 절반은 +1 (하차)로 설정한다.

다음 저녁 운행에 대한 승하차 위치를 열로 붙인 (2n, 2) 배열을 생성한다. np.r_ 문법은 첫 번째 열 차원을 따라 배열을 붙이게 한다. np.vstack( )  나 np.concatenate( ) 와 같은 조작 함수를 좀 더 명시적으로 사용할 수 있다.

1
2
3
4
5
6
7
8
9
>>> weights = np.zeros(2 * n)
>>> weights[:n] = -1
>>> weights[n:] = 1
>>> 
>>> points = np.r_[pickup[evening], dropoff[evening]]
>>>
>>> points.shape
 
(4856362)
cs


지리적 좌표를 픽셀 좌표로 변환해보자.

1
2
3
4
5
6
7
8
9
>>> def lat_lon_to_pixels(lat, lon):
>>>     lat_rad = lat * np.pi / 180.0
>>>     lat_rad = np.log(np.tan((lat_rad + np.pi / 2.0/ 2.0))
>>>     x = 100 * (lon + 180.0/ 360.0
>>>     y = 100 * (lat_rad - np.pi) / (2.0 * np.pi)
>>>     return (x, y)
>>> 
>>> lon, lat = points.T
>>> x, y = lat_lon_to_pixels(lat, lon)
cs


밀도 지도에 사용할 2D 히스토그램에 대한 구간bin 을 정의한다. 히스토그램을 계산하는 2D 격자를 정의한다. 이 두 배열은 수직과 수평 구간을 포함한다.

1
2
3
4
5
6
7
>>> lon_min, lat_min = -74.021440.6978
>>> lon_max, lat_max = -74.021440.6978
>>> x_min, y_min = lat_lon_to_pixels(lat_min, lon_min)
>>> x_max, y_max = lat_lon_to_pixels(lat_max, lon_max)
>>> bin = .00003
>>> bins_x = np.arange(x_min, x_max, bin)
>>> bins_y = np.arange(y_min, y_max, bin)
cs


np.histogram2d( ) 함수로 히스토그램을 계산한다. 점의 y, x 좌표 (격자의 첫 번째 축은 y 좌표를 나타내기 위해), weights, bin 을 인수로 전달한다. 이 함수는 모든 구간에서 점의 가중치 합을 계산한다. 몇 개의 객체를 반환하는 데 첫 번째가 우리가 원하는 밀도 지도이다.

1
>>> grid, _, _ = np.histogram2d(y, x, weights=weights, bins = (bins_y, bins_x))
cs


밀도 함수를 출력하기 전에 부드럽게 하기 위해 로지스틱 함수를 적용한다. 이 로지스틱 함수를 expit 함수라고 한다.

1
>>> density = 1. / (1. + np.exp(-.5 * grid))
cs


마지막으로 plt.imshow( ) 로 밀도 지도를 출력한다. 아래 그림에서 흰 지역은 하차 위치에 해당하며, 검은 지역은 승차 위치에 해당한다.

1
2
3
4
5
>>> plt.imshow(density, origin='lower', interpolation='bicubic')
>>>
>>> plt.axis('off')
 
(-0.5637.5-0.51226.5)
cs