Pythonを用いたAISデータ分析

 AISデータは表形式のデータなのでpandasで取り扱えると便利
 Pythonを用いたAISデータハンドリングのtips

## turfpy+ポリゴンを用いた分析
- ある領域内に船が何隻滞在しているか、何時間滞在しているのか調べたいことがままある
- Pythonのパッケージでそれが可能になる

### tufpyとは？
- [公式ドキュメント](https://turfpy.readthedocs.io/en/latest/)
- 地理情報分析に便利
- Pythonで使える
- 2つの座標間の距離の算出やポリゴンを用いた操作などができる
```python
#pointがPolygon内にあるかどうかを真偽値で返す
from turfpy.measurement import boolean_point_in_polygon
from geojson import Point, MultiPolygon, Feature
point = Feature(geometry=Point([-77, 44]))
polygon = Feature(geometry=MultiPolygon([([(-81, 41), (-81, 47), (-72, 47), (-72, 41), (-81, 41)],),
([(3.78, 9.28), (-130.91, 1.52), (35.12, 72.234), (3.78, 9.28)],)]))
boolean_point_in_polygon(point, polygon)
>> True
```

### 分析の手順
#### 1. ポリゴンの作成 
- 座標を取得するために[Google Earth](https://earth.google.com/web)を用いる
- 左のメニューバーからこの画像の青く光っているアイコンをクリック
    - `新しいプロジェクトの作成` > `KMLファイルを作成`

![](https://i.imgur.com/Amvb8dx.png)

- `アイテムを追加` > `ラインやシェイプを描画`

![](https://i.imgur.com/prmYeVG.jpg)

- するとこのようにポリゴンの描画ができる

![](https://i.imgur.com/BrIvO7Q.jpg)

![](https://i.imgur.com/jpJVhPT.jpg)

- 左上の戻るボタンを押してゴミ箱横の3点を押し、`KMLファイルとしてエクスポート`をクリックするとkmlファイルがダウンロードされる

![](https://i.imgur.com/cFaPLPH.png)

#### 2. kmlファイルの操作
- まず取得したkmlファイルから座標を抽出する
    - pykmlを使うのが便利
    - 座標データが含まれる場所
        - `Document.Placemark.Polygon.outerBoundaryIs.LinearRing.coordinates`
    - ポリゴンの頂点の座標のx座標, y座標, z座標が空白区切りで保存されている
        - `x1,y1,z1 x2,y2,z2 ... xn,yb,zn`
    - z座標は使わないので除く
```python
from pykml import parser
#pathはkmlファイルのパス
with open(path, 'r', encoding='utf-8') as f: 
    doc = parser.parse(f).getroot()
#ポリゴンの頂点の座標を取得
location = str(doc.Document.Placemark.Polygon.outerBoundaryIs.LinearRing.coordinates).replace('\n', '').replace('\t', '')
#分割してx,y座標のみ取得, 各座標のtupleのリストを得る
location = location.split(' ')[:-1]
for i in range(len(location)):
    location[i] = tuple(map(float, location[i].split(',')[:2]))
location

>> [(121.2439643593101, 30.58827889879364),
    (122.0651397675725, 30.86645010206499),
    (122.160478178912, 31.18354126504846),
    (121.9142857018593, 31.8072456187944),
    (121.1425848528823, 31.81215480212967),
    (121.2439643593101, 30.58827889879364)]
```

#### 3. turfpyによる分析
- 取得した座標からポリゴンを作成し、ポリゴン内にその点が含まれるかどうかを調べたい
- `turfpy.measurement`の`boolean_point_in_polygon`が有効

```python
from turfpy.measurement import boolean_point_in_polygon
from geojson import Point, MultiPolygon, Feature

#pointがpolygonに含まれるかどうかを真偽値で返す
point = Feature(geometry=Point([121.5, 30.8]))
polygon = Feature(geometry=Polygon([location]))
boolean_point_in_polygon(point, polygon)
>> True
```

- `geojson`の`Multipolygon`を用いることで複数個のポリゴンを渡すこともできる


#### 4. AISデータとの合わせ技
- AISデータは各行に緯度経度の情報が含まれている
    - latitude:緯度, longtitude:経度
    
![](https://i.imgur.com/Okl395E.png)

- 各行の座標がポリゴン内に含まれているかどうか調べる
    - is_in_polygonという列を追加する
    - ポリゴン内にあればTrue, なければFalse

```python
import pandas as pd

df_ais = pd.read_csv(path)
df_ais['is_in_polygon'] = [boolean_point_in_polygon(Feature(geometry=Point(list(point))), polygon) for point in list(df[['Lng', 'Lat']].values)]
#リスト内包表記を用いると高速化
df
```
![](https://i.imgur.com/DzAy3Lm.png)

- これを用いて2020/3/1にポリゴン内に滞在した形跡のある船舶数を数えてみる
    - 各船で2020/3/1 0:00 ~ 2020/3/2 0:00にTrueとなっているかどうかを調べれば良い

```python
from datetime import datetime
import re

#日付をstr型からdatetimeに変換
def convert_datetime(date):
    year, month, day, hour, minute, second = list(map(int, re.findall(r'\d+', date)))
    return datetime(year, month, day, hour, minute, second)
 
df['Date/Time'] = [convert_datetime(date) for date in df['Date/Time'].values]

#ポリゴン内に滞在した形跡のある船のリスト
vessels_in_polygon = []

#日付で絞込
df2 = df[(datetime(2020,3,1)<=df['Date/Time']) & (df['Date/Time']<datetime(2020,3,2))]

#各船で滞在履歴を見る
for imo in df2['IMO'].unique():
    df_imo = df2[df2['IMO']==imo]
    if df_imo['is_in_polygon'].any():
        vessels_in_polygon.append(imo)
        
>> [11016629, 11107608, 11185908, 12247559,
    12590648, 367915, 372253, 373010, 
    411031, 9997379]
```

## おまけ
- ダウウンロードしたkmlファイルはseasearcherにも取り込める
- `Draw polygon area`をクリックすると`Last Known`の横にアップロードを表すアイコンが現れる

![](https://i.imgur.com/HYSY6s1.png)

- アップロード画面が現れるのでここにkmlファイルをアップロードすればよい

![](https://i.imgur.com/DMds3bf.png)

- seasearcherからkmlファイルをダウンロードすることはできないので注意