Pythonを用いたAISデータ分析

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

turfpy+ポリゴンを用いた分析

• ある領域内に船が何隻滞在しているか、何時間滞在しているのか調べたいことがままある
• Pythonのパッケージでそれが可能になる

tufpyとは？

• 公式ドキュメント
• 地理情報分析に便利
• Pythonで使える
• 2つの座標間の距離の算出やポリゴンを用いた操作などができる

#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を用いる
• 左のメニューバーからこの画像の青く光っているアイコンをクリック
  • 新しいプロジェクトの作成 > KMLファイルを作成

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

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

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

2. kmlファイルの操作

• まず取得したkmlファイルから座標を抽出する
  • pykmlを使うのが便利
  • 座標データが含まれる場所
    • Document.Placemark.Polygon.outerBoundaryIs.LinearRing.coordinates
  • ポリゴンの頂点の座標のx座標, y座標, z座標が空白区切りで保存されている
    • x1,y1,z1 x2,y2,z2 ... xn,yb,zn
  • z座標は使わないので除く

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が有効

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:経度

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

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

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

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の横にアップロードを表すアイコンが現れる

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

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