はじめに

バスの現在位置は、Google Mapからでも見られる時代ですが、自分の手でマッピングしてみましょう。pythonにfoliumというWeb地図表示プラグインを用いて表示してみます。
地図のアイコンはクリッカブルで、詳細な情報を表示することも可能です。

実装

realtimebus.py3

from google.transit import gtfs_realtime_pb2
import urllib.request, urllib.error
import folium
import pandas as pd

#GTFS-RTの公開されているURL(宇野バスさん）
url ='http://www3.unobus.co.jp/GTFS/GTFS_RT-VP.bin'
#カラム名の宣言
list_df = pd.DataFrame(columns=['id' , 'vehicle_id', 'trip_id','vehicle_timestamp','longitude','latitude','occupancy_status'])
feed = gtfs_realtime_pb2.FeedMessage()

#データダウンロードとフォーマット変換
with urllib.request.urlopen(url) as res: #データダウンロード
    feed.ParseFromString(res.read()) #プロトコルバッファをDeserialize
    for entity in feed.entity:
        tmp_se = pd.Series( [
              entity.id,                         #車両ID
              entity.vehicle.vehicle.id,         #車両ナンバー
              entity.vehicle.trip.trip_id,       #路線番号？
              entity.vehicle.timestamp,          #車両時刻
              entity.vehicle.position.longitude, #車両緯度
              entity.vehicle.position.latitude,  #車両経度
              entity.vehicle.occupancy_status   #混雑度
        ], index=list_df.columns )
        list_df = list_df.append( tmp_se, ignore_index=True )

#地図の中心地点を算出
average_pos = list_df.mean()
# 中心座標付近の地図を作成
m = folium.Map(location=[average_pos['latitude'], average_pos['longitude']], zoom_start=11)
# 各バスをプロット(色は青)
list_df.apply(lambda row:
             folium.Marker(
                 location=[row['latitude'], row['longitude']],
                 popup='<table border="1"><tr><th>occupancy_status</th></tr><tr><td>'
                 +str(row['occupancy_status'])
                 +'</td></tr></table>',
                 icon=folium.Icon(color='blue',icon='bus', prefix='fa')
             ).add_to(m), axis=1)
m

結果

ダウンロードされたリアルタイム位置情報が地図上にプロットされました

実行したタイミングは２２時半過ぎてますが、まだバス動いているのですね。。。

Youtube

動画解説もしています。

問題

P-022: レシート明細データフレーム（df_receipt）の顧客ID（customer_id）に対し、ユニーク件数をカウントせよ。

解答

コード

len(df_receipt['customer_id'].unique())

出力

8307

解説

・PandasのDataFrame/Seriesにて、ユニークな行数をカウントする方法です。
・ユニークな行数が全部で何件あるのかを確認したい時に使用します。

※ちなみに、以下のコードもユニークな件数をカウントするコードであり、同じ結果を出力します

コード

len(df_receipt['customer_id'].value_counts())

はじめに

「Djangoの管理者ページのカスタマイズなんて、今までしたことないよ〜」という方向け。
初級編です。

カスタマイズしなくても十分使えそうだが、エンジニア以外の人が触るようになると、デフォルトだと少し寂しい？

詳細な解説というより、ググると良い用語をパパッと上げていくスタイルです。
ちゃんとまとめると大変ですからね。。。

Djangoのドキュメント
https://docs.djangoproject.com/en/3.0/ref/contrib/admin/

ModelAdminクラス

基本的には、ModelAdminクラスのサブクラスを作成して、カスタマイズ内容を記述する。
例えばSampleModelという名前のモデルをカスタマイズしようと思ったら、以下のようにする。

admin.py

from django.contrib import admin

# 通常は「モデル名＋Admin」という名前にする
class SampleModelAdmin(admin.ModelAdmin):
    '''
    このクラスの中にカスタマイズ内容を記述
    '''

# SampleModelにSampleModelAdminを適用
admin.site.register(SampleModel, SampleModelAdmin)

これ以降は、ModelAdminクラスの中の記述のみを書く。
例えば

def hoge():
    return

と書いたら、実際は

admin.py

class SampleModelAdmin(admin.ModelAdmin):
    def hoge():
        return

と読み替えてください。

どんなカラムを表示するか

レコード一覧の一番上にカラム名が表示されるが、どんなカラムを表示したいか設定しよう。
以下はモデルがtext,updated_atというフィールドを持っていた場合の例。

list_display = ('__str__', 'text', 'updated_at')  # この３つのカラムを表示する

次に、フィールドにはない値を表示する場合。

list_display = ('__str__', 'text', 'updated_at', 'get_owner_name')

def get_owner_name(self, obj):
    # 第二引数のobjはレコード自身を表す
    return "I am owner"
get_owner_name.short_description = 'オーナー名'

こうすると「オーナー名」というカラムが追加され、この場合は全てのレコードに「I am owner」と表示される。

検索対象を設定する

上部に検索ボックスが存在するが、その検索対象を絞ることができる。

# SampleModel.text, SampleModel.updated_at, SampleModel.owner.username を検索対象にする
search_fields = ('text', 'updated_at', 'owner__username`)

操作を追加する

デフォルトでは「選択された<モデル名>の削除」という操作が存在するはず。
自作のアクションは以下のようにactionsにリスト形式で記述するが、ここに明示的に含めなくても「選択された<モデル名>の削除」は存在し続ける。

actions = ['my_custom_action']

def my_custom_action(self, request, queryset):
    # なんか処理
    return None
my_custom_action.short_description = '私だけのアクション'

return Noneとすれば操作完了時に何もしないが、例えばテンプレートを返してページ遷移なんかもできる。

モデルの表示名を変更する

これはModelAdminではなく元のModelクラスに書いておく。

class SampleModel(models.Model):
     class Meta:
         verbose_name = 'サンプルモデル'
         verbose_name_plural = 'サンプルモデル一覧'

フィールドにリンクを貼るか

list_display_links = ('text',)

応用編：検索処理のカスタマイズ

検索対象のカスタマイズを紹介しましたが、検索ワード（クエリ）が与えられた時に、実際にどのような検索を行うかもカスタマイズできます。
これは長くなるので概要だけにしますが、以下の関数を定義して検索処理を上書きします。

def get_changelist(self, request, **kwargs):

最後に

以上のカスタマイズができれば、はじめはそこそこ満足ではないでしょうか？
adminページのテンプレート自体のカスタマイズや、admi内でのページ遷移については別の記事で書こうかと思います。

• 決定木における分類の分割基準はいくつかあります。
• 一般的に使用されているものに、ジニ不純度（gini impurity）、エントロピー（entropy）という２つの不純物測定と、それから分類誤差（misclassification error）という分割基準があります。
• これらの違いを見ていきます。

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

２クラス分類での各指標の比較

• ２クラス分類で、一方のクラスの割合を p とすると、各指標はそれぞれ次のように求められます。
• ジニ(gini): $2p * (1-p)$
• エントロピー(entropy): $-p * \log p-(1-p) * log(1-p)$
• 誤分類率(misclassification rate): $1-max(p, 1-p)$

# p に相当する等差数列を生成
xx = np.linspace(0, 1, 50) # 開始値0、終了値1、要素数50

plt.figure(figsize=(10, 8))

# 各指標を計算
gini = [2 * x * (1-x) for x in xx]
entropy = [-x * np.log2(x) - (1-x) * np.log2(1-x)  for x in xx]
misclass = [1 - max(x, 1-x) for x in xx]

# グラフを表示
plt.plot(xx, gini, label='gini', lw=3, color='b')
plt.plot(xx, entropy, label='entropy', lw=3, color='r')
plt.plot(xx, misclass, label='misclass', lw=3, color='g')

plt.xlabel('p', fontsize=15)
plt.ylabel('criterion', fontsize=15)
plt.legend(fontsize=15)
plt.xticks(fontsize=12) 
plt.yticks(fontsize=12) 
plt.grid()

• ２クラスの割合が均等のとき(p=0.5)、いずれも最大値となりますが、
• エントロピーの最大値は 1.0、ジニは 0.5 なので、比較しやすいようにエントロピーの方を1/2の係数でスケーリングします。

# p に相当する等差数列を生成
xx = np.linspace(0, 1, 50) # 開始値0、終了値1、要素数50

plt.figure(figsize=(10, 8))

# 各指標を計算
gini = [2 * x * (1-x) for x in xx]
entropy = [(x * np.log((1-x)/x) - np.log(1-x)) / (np.log(2))  for x in xx]
entropy_scaled = [(x * np.log((1-x)/x) - np.log(1-x)) / (2*np.log(2))  for x in xx]
misclass = [1 - max(x, 1-x) for x in xx]

# グラフを表示
plt.plot(xx, gini, label='gini', lw=3, color='b')
plt.plot(xx, entropy, label='entropy', lw=3, color='r', linestyle='dashed')
plt.plot(xx, entropy_scaled, label='entropy(scaled)', lw=3, color='r')
plt.plot(xx, misclass, label='misclass', lw=3, color='g')

plt.xlabel('p', fontsize=15)
plt.ylabel('criterion', fontsize=15)
plt.legend(fontsize=15)
plt.xticks(fontsize=12) 
plt.yticks(fontsize=12) 
plt.grid()

• ジニ(gini)とエントロピー(entropy)はよく似ていて、どちらも二次曲線を描いて p = 1/2 で最大化しています。
• 誤分類率(misclass)は、線形であるという点で明らかに異なりますが、多くの同じ特性を共有しています。p = 1/2 で最大化し、p = 0, 1 でゼロに等しく、山なりに変化しています。

情報利得にみる各指標の違い

1. 情報利得とは

• 情報利得(information gain)とは、ある変数を使ってデータを分割したとき、分割の前後でどれだけ不純度が減少したかを表す指標です。
• 不純度が減れば減るほど、その変数は分類条件として有益だということになります。
• その意味で、いかに不純度を低減できるか、ということは各分岐における情報利得の最大化と同義といえます。
• ２クラス分類の場合、情報利得(IG)は次の式に定義されます。
• $\displaystyle IG(D_{p}, a) = I(D_{p}) - \frac{N_{left}}{N} I(D_{left}) - \frac{N_{right}}{N} I(D_{right})$
• $I$は不純度を意味し、$D_{p}$は親(parent)ノードのデータ、子ノードのデータの左右が$D_{left}$と$D_{right}$で、また親ノードのサンプル数$N$を分母として子ノードの左右それぞれの割合となっています。

2. 各指標での情報利得を計算する

• 次のような２通りの分岐条件A, Bを仮定します。
• ジニ、エントロピー、誤分類率のそれぞれについて、分岐条件A, B間での情報利得を比較します。

➀ジニ不純度による情報利得

# 親ノードのジニ不純度
IGg_p = 2 * 1/2 * (1-(1/2))

# 子ノードAのジニ不純度
IGg_A_l = 2 * 3/4 * (1-(3/4)) # 左
IGg_A_r = 2 * 1/4 * (1-(1/4)) # 右
# 子ノードBのジニ不純度
IGg_B_l = 2 * 2/6 * (1-(2/6)) # 左
IGg_B_r = 2 * 2/2 * (1-(2/2)) # 右

# 各分岐の情報利得
IG_gini_A = IGg_p - 4/8 * IGg_A_l - 4/8 * IGg_A_r
IG_gini_B = IGg_p - 6/8 * IGg_B_l - 2/8 * IGg_B_r

print("分岐Aの情報利得：", IG_gini_A)
print("分岐Bの情報利得：", IG_gini_B)

➁エントロピーによる情報利得

# 親ノードのエントロピー
IGe_p = (4/8 * np.log((1-4/8)/(4/8)) - np.log(1-4/8)) / (np.log(2))

# 子ノードAのエントロピー
IGe_A_l = (3/4 * np.log((1-3/4)/(3/4)) - np.log(1-3/4)) / (np.log(2)) # 左
IGe_A_r = (1/4 * np.log((1-1/4)/(1/4)) - np.log(1-1/4)) / (np.log(2)) # 右
# 子ノードBのエントロピー
IGe_B_l = (2/6 * np.log((1-2/6)/(2/6)) - np.log(1-2/6)) / (np.log(2)) # 左
IGe_B_r = (2/2 * np.log((1-2/2+1e-7)/(2/2)) - np.log(1-2/2+1e-7)) / (np.log(2)) # 右、+1e-7は0除算回避のために微小値を加算

# 各分岐の情報利得
IG_entropy_A = IGe_p - 4/8 * IGe_A_l - 4/8 * IGe_A_r
IG_entropy_B = IGe_p - 6/8 * IGe_B_l - 2/8 * IGe_B_r
print("分岐Aの情報利得：", IG_entropy_A)
print("分岐Bの情報利得：", IG_entropy_B)

➂誤分類率による情報利得

# 親ノードの誤分類率
IGm_p = 1 - np.maximum(4/8, 1-4/8)

# 子ノードAの誤分類率
IGm_A_l = 1 - np.maximum(3/4, 1-3/4) # 左
IGm_A_r = 1 - np.maximum(1/4, 1-1/4) # 右
# 子ノードBの誤分類率
IGm_B_l = 1 - np.maximum(2/6, 1-2/6) # 左
IGm_B_r = 1 - np.maximum(2/2, 1-2/2) # 右

# 各分岐の情報利得
IG_misclass_A = IGm_p - 4/8 * IGm_A_l - 4/8 * IGm_A_r
IG_misclass_B = IGm_p - 6/8 * IGm_B_l - 2/8 * IGm_B_r
print("分岐Aの情報利得：", IG_misclass_A)
print("分岐Bの情報利得：", IG_misclass_B)

まとめ

	分類条件A	分類条件B
ジニ不純度	0.125	0.167
エントロピー	0.189	0.311
誤分類率	0.250	0.250

• ジニとエントロピーでは共に、分類条件Bの方が優先されるのは明白です。実際、結果も非常に似たものになります。
• 一方、誤分類率では、分類条件AとBがほぼ同率となっています。
• このように、誤分類率は変数間で明らかな差分をとらない傾向があることから、sklearnの決定木モデルでは、分割基準をジニとエントロピーの２つとしている（らしいです）。
• なお、ジニの方が、エントロピーのような対数の計算がない分だけ少し処理が速いかも知れませんね。

• 決定木における分類の分割基準はいくつかあります。
• 一般的に使用されているものに、ジニ不純度（gini impurity）、エントロピー（entropy）という２つの不純物測定と、それから分類誤差（misclassification error）という分割基準があります。
• これらの違いを見ていきます。

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

２クラス分類での各指標の比較

• ２クラス分類で、一方のクラスの割合を p とすると、各指標はそれぞれ次のように求められます。
• ジニ(gini): $2p * (1-p)$
• エントロピー(entropy): $-p * \log p-(1-p) * log(1-p)$
• 誤分類率(misclassification rate): $1-max(p, 1-p)$

# p に相当する等差数列を生成
xx = np.linspace(0, 1, 50) # 開始値0、終了値1、要素数50

plt.figure(figsize=(10, 8))

# 各指標を計算
gini = [2 * x * (1-x) for x in xx]
entropy = [-x * np.log2(x) - (1-x) * np.log2(1-x)  for x in xx]
misclass = [1 - max(x, 1-x) for x in xx]

# グラフを表示
plt.plot(xx, gini, label='gini', lw=3, color='b')
plt.plot(xx, entropy, label='entropy', lw=3, color='r')
plt.plot(xx, misclass, label='misclass', lw=3, color='g')

plt.xlabel('p', fontsize=15)
plt.ylabel('criterion', fontsize=15)
plt.legend(fontsize=15)
plt.xticks(fontsize=12) 
plt.yticks(fontsize=12) 
plt.grid()

• ２クラスの割合が均等のとき(p=0.5)、いずれも最大値となりますが、
• エントロピーの最大値は 1.0、ジニは 0.5 なので、比較しやすいようにエントロピーの方を1/2の係数でスケーリングします。

# p に相当する等差数列を生成
xx = np.linspace(0, 1, 50) # 開始値0、終了値1、要素数50

plt.figure(figsize=(10, 8))

# 各指標を計算
gini = [2 * x * (1-x) for x in xx]
entropy = [(x * np.log((1-x)/x) - np.log(1-x)) / (np.log(2))  for x in xx]
entropy_scaled = [(x * np.log((1-x)/x) - np.log(1-x)) / (2*np.log(2))  for x in xx]
misclass = [1 - max(x, 1-x) for x in xx]

# グラフを表示
plt.plot(xx, gini, label='gini', lw=3, color='b')
plt.plot(xx, entropy, label='entropy', lw=3, color='r', linestyle='dashed')
plt.plot(xx, entropy_scaled, label='entropy(scaled)', lw=3, color='r')
plt.plot(xx, misclass, label='misclass', lw=3, color='g')

plt.xlabel('p', fontsize=15)
plt.ylabel('criterion', fontsize=15)
plt.legend(fontsize=15)
plt.xticks(fontsize=12) 
plt.yticks(fontsize=12) 
plt.grid()

• ジニ(gini)とエントロピー(entropy)はよく似ていて、どちらも二次曲線を描いて p = 1/2 で最大化しています。
• 誤分類率(misclass)は、線形であるという点で明らかに異なりますが、多くの同じ特性を共有しています。p = 1/2 で最大化し、p = 0, 1 でゼロに等しく、山なりに変化しています。

情報利得にみる各指標の違い

1. 情報利得とは

• 情報利得(information gain)とは、ある変数を使ってデータを分割したとき、分割の前後でどれだけ不純度が減少したかを表す指標です。
• 不純度が減れば減るほど、その変数は分類条件として有益だということになります。
• その意味で、いかに不純度を低減できるか、ということは各分岐における情報利得の最大化と同義といえます。
• ２クラス分類の場合、情報利得(IG)は次の式に定義されます。
• $\displaystyle IG(D_{p}, a) = I(D_{p}) - \frac{N_{left}}{N} I(D_{left}) - \frac{N_{right}}{N} I(D_{right})$
• $I$は不純度を意味し、$D_{p}$は親(parent)ノードのデータ、子ノードのデータの左右が$D_{left}$と$D_{right}$で、また親ノードのサンプル数$N$を分母として子ノードの左右それぞれの割合となっています。

2. 各指標での情報利得を計算する

• 次のような２通りの分岐条件A, Bを仮定します。
• ジニ、エントロピー、誤分類率のそれぞれについて、分岐条件A, B間での情報利得を比較します。

➀ジニ不純度による情報利得

# 親ノードのジニ不純度
IGg_p = 2 * 1/2 * (1-(1/2))

# 子ノードAのジニ不純度
IGg_A_l = 2 * 3/4 * (1-(3/4)) # 左
IGg_A_r = 2 * 1/4 * (1-(1/4)) # 右
# 子ノードBのジニ不純度
IGg_B_l = 2 * 2/6 * (1-(2/6)) # 左
IGg_B_r = 2 * 2/2 * (1-(2/2)) # 右

# 各分岐の情報利得
IG_gini_A = IGg_p - 4/8 * IGg_A_l - 4/8 * IGg_A_r
IG_gini_B = IGg_p - 6/8 * IGg_B_l - 2/8 * IGg_B_r

print("分岐Aの情報利得：", IG_gini_A)
print("分岐Bの情報利得：", IG_gini_B)

➁エントロピーによる情報利得

# 親ノードのエントロピー
IGe_p = (4/8 * np.log((1-4/8)/(4/8)) - np.log(1-4/8)) / (np.log(2))

# 子ノードAのエントロピー
IGe_A_l = (3/4 * np.log((1-3/4)/(3/4)) - np.log(1-3/4)) / (np.log(2)) # 左
IGe_A_r = (1/4 * np.log((1-1/4)/(1/4)) - np.log(1-1/4)) / (np.log(2)) # 右
# 子ノードBのエントロピー
IGe_B_l = (2/6 * np.log((1-2/6)/(2/6)) - np.log(1-2/6)) / (np.log(2)) # 左
IGe_B_r = (2/2 * np.log((1-2/2+1e-7)/(2/2)) - np.log(1-2/2+1e-7)) / (np.log(2)) # 右、+1e-7は0除算回避のために微小値を加算

# 各分岐の情報利得
IG_entropy_A = IGe_p - 4/8 * IGe_A_l - 4/8 * IGe_A_r
IG_entropy_B = IGe_p - 6/8 * IGe_B_l - 2/8 * IGe_B_r
print("分岐Aの情報利得：", IG_entropy_A)
print("分岐Bの情報利得：", IG_entropy_B)

➂誤分類率による情報利得

# 親ノードの誤分類率
IGm_p = 1 - np.maximum(4/8, 1-4/8)

# 子ノードAの誤分類率
IGm_A_l = 1 - np.maximum(3/4, 1-3/4) # 左
IGm_A_r = 1 - np.maximum(1/4, 1-1/4) # 右
# 子ノードBの誤分類率
IGm_B_l = 1 - np.maximum(2/6, 1-2/6) # 左
IGm_B_r = 1 - np.maximum(2/2, 1-2/2) # 右

# 各分岐の情報利得
IG_misclass_A = IGm_p - 4/8 * IGm_A_l - 4/8 * IGm_A_r
IG_misclass_B = IGm_p - 6/8 * IGm_B_l - 2/8 * IGm_B_r
print("分岐Aの情報利得：", IG_misclass_A)
print("分岐Bの情報利得：", IG_misclass_B)

まとめ

	分類条件A	分類条件B
ジニ不純度	0.125	0.167
エントロピー	0.189	0.311
誤分類率	0.250	0.250

• ジニとエントロピーでは共に、分類条件Bの方が優先されるのは明白です。実際、結果も非常に似たものになります。
• 一方、誤分類率では、分類条件AとBがほぼ同率となっています。
• このように、誤分類率は変数間で明らかな差分をとらない傾向があることから、sklearnの決定木モデルでは、分割基準をジニとエントロピーの２つとしている（らしいです）。
• なお、ジニの方が、エントロピーのような対数の計算がない分だけ少し処理が速いかも知れませんね。

偉大な先人

既に伊藤先生が、GTFSファイルをマージしてPostgreSQL+PostGISに投入する方法、ツールを活用してPostgreSQLにデータを投入する方法を掲載しております。

今回のモチベーション

日本のGTFSデータが、GCSに保存されているとの情報を入手したので、そのデータで遊ぶための準備です。

今回やること

Google Clound Platform のファイルサービスである、Google Cloud Strage 保存されている protocol buffers 形式の大量のファイルを、一気にダウンロードして一纏めのデータに変換します.

STEP1 : gsutilを用いてGCSからデータをダウンロードします
STEP2 : pythonを用いてファイルから構造化データに変換します。 変換対象は、Data Frame形式とCSV種類

実装

環境について

今回のコードはJupyter labなどに貼り付けて利用することを想定しています。
単体でツールとしても使えませんし、ライブラリとしても利用できません。
実装した環境はMacです。

STEP1 gsutilを用いたデータコピー

GCPにはpython用のツールもありますが、アカウント設定等が煩雑な為、コマンドラインツールであるgsutilを用いて、フォルダの中にあるファイルをまとめてローカルにダウンロードします。

dowload.sh

gsutil cp -r gs://BACKET_NAME/realtime/unobus.co.jp/vehicle_position/2020/20200801/ ~/dev/unobus/

解説

• gsutil cp はUNIXのcpコマンドと似ています。
• -r 再帰的に処理するオプション
• gs://BACKET_NAME/realtime/un.... コピー元です。今回はダミーの文字列
• ~/dev/unobus/ コピー先です

一般的に、gsutil -m cp とするとマルチスレッドで動作し高速化できます。しかし、今回のバケットのアクセス権限の問題なのか上手く動作しませんでした。

STEP2 ファイルの読み込みと構造化

データフレームに変換

pq2df.py

from google.transit import gtfs_realtime_pb2
import pandas as pd
import numpy as np
import os
import datetime
import time
path = '/Users/USER_NAME/dev/unobus/20200801'
files = os.listdir(path)
feed = gtfs_realtime_pb2.FeedMessage()

start = time.time()
i = 0;
temp_dict = {}
for file in files:
    with open(path+'/'+file, 'rb') as f:
        data = f.read()
        feed.ParseFromString(data)
        for entity in feed.entity:
            temp_dict[i] = [
                  entity.id,                         #車両ID
                  entity.vehicle.vehicle.id,         #車両ナンバー
                  entity.vehicle.trip.trip_id,       #路線番号？
                  entity.vehicle.timestamp,          #車両時刻
                  entity.vehicle.position.longitude, #車両緯度
                  entity.vehicle.position.latitude,  #車両経度
                  entity.vehicle.occupancy_status #混雑度
            ]
            i +=1
df = pd.DataFrame.from_dict(temp_dict, orient='index',columns=['id' , 'vehicle_id', 'trip_id','vehicle_timestamp','longitude','latitude','occupancy_status'])
elapsed_time = time.time() - start
print ("elapsed_time:{0}".format(elapsed_time) + "[sec]")

CSVへ出力

pq2csv.py

from google.transit import gtfs_realtime_pb2
import pandas as pd
import numpy as np
import os
import datetime
import time
path = '/Users/USER_NAME/dev/unobus/20200801'
csvfilename = 'unobus_20200801.csv'
files = os.listdir(path)
feed = gtfs_realtime_pb2.FeedMessage()
with open(csvfilename, 'a') as csv :
    start = time.time()
    for file in files:
        with open(path+'/'+file, 'rb') as f:
            data = f.read()
            feed.ParseFromString(data)
            for entity in feed.entity:
                print(
                    entity.id,                         #車両ID
                    entity.vehicle.vehicle.id,         #車両ナンバー
                    entity.vehicle.trip.trip_id,       #路線番号？
                    entity.vehicle.timestamp,          #車両時刻
                    entity.vehicle.position.longitude, #車両緯度
                    entity.vehicle.position.latitude,  #車両経度
                    entity.vehicle.occupancy_status,   #混雑度
                    sep=',',file=csv)
    elapsed_time = time.time() - start
    print ("elapsed_time:{0}".format(elapsed_time) + "[sec]")

失敗作 データフレームに変換（超低速）

pq2df.pyと比較すると、220倍も遅いコードです。
df.appendがその原因です。

pq2df_VeryLowSpeed.py

from google.transit import gtfs_realtime_pb2
import pandas as pd
import numpy as np
import os
import datetime
import time
path = '/Users/USER_NAME/dev/unobus/20200801'
files = os.listdir(path)
feed = gtfs_realtime_pb2.FeedMessage()

df = pd.DataFrame(columns=['id' , 'vehicle_id', 'trip_id','vehicle_timestamp','longitude','latitude','occupancy_status'])

start = time.time()
for file in files:
    with open(path+'/'+file, 'rb') as f:
        data = f.read()
        feed.ParseFromString(data)
        for entity in feed.entity:
            tmp_se = pd.Series( [
                  entity.id,                         #車両ID
                  entity.vehicle.vehicle.id,         #車両ナンバー
                  entity.vehicle.trip.trip_id,       #路線番号？
                  entity.vehicle.timestamp,          #車両時刻
                  entity.vehicle.position.longitude, #車両緯度
                  entity.vehicle.position.latitude,  #車両経度
                  entity.vehicle.occupancy_status #混雑度
            ], index=df.columns )
            df = df.append( tmp_se, ignore_index=True ) #ここがだめ！！

elapsed_time = time.time() - start
print ("elapsed_time:{0}".format(elapsed_time) + "[sec]")

終わりに

クラウド上に保存されているGTFSの過去データを、ダウンロードから構造化データへの変換までを実施しました。
DataFrameでデータを追記する処理は書き方によっては非常に処理が遅くなり、苦労しました。

微力ではありますが、バスデータ活用に貢献できれば幸いです。

偉大な先人

既に伊藤先生が、GTFSファイルをマージしてPostgreSQL+PostGISに投入する方法にて、ツールを活用してPostgreSQLにデータを投入する方法を掲載しております。

今回のモチベーション

日本のGTFSデータが、GCSに保存されているとの情報を入手したので、そのデータで遊ぶための準備です。

今回やること

Google Clound Platform のファイルサービスである、Google Cloud Strage 保存されている protocol buffers 形式の大量のファイルを、一気にダウンロードして一纏めのデータに変換します.

STEP1 : gsutilを用いてGCSからデータをダウンロードします
STEP2 : pythonを用いてファイルから構造化データに変換します。 変換対象は、Data Frame形式とCSVの２種類

実装

環境について

今回のコードはJupyter labなどに貼り付けて利用することを想定しています。
単体でツールとしても使えませんし、ライブラリとしても利用できません。
実装した環境はMacです。

STEP1 gsutilを用いたデータコピー

GCPにはpython用のツールもありますが、アカウント設定等が煩雑な為、コマンドラインツールであるgsutilを用いて、フォルダの中にあるファイルをまとめてローカルにダウンロードします。

dowload.sh

gsutil cp -r gs://BACKET_NAME/realtime/unobus.co.jp/vehicle_position/2020/20200801/ ~/dev/unobus/

解説

• gsutil cp はUNIXのcpコマンドと似ています。
• -r 再帰的に処理するオプション
• gs://BACKET_NAME/realtime/un.... コピー元です。今回はダミーの文字列
• ~/dev/unobus/ コピー先です

一般的に、gsutil -m cp とするとマルチスレッドで動作し高速化できます。しかし、今回のバケットのアクセス権限の問題なのか上手く動作しませんでした。

STEP2 ファイルの読み込みと構造化

データフレームに変換

pq2df.py

from google.transit import gtfs_realtime_pb2
import pandas as pd
import numpy as np
import os
import datetime
import time
path = '/Users/USER_NAME/dev/unobus/20200801'
files = os.listdir(path)
feed = gtfs_realtime_pb2.FeedMessage()

start = time.time()
i = 0;
temp_dict = {}
for file in files:
    with open(path+'/'+file, 'rb') as f:
        data = f.read()
        feed.ParseFromString(data)
        for entity in feed.entity:
            temp_dict[i] = [
                  entity.id,                         #車両ID
                  entity.vehicle.vehicle.id,         #車両ナンバー
                  entity.vehicle.trip.trip_id,       #路線番号？
                  entity.vehicle.timestamp,          #車両時刻
                  entity.vehicle.position.longitude, #車両緯度
                  entity.vehicle.position.latitude,  #車両経度
                  entity.vehicle.occupancy_status #混雑度
            ]
            i +=1
df = pd.DataFrame.from_dict(temp_dict, orient='index',columns=['id' , 'vehicle_id', 'trip_id','vehicle_timestamp','longitude','latitude','occupancy_status'])
elapsed_time = time.time() - start
print ("elapsed_time:{0}".format(elapsed_time) + "[sec]")

CSVへ出力

pq2csv.py

from google.transit import gtfs_realtime_pb2
import pandas as pd
import numpy as np
import os
import datetime
import time
path = '/Users/USER_NAME/dev/unobus/20200801'
csvfilename = 'unobus_20200801.csv'
files = os.listdir(path)
feed = gtfs_realtime_pb2.FeedMessage()
with open(csvfilename, 'a') as csv :
    start = time.time()
    for file in files:
        with open(path+'/'+file, 'rb') as f:
            data = f.read()
            feed.ParseFromString(data)
            for entity in feed.entity:
                print(
                    entity.id,                         #車両ID
                    entity.vehicle.vehicle.id,         #車両ナンバー
                    entity.vehicle.trip.trip_id,       #路線番号？
                    entity.vehicle.timestamp,          #車両時刻
                    entity.vehicle.position.longitude, #車両緯度
                    entity.vehicle.position.latitude,  #車両経度
                    entity.vehicle.occupancy_status,   #混雑度
                    sep=',',file=csv)
    elapsed_time = time.time() - start
    print ("elapsed_time:{0}".format(elapsed_time) + "[sec]")

失敗作 データフレームに変換（超低速）

pq2df.pyと比較すると、220倍も遅いコードです。
df.appendがその原因です。

pq2df_VeryLowSpeed.py

from google.transit import gtfs_realtime_pb2
import pandas as pd
import numpy as np
import os
import datetime
import time
path = '/Users/USER_NAME/dev/unobus/20200801'
files = os.listdir(path)
feed = gtfs_realtime_pb2.FeedMessage()

df = pd.DataFrame(columns=['id' , 'vehicle_id', 'trip_id','vehicle_timestamp','longitude','latitude','occupancy_status'])

start = time.time()
for file in files:
    with open(path+'/'+file, 'rb') as f:
        data = f.read()
        feed.ParseFromString(data)
        for entity in feed.entity:
            tmp_se = pd.Series( [
                  entity.id,                         #車両ID
                  entity.vehicle.vehicle.id,         #車両ナンバー
                  entity.vehicle.trip.trip_id,       #路線番号？
                  entity.vehicle.timestamp,          #車両時刻
                  entity.vehicle.position.longitude, #車両緯度
                  entity.vehicle.position.latitude,  #車両経度
                  entity.vehicle.occupancy_status #混雑度
            ], index=df.columns )
            df = df.append( tmp_se, ignore_index=True ) #ここがだめ！！

elapsed_time = time.time() - start
print ("elapsed_time:{0}".format(elapsed_time) + "[sec]")

終わりに

クラウド上に保存されているGTFSの過去データを、ダウンロードから構造化データへの変換までを実施しました。
DataFrameでデータを追記する処理は書き方によっては非常に処理が遅くなり、苦労しました。

微力ではありますが、バスデータ活用に貢献できれば幸いです。

はじめに

注意点などはこちらを参照してください。

この記事を描こうと思ったきっかけ

Pythonの計算に関して、適当な数値を入れてみてみたい。
そこでランダムな数値を精製してくれるrandomモジュールを使用してみることにしました。
そんなrandomモジュールについての備忘録のようなものです。
（更新予定の途中記事です。）

ようこそ、乱数の世界へ

randomモジュールをインポート

randomモジュールはPythonの標準ライブラリに含まれているため、ダウンロードは入りません。
ただし、インポートは必要です。

import random # 上部の段にこれを記述

randomモジュールの関数

random.random()

0.0以上1.0未満の浮動小数点数float型の乱数を生成します。

import random

print(random.random())

# 実行結果↓
# 0.830963714448
# 0.253643551299
# 0.906046050863 など

random.uniform(a,b)

任意の範囲（a から b の間）の浮動小数点数float型の乱数を生成する。

import random

print(random.uniform(10,20))

# 実行結果↓
# 10.4740789672
# 14.8562079851
# 10.2516560825 など

・2つの数値は大きさ順不同。

import random

print(random.uniform(1,-2))

# 実行結果↓
# -14.6138465193
# -1.29546679922
# 5.63021229752 など

・もしも同じ数だったら、その値のみ返す。

import random

print(random.uniform(10,10))

# 実行結果↓
# 10.0

引数は浮動小数点数float型でも可能。
※ bの値が範囲に含まれるかはドキュメント引用部分より、a + (b-a) * random.random()の丸めに依存。

端点の値 b が範囲に含まれるかどうかは、等式 a + (b-a) * random() における浮動小数点の丸めに依存します。 ー Pythonドキュメント：random --- 擬似乱数を生成する より引用

import random

print(random.uniform(1.234,4.321))

# 実行結果↓
# 1.93441022683
# 2.75716839399
# 3.91940634098 など

参考

https://docs.python.org/ja/3/library/random.html#functions-for-integers

今日から、「東京大学のデータサイエンティスト育成講座」を読んで、ちょっと疑問を持ったところや有用だと感じた部分をまとめて行こうと思う。
したがって、あらすじはまんまになると思うが内容は本書とは関係ないと思って読んでいただきたい。
また、環境は問題が発生するまで、Python3.6で実施する。
Python 3.6.10 |Anaconda, Inc.| (default, Jan 7 2020, 15:18:16) [MSC v.1916 64 bit (AMD64)] on win32

Chapter1-1 データサイエンティストの仕事

図1-1-1は、．．．会計士というのが入っているけど、何だろう？？
ドメイン知識とかでいいような気がした。

Chapter1-2 Pythonの基礎

Jupyter Notebookの使い方から、入るけど割愛します。

1-2-2 Pythonの基礎

1-2-2-1　変数

以下のコードは、

msg = 'test'
print(msg)
print(msg[0])
print(msg[1])
print(msg[5])

結果

test
t
e
Traceback (most recent call last):
  File "msg-ex.py", line 5, in <module>
    print(msg[5])
IndexError: string index out of range

errorが出たらググれとあるので、ググってみます。
そのものずばりの回答もありますが、今回はよく使うStackOverflowのサイトのものを貼っておきます。
この程度のエラーでも、躓くことはありそうです。
IndexError: string index out of range in Python [closed]
回答は以下のとおり

As other people have indicated, s[s.Length] isn't actually a valid index; indices are in the closed interval [0, length - 1] (i.e. the last valid index is length - 1 and the first index is 0). Note that this isn't true for every language (there are languages where the first index is 1), but it's certainly true for Python.

Google翻訳してみると以下のとおりです。

他の人が示したように、s [s.Length]は実際には有効なインデックスではありません。インデックスは閉じた間隔[0、length-1]にあります（つまり、最後の有効なインデックスはlength-1で、最初のインデックスは0です）。これはすべての言語に当てはまるわけではないことに注意してください（最初のインデックスが1である言語があります）が、Pythonには確かに当てはまります。

1-2-2-2　演算

data=1
print(data)
data = data+10
print(data)

結果

1
11

ここで大切なお話は、「変数名は大事」というメッセージが大切
以下で確認しておくといいと思います。
Python命名規則一覧
コーディング規約は、一度読んでおくといいと思います。
pythonの標準的なコーディング規約（参考URLの[B-4];p413）
PEP: 8 Python の標準ライブラリに含まれているPythonコードのコーディング規約
PEP 命名規約

1-2-2-3　予約語

予約語と組み込み関数名は、変数名として使えない。
Python3.7予約語またはkeyword一覧と組み込み関数取得・確認方法

1-2-3　リストと辞書型

data_list = [1,2,3,4,5,6,7,8,9]
print(data_list)
print(type(data_list))
print(data_list[1])
print(len(data_list))

結果

[1, 2, 3, 4, 5, 6, 7, 8, 9]
<class 'list'>
2
9

listは、以下のような変更ができる。

data_list = [1,2,3,4,5,6,7,8,9]
print(data_list*2)
data_list.append(4)
print(data_list)
data_list.remove(5)
print(data_list)
print(data_list.pop(6))
print(data_list)
del data_list[3:6]
print(data_list)

結果

[1, 2, 3, 4, 5, 6, 7, 8, 9, 1, 2, 3, 4, 5, 6, 7, 8, 9]
[1, 2, 3, 4, 5, 6, 7, 8, 9, 4]
[1, 2, 3, 4, 6, 7, 8, 9, 4]
8
[1, 2, 3, 4, 6, 7, 9, 4]
[1, 2, 3, 9, 4]

【参考】
Pythonでリスト（配列）の要素を削除するclear, pop, remove, del

tupleは、listに似ているが変更できない。

data_tuple = (1,2,3,4,5,6,7,8,9)
print(data_tuple)
print(type(data_tuple))
print(data_tuple[1])
print(len(data_tuple))

結果

(1, 2, 3, 4, 5, 6, 7, 8, 9)
<class 'tuple'>
2
9

setは集合であり、要素の重複を許さない

data_set = {1,2,3,4,5,6,7,8,9}
print(data_set)
print(type(data_set))
#print(data_set[1])
print(len(data_set))

結果
コメントアウトの行
TypeError: 'set' object does not support indexing
setはインデックスをサポートしていない

{1, 2, 3, 4, 5, 6, 7, 8, 9}
<class 'set'>
9

class 'list' 'tuple''set'の使い方

listは、自由度高く、要素を追加削除できる。
tupleは、自由度無く、変更できない
setは、重複しない要素（同じ値ではない要素、ユニークな要素）のコレクションで、和集合、積集合、差集合などの集合演算を行うことができる。
要素の重複を許さないので、以下のような変更でユニークな要素に出来る。
組み合わせると、以下のようなことが出来る。

data_list = [1,2,3,4,5,6,7,2,3,4,8,9]
print(data_list)
data_tuple = tuple(data_list)
print(data_tuple)
data_set = set(data_list)
print(data_set)
data_list_new = list(data_set)
print(data_list_new)
data_tuple_new = tuple(data_set)
print(data_tuple_new)
data_set.discard(6)
print(data_set)
data_set.add(11)
print(data_set)

結果

[1, 2, 3, 4, 5, 6, 7, 2, 3, 4, 8, 9]
(1, 2, 3, 4, 5, 6, 7, 2, 3, 4, 8, 9)
{1, 2, 3, 4, 5, 6, 7, 8, 9}
[1, 2, 3, 4, 5, 6, 7, 8, 9]
(1, 2, 3, 4, 5, 6, 7, 8, 9)
{1, 2, 3, 4, 5, 7, 8, 9}
{1, 2, 3, 4, 5, 7, 8, 9, 11}

【参考】
Pythonのsetを知る
Python, set型で集合演算（和集合、積集合や部分集合の判定など）

辞書型

dict_data = {'apple':100,'banana':100,'orange':300,'mango':400,'melon':500}
print(dict_data['apple'])
print(dict_data['apple']+dict_data['orange'])

結果

100
400

辞書型追加

dict_data['りんご'] = 100
print(dict_data)

結果

{'apple': 100, 'banana': 100, 'orange': 300, 'mango': 400, 'melon': 500, 'りんご': 100}

辞書型削除

removed_value = dict_data.pop('banana')
print(removed_value)
print(dict_data)

結果

100
{'apple': 100, 'orange': 300, 'mango': 400, 'melon': 500, 'りんご': 100}

もう一つの削除

del dict_data['mango']
print(dict_data)

結果

{'apple': 100, 'orange': 300, 'melon': 500, 'りんご': 100}

辞書同士を連結（結合）、複数の要素を追加・更新: update()

連結（結合）

d1 = {'apple': 100, 'banana': 100}
d2 = {'orange': 300, 'mango': 400}
d3 = {'melon': 500, 'りんご': 100}
d1.update(d2)
print(d1)
d1.update(**d2, **d3)
print(d1)

結果

{'apple': 100, 'banana': 100, 'orange': 300, 'mango': 400}
{'apple': 100, 'banana': 100, 'orange': 300, 'mango': 400, 'melon': 500, 'りんご': 100}

複数の要素を追加・更新

d={'apple': 100, 'banana': 100}
print(d)
d.update([('orange', 300), ('mango', 400), ('melon', 500), ('りんご', 100)])
print(d)

結果

{'apple': 100, 'banana': 100}
{'apple': 100, 'banana': 100, 'orange': 300, 'mango': 400, 'melon': 500, 'りんご': 100}

【参考】
Pythonで辞書に要素を追加、辞書同士を連結（結合）

まとめ

・変数
・演算
・予約語と組み込み関数
・リストと集合、そして辞書型
を見た

基本に立ち返って、並べてみると分かり易い。

Optunaを使ったxgboostの設定方法

xgboostの回帰について設定してみる。

xgboostについては、他のHPを参考にしましょう。
「ザックリとした『Xgboostとは』& 主要なパラメータについてのメモ」
https://qiita.com/2357gi/items/913af8b813b069617aad

後、公式HPのパラメーターのところを参考にしました。
https://xgboost.readthedocs.io/en/latest/parameter.html

いろいろ入れたけど、決定木系は過学習になりやすいので、それを制御するパラメーターをしっかり設定した方が良いと思ってます。
xgboostでは、lambdaとalphaですが、pythonで設定するときは、reg_lambdaとreg_alphaのように、reg_をつけて指定します。

# optunaの目的関数を設定する
#gtreeのパラメーター設定です。
def objective(trial):
    eta =  trial.suggest_loguniform('eta', 1e-8, 1.0)
    gamma = trial.suggest_loguniform('gamma', 1e-8, 1.0)
    max_depth = trial.suggest_int('max_depth', 1, 20)
    min_child_weight = trial.suggest_loguniform('min_child_weight', 1e-8, 1.0)
    max_delta_step = trial.suggest_loguniform('max_delta_step', 1e-8, 1.0)
    subsample = trial.suggest_uniform('subsample', 0.0, 1.0)
    reg_lambda = trial.suggest_uniform('reg_lambda', 0.0, 1000.0)
    reg_alpha = trial.suggest_uniform('reg_alpha', 0.0, 1000.0)


    regr =xgb.XGBRegressor(eta = eta, gamma = gamma, max_depth = max_depth,
                           min_child_weight = min_child_weight, max_delta_step = max_delta_step,
                           subsample = subsample,reg_lambda = reg_lambda,reg_alpha = reg_alpha)

    score = cross_val_score(regr, X_train, y_train, cv=5, scoring="r2")
    r2_mean = score.mean()
    print(r2_mean)

    return r2_mean

# optunaで最適値を見つける
study = optuna.create_study(direction='maximize')
study.optimize(objective, n_trials=500)

# チューニングしたハイパーパラメーターをフィット
optimised_rf = xgb.XGBRegressor(eta = study.best_params['eta'],gamma = study.best_params['gamma'],
                                max_depth = study.best_params['max_depth'],min_child_weight = study.best_params['min_child_weight'],
                                max_delta_step = study.best_params['max_delta_step'],subsample = study.best_params['subsample'],
                                reg_lambda = study.best_params['reg_lambda'],reg_alpha = study.best_params['reg_alpha'])

optimised_rf.fit(X_train ,y_train)

ボストンの結果です。

全部です。

# -*- coding: utf-8 -*-

from sklearn import datasets
import xgboost as xgb
from sklearn.model_selection import train_test_split
from sklearn.metrics import r2_score
import pandas as pd
import optuna
import numpy as np
import matplotlib.pyplot as plt
from sklearn.model_selection import cross_val_score

#ボストンのデータセットを読み込む 
boston = datasets.load_boston()

#print(boston['feature_names'])
#特徴量と目的変数をわける
X = boston['data']
y = boston['target']
X_train, X_test, y_train, y_test = train_test_split(X, y, train_size=0.8)

# optunaの目的関数を設定する
#gtreeのパラメーター設定です。
def objective(trial):
    eta =  trial.suggest_loguniform('eta', 1e-8, 1.0)
    gamma = trial.suggest_loguniform('gamma', 1e-8, 1.0)
    max_depth = trial.suggest_int('max_depth', 1, 20)
    min_child_weight = trial.suggest_loguniform('min_child_weight', 1e-8, 1.0)
    max_delta_step = trial.suggest_loguniform('max_delta_step', 1e-8, 1.0)
    subsample = trial.suggest_uniform('subsample', 0.0, 1.0)
    reg_lambda = trial.suggest_uniform('reg_lambda', 0.0, 1000.0)
    reg_alpha = trial.suggest_uniform('reg_alpha', 0.0, 1000.0)


    regr =xgb.XGBRegressor(eta = eta, gamma = gamma, max_depth = max_depth,
                           min_child_weight = min_child_weight, max_delta_step = max_delta_step,
                           subsample = subsample,reg_lambda = reg_lambda,reg_alpha = reg_alpha)

    score = cross_val_score(regr, X_train, y_train, cv=5, scoring="r2")
    r2_mean = score.mean()
    print(r2_mean)

    return r2_mean

# optunaで最適値を見つける
study = optuna.create_study(direction='maximize')
study.optimize(objective, n_trials=500)

# チューニングしたハイパーパラメーターをフィット
optimised_rf = xgb.XGBRegressor(eta = study.best_params['eta'],gamma = study.best_params['gamma'],
                                max_depth = study.best_params['max_depth'],min_child_weight = study.best_params['min_child_weight'],
                                max_delta_step = study.best_params['max_delta_step'],subsample = study.best_params['subsample'],
                                reg_lambda = study.best_params['reg_lambda'],reg_alpha = study.best_params['reg_alpha'])

optimised_rf.fit(X_train ,y_train)
#結果の表示
print("訓練データにフィット")
print("訓練データの精度 =", optimised_rf.score(X_train, y_train))
pre_train = optimised_rf.predict(X_train)
print("テストデータにフィット")
print("テストデータの精度 =", optimised_rf.score(X_test, y_test))
pre_test = optimised_rf.predict(X_test)

#グラフの表示
plt.scatter(y_train, pre_train, marker='o', cmap = "Blue", label="train")
plt.scatter(y_test ,pre_test, marker='o', cmap= "Red", label="test")
plt.title('boston')
plt.xlabel('measurment')
plt.ylabel('predict')
#ここでテキストは微調整する
x = 30  
y1 = 12
y2 = 10
s1 =  "train_r2 =" + str(optimised_rf.score(X_train, y_train))
s2 =  "test_r2 =" + str(optimised_rf.score(X_test, y_test))
plt.text(x, y1, s1)
plt.text(x, y2, s2)

plt.legend(loc="upper left", fontsize=14)
plt.show()

はじめに

これも自分用のメモです。python を使うと楽にできます。注意点として、

• 緯度経度は度分秒なので６０進数
• 値は実数ではなく、分数で。（有理数、というのかな）

いくつか関係するライブラリがあります。どれか一つに統一したほうがよいのですが、とりあえず動いたものメモです。（おい）

仕様書?

下記の文献を見つけました。

• デジタルスチルカメラ用 画像ファイルフォーマット規格 Exif 2.3

これによるとGPS Info IFDは以下のように定義されているようです。

たくさんありますが、私が扱ったファイルでは0 - 6 までのタグのみでした。GPSの緯度、経度、高さが読み書きできたので、とりあえずは十分でした。ここでGPSLatitudeRefは一文字の'N'か'S'で北緯、南緯を表します。GPSLatitudeは分数で表します。

python で読み書きする

Pillow で緯度経度を読む

Pillow を使用しました。準備で苦労したかは覚えていません。

pip3 install pillow

下記で緯度経度を読めました。

read_gpsexif.py

from PIL import Image
infile = "test.JPG"
img = Image.open( infile )
exif = img._getexif()
for k, v in exif.items():
    if k==34853:
        print(v)

これで

{0: b'\x02\x00\x00\x00',
 1: 'N',
 2: (35.0, 20.0, 53.51123),
 3: 'E',
 4: (137.0, 9.0, 17.83123),
 5: b'\x00', 6: 256.123}

のような結果を得ます。先の仕様書のGPS INfo IFD の定義と見比べると意味が分かると思います。度分秒なので、１０進法にするには変換が必要です。例えば

deg,minu,sec = 35.0, 20.0, 53.51366
deg + minu/60.0 + sec/3600.0

で 35.34819823888889 を得ます。

piexif でGPS Exif の緯度経度を書き換える

piexifを利用しました。
度分秒と分数への変換を関数にし、Gps exif を読んで書き換える作戦です。

def _to_deg(value, loc):                                                                                                                                               
    if value < 0:                                                                                                                                                    
       loc_value = loc[0]                                                                                                                                       
    elif value > 0:                                                                                                                                                  
       loc_value = loc[1]
    else:                                                                                                                                                            
       loc_value = ""                                                                                                                                           
    abs_value = abs(value)
    deg =  int(abs_value)                                                                                                                                        
    t1 = (abs_value-deg)*60                                                                                                                                      
    minu = int(t1)                                                                                                                                                
    sec = round((t1 - minu)* 60, 5)
    return (deg, min, sec, loc_value)                                                                                                                                                                                                                                                                                     

import Fraction
def _change_to_rational(number):
    f = Fraction(str(number))                                                                                                                                    
    return (f.numerator, f.denominator)     

を用いて、以下を実行します。

write_gpsexif.py

import piexif
lat, lng, altitude = 39.123, 139.123, 50.123
jpg_filepath = "test.JPG"

lat_deg = to_deg(lat, ["S", "N"])
lng_deg = to_deg(lng, ["W", "E"])
                                                                                                                            exif_lat = (change_to_rational(lat_deg[0]), change_to_rational(lat_deg[1]), change_to_rational(lat_deg[2]))
exif_lng = (change_to_rational(lng_deg[0]), change_to_rational(lng_deg[1]), change_to_rational(lng_deg[2]))
gps_ifd = {
    piexif.GPSIFD.GPSVersionID: (2, 0, 0, 0),
    piexif.GPSIFD.GPSAltitudeRef: 0,
    piexif.GPSIFD.GPSAltitude: change_to_rational(round(altitude, 3)),
    piexif.GPSIFD.GPSLatitudeRef: lat_deg[3],
    piexif.GPSIFD.GPSLatitude: exif_lat,
    piexif.GPSIFD.GPSLongitudeRef: lng_deg[3],
    piexif.GPSIFD.GPSLongitude: exif_lng,
}

gps_exif = {"GPS": gps_ifd}
exif_data = piexif.load(jpg_filepath)
exif_data.update(gps_exif)
exif_bytes = piexif.dump(exif_data)
piexif.insert(exif_bytes, file_name)

これは確かStackOverflow を参考にさせていただいたと思うのですが、linkを無くしてしまいました。見つけたら追加します。

まとめ・雑感

必要にせまらて写真のGPSExif の読み書きをした時のメモです。

• Pillow, piexif を利用したが、何となくどちらかに統一して正しく使うのがよさそう。

えいやの作業メモなので、拡充させていきたい、、、という投稿でした。最近雑だな。。^^; (2020/08/08)

これは何度も自分でググっているので、自分用のメモです。

dtype 変更

まとめて変更したいとき

df = df.astype({'a': 'float32', 'c': 'int8'})

備考

pandas 関係でいつも何度も調べなおしていることがあり、思い出したらここに追加していこう。(2020/08/08)

環境

• Windows 10
• Python 3.8

プログラム

↓を参考にさせていただきました。
https://cercopes-z.com/Python/tips-mail-py.html#yahoo

しかし、なぜか↓のエラーが出た。
smtplib.SMTPServerDisconnected: Connection unexpectedly closed

Yahooメール側で外部からのメールを受信できる設定にする必要があった。
https://knowledge.support.yahoo-net.jp/PccMail/s/article/H000007321

Pythonの多次元リストのコピーについて、その挙動が自分にとっては理解しづらいので実験結果とともに整理。
結論から言うと、copy()、リスト[:]、copyライブラリのcopy.copy()は浅いコピー(shallow copy)であり、新たな多次元リストを作成し、その後元のリスト中に見つかったリストに対する参照を挿入する。
copyライブラリのcopy.deepcopy()は深いコピー(deep copy)であり、新たな多次元リストを作成し、その後元のリスト中に見つかったリストのコピーを挿入する。

copy()

「1層目」というのがいいのか、「1次元目」というのがいいのかよくわからないが…。コピーによって別の多次元リストができるので、IDは当然別のものになる。しかしそのリストの中身は参照なので、コピーの1層目を書き換えてもオリジナルのリストにその影響はないが、2層目以降の書き換えはIDが違うのにもかかわらずオリジナルにも波及する。

copy()

lst = [[0]*3 for _ in range(3)]

lst1 = lst.copy()

lst1[1] = ['r'] * 3
print('lst1 = ', lst1)
print('id(lst1) = %d' % id(lst1))
print('lst = ', lst)
print('id(lst) = %d' % id(lst))

lst1[0][1] = 'r'
print('lst1 = ', lst1)
print('id(lst1) = %d' % id(lst1))
print('lst = ', lst)
print('id(lst) = %d' % id(lst))

出力

lst1 =  [[0, 0, 0], ['r', 'r', 'r'], [0, 0, 0]]
id(lst1) = 139730947400512
lst =  [[0, 0, 0], [0, 0, 0], [0, 0, 0]]
id(lst) = 139730946890944
lst1 =  [[0, 'r', 0], ['r', 'r', 'r'], [0, 0, 0]]
id(lst1) = 139730947400512
lst =  [[0, 'r', 0], [0, 0, 0], [0, 0, 0]]
id(lst) = 139730946890944

リスト[:]

copy()と同様。

リスト[

lst = [[0]*3 for _ in range(3)]

lst1 = lst[:]

lst1[1] = ['r'] * 3
print('lst1 = ', lst1)
print('id(lst1) = %d' % id(lst1))
print('lst = ', lst)
print('id(lst) = %d' % id(lst))

lst1[0][1] = 'r'
print('lst1 = ', lst1)
print('id(lst1) = %d' % id(lst1))
print('lst = ', lst)
print('id(lst) = %d' % id(lst))

出力

lst1 =  [[0, 0, 0], ['r', 'r', 'r'], [0, 0, 0]]
id(lst1) = 140313735131456
lst =  [[0, 0, 0], [0, 0, 0], [0, 0, 0]]
id(lst) = 140313734622336
lst1 =  [[0, 'r', 0], ['r', 'r', 'r'], [0, 0, 0]]
id(lst1) = 140313735131456
lst =  [[0, 'r', 0], [0, 0, 0], [0, 0, 0]]
id(lst) = 140313734622336

copyライブラリcopy.copy()

copy()と同様。

copy.copy()

import copy

lst = [[0]*3 for _ in range(3)]

lst1 = copy.copy(lst)

lst1[1] = ['r'] * 3
print('lst1 = ', lst1)
print('id(lst1) = %d' % id(lst1))
print('lst = ', lst)
print('id(lst) = %d' % id(lst))

lst1[0][1] = 'r'
print('lst1 = ', lst1)
print('id(lst1) = %d' % id(lst1))
print('lst = ', lst)
print('id(lst) = %d' % id(lst))

出力

lst1 =  [[0, 0, 0], ['r', 'r', 'r'], [0, 0, 0]]
id(lst1) = 139783910001664
lst =  [[0, 0, 0], [0, 0, 0], [0, 0, 0]]
id(lst) = 139783910057280
lst1 =  [[0, 'r', 0], ['r', 'r', 'r'], [0, 0, 0]]
id(lst1) = 139783910001664
lst =  [[0, 'r', 0], [0, 0, 0], [0, 0, 0]]
id(lst) = 139783910057280

copyライブラリcopy.deepcopy()

コピーによって別の多次元リストができるのでIDは別のものになり、そのリストの中身もコピーされる。いわば、クローンが作られる。全く同じ内容だが別のものなので、それを書き換えてもオリジナルには影響がない。

copy.deepcopy()

import copy

lst = [[0]*3 for _ in range(3)]

lst1 = copy.deepcopy(lst)

lst1[1] = ['r'] * 3
print('lst1 = ', lst1)
print('id(lst1) = %d' % id(lst1))
print('lst = ', lst)
print('id(lst) = %d' % id(lst))

lst1[0][1] = 'r'
print('lst1 = ', lst1)
print('id(lst1) = %d' % id(lst1))
print('lst = ', lst)
print('id(lst) = %d' % id(lst))

出力

lst1 =  [[0, 0, 0], ['r', 'r', 'r'], [0, 0, 0]]
id(lst1) = 139629482679488
lst =  [[0, 0, 0], [0, 0, 0], [0, 0, 0]]
id(lst) = 139629482735040
lst1 =  [[0, 'r', 0], ['r', 'r', 'r'], [0, 0, 0]]
id(lst1) = 139629482679488
lst =  [[0, 0, 0], [0, 0, 0], [0, 0, 0]]
id(lst) = 139629482735040

参考

Python - リストを複製する
copy --- 浅いコピーおよび深いコピー操作
【Python】 自分で作った関数の引数にリストを代入すると、そのオリジナルのリストが改変されてしまうのを防ぐ方法

はじめに

測量では成果を緯度経度でなく、平面直角座標系という２次元の平面の座標値で表現します。その相互変換を行う方法のメモです。

日本では現在、19の直角座標系が用いられています。(「分かりやすい平面直角座標系」from 国土地理院)

よくある(?)誤解: 慣性航法で使う座標系とは違う！

座標系の変換というと簡単に聞こえますが、これは

• WGS84の楕円体の極座標(緯度,経度,楕円体高)を地球中心の座標系(ECEF)の直交座標系の座標値XYZに変換
• ECEFの直交座標系の値を、該当する楕円体における接平面に対応するように回転させる

という話とは『異なり』ます。ガウス・クリューゲル投影法という計算方法を使うようなのですが、私は数式が複雑で着いていけませんでした。

当初、自分は上記のような単純な回転計算と誤解していて、周りに「そんなの簡単な計算だよ」などと話していて、思い返すと恥ずかしいです。

国土地理院のサポート

変換するために、国土地理院から提供されているものを利用できます。（さすが公共機関！）

まず、ブラウザに直接座標値を入力して変換することができます。

https://vldb.gsi.go.jp/sokuchi/surveycalc/surveycalc/bl2xyf.html

また、WEB APIからも利用できます。こちらに、提供されている機能と説明があります。アクセスが集中しないように同一アドレスから10秒間で10回の制限がある等の注意も書かれています。

https://vldb.gsi.go.jp/sokuchi/surveycalc/main.html

平面直角座標系の番号を指定して取り出す例です。

use_gsi_webapi.py

import requests
from numpy import round
lat, lon = 35.71, 139.74

url = "http://vldb.gsi.go.jp/sokuchi/surveycalc/surveycalc/bl2xy.pl?"
params={'latitude':round(lat, 8), 'longitude': round(lon, 8), "refFrame":2, "zone": 9, 'outputType':'json'}

res =　requests.get(url, params=params)
if res.status_code == requests.codes.ok:
    print(res.json())

ここでは次の辞書が返されます。

{'OutputData': {
  'publicX': '-32170.0990',
  'publicY': '-8445.1361',
  'gridConv': '0.054477778',
  'scaleFactor': '0.99990088'}}

EPSG

日本の平面直角座標系の番号は日本が定めるものですが、国際的にはEPSGとして定めらています。国際的には平面直角座標系の番号ではなくEPSGで参照します。(例えば写真測量のソフトであるmetashape などを使う場合とか。)

EPSG	測地系
EPSG:4326	WGS 84
EPSG:6668	JGD2011
EPSG:6669	JGD2011 / Japan Plane Rectangular CS I
EPSG:6670	JGD2011 / Japan Plane Rectangular CS II
EPSG:6671	JGD2011 / Japan Plane Rectangular CS III
EPSG:6672	JGD2011 / Japan Plane Rectangular CS IV
EPSG:6673	JGD2011 / Japan Plane Rectangular CS V
EPSG:6674	JGD2011 / Japan Plane Rectangular CS VI
EPSG:6675	JGD2011 / Japan Plane Rectangular CS VII
EPSG:6676	JGD2011 / Japan Plane Rectangular CS VIII
EPSG:6677	JGD2011 / Japan Plane Rectangular CS IX
EPSG:6678	JGD2011 / Japan Plane Rectangular CS X
EPSG:6679	JGD2011 / Japan Plane Rectangular CS XI
EPSG:6680	JGD2011 / Japan Plane Rectangular CS XII
EPSG:6681	JGD2011 / Japan Plane Rectangular CS XIII
EPSG:6682	JGD2011 / Japan Plane Rectangular CS XIV
EPSG:6683	JGD2011 / Japan Plane Rectangular CS XV
EPSG:6684	JGD2011 / Japan Plane Rectangular CS XVI
EPSG:6685	JGD2011 / Japan Plane Rectangular CS XVII
EPSG:6686	JGD2011 / Japan Plane Rectangular CS XVIII
EPSG:6687	JGD2011 / Japan Plane Rectangular CS XIX

情報源ですが、本家本元や管理団体などよく理解していません。以下がわりとauthority のような気がします。

• epsg.io で説明を読むことができます。例えばWGS84とか。
• EPSG Geodetic Parameter Registryで検索できます
• spatial reference というページにもリストや検索機能があります。

pyproj　を使う

さきほどはWEB API をりようしましたが、これはpython ライブラリを用いてもできます。
インストールで苦労した記憶はありません。

test_webapi.py

from pyproj import Transformer

lat, lon = 35.71, 139.74
wgs84_epsg, rect_epsg = 4326, 6677
tr = Transformer.from_proj(wgs84_epsg, rect_epsg)
x, y = tr.transform(lat, lon)
print(x, y)

実行すると以下を得ます。

-32170.09900022997 -8445.136058616452

有効数字は分かりませんが、地理院の結果とは一致しています。（小数点以下４桁、0.1ミリメートル！）

CRSの利用

巷では、CRSを使う方法が紹介されていましたが、利用してみると、、、

test_use_pyproj.py

import pyproj
wgs84=pyproj.CRS("EPSG:4326")
jgd2011_9 = pyproj.CRS("EPSG:6677")
lat, lon = 35.71, 139.74
print( pyproj.transform(wgs84, jgd2011_9, lat, lon) )

実行すると警告が出ますが、同じ結果を得ます。

/usr/bin/ipython3:1: DeprecationWarning: This function is deprecated. See: https://pyproj4.github.io/pyproj/stable/gotchas.html#upgrading-to-pyproj-2-from-pyproj-1
(-32170.09900022997, -8445.136058616452)

これは多分、古い方法なのだと思います。しっかりチュートリアルを読まなければ！
https://pyproj4.github.io/pyproj/stable/examples.html

まとめ・雑感

• WGS84と平面直角座標系の変換は地理院提供のＡＰＩでできるが、pyproj でも同じことはできそう
• pyproj は多機能で4Dtansform など気になる機能もあり、余力があれば調べてみたい

等々(2020/08/08)

Django==2.2.12でのお話です。
ログインした状態で、自作の画面からユーザパスワードを変更したあと

・別のページに遷移しようとすると、まずログイン画面に飛んでしまう
・普段なら問題のないサーバーへのリクエストがエラーを返す

といった状況を不思議に思ったアナタ！

実はカスタムのパスワード変更画面を利用している場合、パスワード変更後に強制的にログアウト状態になってしまうんですね。
それまで使用していたセッションが使われなくなるのが原因です。

Djangoが用意しているPasswordChangeViewなんかを使用している場合は、パスワード変更後に自動で新しいセッションを作成してくれるので、ログイン状態が続いてくれます。

カスタム画面でもそうしたい！という場合はupdate_session_auth_hash()関数を使うといいよ。

といったことが、ドキュメント
https://docs.djangoproject.com/en/3.0/topics/auth/default/#session-invalidation-on-password-change
に書かれていました...

ついでに、このドキュメントに書かれている使用例も貼っておきます。

from django.contrib.auth import update_session_auth_hash

def password_change(request):
    if request.method == 'POST':
        form = PasswordChangeForm(user=request.user, data=request.POST)
        if form.is_valid():
            form.save()
            update_session_auth_hash(request, form.user)
    else:
        ...

何をしたいのか？

• PythonでGrADS用のファイルを読み込んで後処理をしたい
• しかし共有計算資源の環境はクリーンでなければならない
• だがしかし、変換ソフト、ライブラリなどなどインストールしたいものがいっぱいある

ということで、コンテナ仮想化を使おうということです。

Docker

Dockerを使うと何がいいかというと、計算資源を汚すことなく、必要な後処理を行う専用の仮想計算機を作り上げることができます。使うときにコンテナを立ち上げ、処理を行った後は破棄してまた必要になったら立ち上げることができます。
Dockerの場合は立ち上げのコストが非常に低いので、気軽に立ち上げて一つの処理をした後また破棄することが気軽にできます。

環境整備手順

Dockerfileを書くと自動化できていいのですが、ここでは単にシェルを立ち上げて整備して、それを単純にコミットしてイメージを作ることにします。
長期的に使いたいなら、バージョンアップも可能になるようにDockerfileを作るべきですが、最初の基本の機能の確認はシェルで行ったほうが理解が深まります。

今回はUbuntuをもとに、cdoとpip3, netcdf4をインストールします。

1. もとになるイメージをpull(拾ってくる)して走らせる

今回はUbuntuを使います。

$ docker pull ubuntu/ubuntu

$ docker images

で、イメージが作れたことが確認できます。

イメージを走らせるためには、

$ docker run -it ubuntu/ubuntu

とします。

2. 必要なソフトウェアのインストール

pip3とcdoをインストールします。
pip3でnetCDF4もインストールします。

# apt update
# apt install python3 python3-pip cdo
# pip3 install netCDF4

3. 動作確認

きちんとcdoで変換できてnetCDF4もインポートできるか確認します。

# python3
>>> import netCDF4 as nc
>>> ^D
# cdo
....

4. コミットしてimageを作成する

exitコマンドかCtrl+dで抜けて、イメージを作りましょう。
イメージを作るために必要なコンテナIDは、終了したコンテナを確認するために-aが必要です。

# exit
(host) $ docker ps -a

このような画面になると思います。
CONTAINER IDはコミットの時に指定するコンテナIDになります。
IMAGEのところにimageの名前があって、STATUSのところでexited ** agoとなっているので、大体どれかはそこを見ればわかると思います。

$ docker commit [container id] [author]/[image名]

というようにコミットしてあげましょう。

$ docker images

で、イメージが作れたことが確認できます。

5. 実際に計算させてみる

実際に計算させるとき、２通りあって、runしてからインタラクティブにシェルから実行することもできます。ただ、今回はrunの際にコマンドを指定して実行させてみます。

プログラム、データ群のディレクトリをDockerへとバインドすることで、実際にコンテナから容易にアクセスることができます。また、グラフの確認にXフォワーディングをします。

$ sudo docker run -e DISPLAY=$DISPLAY --net host -v /tmp/.X11-unix:/tmp/.X11-unix -v $HOME/.Xauthority:/root/.Xauthority -v /home/xxx/:/home/ --shm-size 16g xxx/ubuntu-cdo2 /bin/bash -c "cd /home/speedy-epyc/speedy/python-script; python3 rmse.py"

--shm-sizeは、足りなくなるプログラムがあれば必要に応じてつけてください。

Flaskモジュールのインストール

(base) root@e8cf64ce12e9:/home/continuumio#pip install flask 

フォルダ構造の作成

Flaskで使用するHtmlファイルは、プロジェクトフォルダに「templates」という名称のフォルダを作成して格納する必要がある。
注：テンプレートという名前だけど、要はFlaskで利用するHtmlファイルはすべてここに格納すれば良い。

.
├── __pycache__
├── hello.py
└── templates
    ├── hello.html
    └── layout.html

余談：上記のようなディレクトリ構造を取得するためには、「tree」コマンドを別途インストールする必要がある。
Docker環境がUbuntuだったので「apt install tree」でインストールできた。Macローカルの場合のインストールは別途作業が必要みたい。

各ファイルの役割は下記の通り
- hello.py: アプリケーションのルーティングファイル。
- template/hello.html: 画面の一部を構成するhtmlファイル。とどのつまり画面を構成するパーツ。
- templates/layout.html: すべてのHtmlファイルのベースとしてのhtmlファイル。このファイルの部分部分にhello.htmlのようなパーツとしてのhtmlフィアルが組み込まれてくイメージ。

各ファイルの内容

layout.html

layout.html

<!DOCTYPE html>
<html>
<head>
    <titile>{{ title }}</titile>
</head>
<body>
    {% block content %}
    <!-- main-->
    {% endblock %}
</body>
</html>

• {{ val }}: val部分には同名の変数の値が反映される。
• {% block content%}~{% endblock %}:このタグ範囲に拡張パーツとしてのhtmlファイルに記載されている{% block content%}~{% endblock %}のスクリプトが差し込まれるイメージ

hello.html

hello.html

<!-- layout.htmlをテンプレートに拡張する-->
{% extends "layout.html" %}
<!-- block content ~ endblock-の範囲がテンプレートの同宣言範囲に差し込まれる -->
{% block content %}
<h3>Hello</h3>
こんにちは。{{ name }}さん。
{% endblock %}

• extends で指定したｈｔｍｌにこのｈｔｍｌ内のblock contentの内容を反映する。継承のイメージ。

hello.py

hello.py

from flask import Flask, render_template

app = Flask(__name__)

@app.route('/')
def hello():
    name = 'TEST'
    # return name
    return render_template('hello.html', tite="flask test", name=name)

if __name__ == "__main__":
    app.run(host='0.0.0.0', port=5000, debug=True)

• 「render_template」をimportすることでJinja2テンプレートの機能が利用できるようになる。
• render_template([呼び出すhtml],[[テンプレートに渡す変数名=値],...]
• 「if name == "main":」　コマンドラインからpythonファイルを実行する際のおまじないという認識でOK。
• app_run([host=0,0,0,0],[port=XXXX],[debug=XXXX])でアプリケーションを起動する。 hostとかポートはここで指定しなくても動かす方法がありそうだが、おいおい調査する。

実行

(base) root@e8cf64ce12e9:/home/continuumio# python hello.py
 * Serving Flask app "hello" (lazy loading)
 * Environment: production
   WARNING: This is a development server. Do not use it in a production deployment.
   Use a production WSGI server instead.
 * Debug mode: on
 * Running on http://0.0.0.0:5000/ (Press CTRL+C to quit)
 * Restarting with inotify reloader
 * Debugger is active!
 * Debugger PIN: 197-268-519

画面が表示できた。

最後に

とりあえず最低限の使い方だけわかったけど、実務でFlask使うことってないだろうし、事業会社はDjangoなんでしょうね。
あくまでも、フレームワークの足がかりとしてなら勉強する価値あるかな。

はじめに

Tech Dive様の記事 を大変参考にさせて頂きました。この場を借りて御礼申し上げます。

この記事は上記内容をWindows環境で実施した際の備忘録になります。

筆者の環境

・Windows 10
・Git Bash

手順

1. Python3 のインストール

2. Pandasのインストール pip install pandas

3. aws cliのインストール
   インストールしてもGit Bash で aws コマンドが使えなかったので、GitBashで
   cmd \\C aws --version
   を叩いてから再起動したら使えるようになりました

4. アクセスキーの確認
   AWS ユーザーメニューの「マイセキュリティ資格情報」から、インポート用のアクセスキーを新規作成

5. アクセスキーとシークレットアクセスキーの登録
   aws configure を叩き、先ほど作ったアクセスキーを登録

   aws configure
   
   AWS Access Key ID [None]: [アクセスキー]
   AWS Secret Access Key [None]: [シークレットアクセスキー]
   Default region name [None]: us-east-2 ※テーブルがあるリージョン
   Default output format [None]: json
   

6. CSVのヘッダーを変更
   先頭行をデータベースのkeyと合わせて、型を()でくくって指定
   例） キーが「UserName」(String)なら UserName (S)

7. import_to_dynamodb をclone
   git clone https://github.com/hidesan-xyz/import_to_dynamodb.git

8. cloneしたディレクトリのルートにインポートするcsvファイルを置く

9. import用shellコマンドの作成
   python create_insert_command.py importdata testtable
   ※第一引数にcsv名(拡張子は記載しない) 第二引数にテーブル名

10. 作成されたシェルコマンドの実行 sh ./testtable_import_20200808164839.sh

上記手順で完了です。

はじめに

写真のJPGファイルに、撮影日時や場所などを記録することができます。ドローン等で撮影された写真にタグ情報が書かれていることがあります。それらは一般的なフォーマットに従って記録されているようです。それらを読み書きしたいのが目的です。

GPS Exif については多くの記述があります。なので、ここでは情報が少ないと思われる xmp についてメモしておきます。なお、読むことはできるのですが、書くところはうまくできていません。(涙)

準備: xmp-toolkit のインストール

XMPとは

どうやら、Adobe が定義し、各ファイルにXML形式のテキストで情報をファイルのヘッダに加えるために用いていたようです。詳しくは分かりません。^^;)

なので、例えば、less コマンドで見ると、テキストのXMLの部分はそのまま読めます。
おじさんはLinux環境(WSL)ですが、
- less でファイルを開き
- 「binaryだけど開くか」と聞かれるのでYes とし
- xml という文字列を検索する(スラッシュ xml リターン）と、

<x:xmpmeta xmlns:x="adobe:ns:meta/">
 <rdf:RDF xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#">
  <rdf:Description rdf:about="DJI Meta Data"
    xmlns:tiff="http://ns.adobe.com/tiff/1.0/"
    xmlns:exif="http://ns.adobe.com/exif/1.0/"
    xmlns:xmp="http://ns.adobe.com/xap/1.0/"
    xmlns:xmpMM="http://ns.adobe.com/xap/1.0/mm/"
    xmlns:dc="http://purl.org/dc/elements/1.1/"
    xmlns:crs="http://ns.adobe.com/camera-raw-settings/1.0/"
    xmlns:drone-dji="http://www.dji.com/drone-dji/1.0/"
   xmp:ModifyDate="2018-12-14"
   xmp:CreateDate="2018-12-14"
   tiff:Make="DJI"
   tiff:Model="FC6310R"
   dc:format="image/jpg"
   drone-dji:AbsoluteAltitude="+255.92"
   ...

のように見えます。これらの情報を読み書きしたい、というのが目的です。

python では

便利なツールはもちろん存在します。

python3 -m pip install python-xmp-toolkit

Adobe からXMP-Toolkit-SDKがソースごと公開もされています。組み込みやスマホアプリで利用できるのかな。
https://github.com/adobe/XMP-Toolkit-SDK

xmp のタグの値を読む

使い方はpython-xmp-toolkit のDocに書かれていて、それをフォローするだけで問題なかったです。メモしておきます。

ファイルを読む

read_xmp01.py

from libxmp import XMPFiles
xmpfile = XMPFiles( file_path="./100_0020_0001.JPG", open_forupdate=True )
xmp = xmpfile.get_xmp()

エラーが出なければ読めています。^^;) いろいろxmp の中身を調べることができます。

値の取得

先のサンプルでは

tiff:Make="DJI"

という行がありました。これを "tiff:Make" というkey に対して value である　"DJI"を直接取得したい。そのためには、２ステップ必要です。

namespace の取得

まず、tiff のnamespace (URI)を取得します。

In[]: xmp.get_namespace_for_prefix("tiff")
Out[]: 'http://ns.adobe.com/tiff/1.0/'

あらかじめadobe が定義されているnamespace もありこれはライブラリからも参照できます。

In []: libxmp.consts.XMP_NS_EXIF
Out[]: 'http://ns.adobe.com/exif/1.0/'

変数名は直接ソースの "Common Namespaces" をみて知ることもできます。

property を取得する

tiff:Make というタグの値を読むには、さきほどのtiffに対応するnamespace を指定してあげます。

In []: xmp.get_property(xmp.get_namespace_for_prefix('tiff'),'Make')
Out[]: 'DJI'

namespace内にそのキーワードがあるかを事前に確認することで、エラーを未然に防げます。

In []: xmp.does_property_exist(xmp.get_namespace_for_prefix('tiff'),'TEST')
Out[]: False

書く

正確には「書き換える」ができる「はず」なのですが、実は。。。

まず、xmp の値をsetします。

In []: xmp.set_property(xmp.get_namespace_for_prefix('tiff'),'Make', 'some/makers')
In []: xmp.get_property(xmp.get_namespace_for_prefix('tiff'),'Make')
Out[]: 'some/makers'

ファイルに書くには、xmpfileにput して file を close する必要があるようです。

In []: if xmpfile.can_put_xmp(xmp):
  ...:     xmpfile.put_xmp(xmp)
  ...:     xmpfile.close_file()

で良いはずなのですが。。。結果は正しく書き換えられていませんでした。手順が違うのかバグなのか分かりませんが、とりあえずここまでメモしておきます。実行時にエラーは発生しません。

メモ

とりあえず読むことだけできる、という状態でした。以下自分の課題：

• 各 namespace にあるタグを事前にすべて取得できるか？
• 書き込みの動作確認

以上。(2020/08/08)

ABC173Cは非常に学びが多く、今後の参考になりそうなので個人的な備忘録代わりに。

ABC173C

import copy  # 多次元リストの値渡し/深いコピー(Deep Copy): copy.deepcopy()のため

h, w, k = map(int, input().split())
grd = [list(i for i in input()) for _ in range(h)]
ans = 0
for i in range(2**h):  # 行の選び方が2**h通りで、そのうちi番目
    grd1 = copy.deepcopy(grd)  # grdの値渡し/深いコピー(Deep Copy)
    for i1 in range(h):  # 2**h通りは２進数にしたときh桁で、そのうちi1桁目
        if i >> i1 & 1:  # 「2**h通りのうちi番目」のiを２進数にしたとき、どの桁が1になっているか
            grd1[i1] = ['r'] * w
    for j in range(2**w):  # 列の選び方が2**w通りで、そのうちj番目
        grd2 = copy.deepcopy(grd1)  # grd1の値渡し/深いコピー(Deep Copy)
        for j1 in range(w):  # 2**w通りは２進数にしたときw桁で、そのうちj1桁目
            if j >> j1 & 1:  # 「2**w通りのうちj番目」のjを２進数にしたとき、どの桁が1になっているか
                for n in range(h):
                    grd2[n][j1] = 'r'
        if sum(grd2,[]).count('#') == k:  # sum()で2次元リストを一次元に平坦化
            ans += 1
print(ans)

bit全探索

bit全探索に関してはいろいろな記事があるが、概して初心者かつPythonしか使えない人には難しい。以下はそんな人にもわかりやすく、まずbit全探索って何をするのかを掴むのにはちょうどよい。
Python de アルゴリズム（bit全探索）

listのコピー

リスト.copy()やリスト[:]を使って、オブジェクトの破壊・非破壊など意識して書いたつもりではあったが、思ったとおりに動いてくれなくて（上のコードではgrdが意図せずに書き換わって）困った。
コピーには参照渡し/浅いコピー(Shallow Copy)と値渡し/深いコピー(Deep Copy)があり、更に値渡し/深いコピー(Deep Copy)でもcopy.copy()は一次元リストと多次元リストで挙動が異なる。
Python - リストを複製する
【Python】多次元リストのコピー

多次元リストの一元化

この問題のように多次元のリストを一元的に集計したいようなときに便利。いくつかやり方はあるようだが、sum()の意外な使い方も。
Pythonでflatten（多次元リストを一次元に平坦化）

ABC173Cは非常に学びが多く、今後の参考になりそうなので個人的な備忘録代わりに。

ABC173C

import copy  # 多次元リストの値渡し/深いコピー(Deep Copy): copy.deepcopy()のため

h, w, k = map(int, input().split())
grd = [list(i for i in input()) for _ in range(h)]
ans = 0
for i in range(2**h):  # 行の選び方が2**h通りで、そのうちi番目
    grd1 = copy.deepcopy(grd)  # grdの値渡し/深いコピー(Deep Copy)
    for i1 in range(h):  # 2**h通りは２進数にしたときh桁で、そのうちi1桁目
        if i >> i1 & 1:  # 「2**h通りのうちi番目」のiを２進数にしたとき、どの桁が1になっているか
            grd1[i1] = ['r'] * w
    for j in range(2**w):  # 列の選び方が2**w通りで、そのうちj番目
        grd2 = copy.deepcopy(grd1)  # grd1の値渡し/深いコピー(Deep Copy)
        for j1 in range(w):  # 2**w通りは２進数にしたときw桁で、そのうちj1桁目
            if j >> j1 & 1:  # 「2**w通りのうちj番目」のjを２進数にしたとき、どの桁が1になっているか
                for n in range(h):
                    grd2[n][j1] = 'r'
        if sum(grd2,[]).count('#') == k:  # sum()で2次元リストを一次元に平坦化
            ans += 1
print(ans)

bit全探索

bit全探索に関してはいろいろな記事があるが、概して初心者かつPythonしか使えない人には難しい。以下はそんな人にもわかりやすく、まずbit全探索って何をするのかを掴むのにはちょうどよい。
Python de アルゴリズム（bit全探索）

listのコピー

リスト.copy()やリスト[:]を使って、オブジェクトの破壊・非破壊など意識して書いたつもりではあったが、思ったとおりに動いてくれなくて（上のコードではgrdが意図せずに書き換わって）困った。
コピーには参照渡し/浅いコピー(Shallow Copy)と値渡し/深いコピー(Deep Copy)があり、copyライブラリのcopy.copy()とcopy.deepcopy()では特に多次元リストで挙動が異なる。
Python - リストを複製する
【Python】多次元リストのコピー

多次元リストの一次元化

この問題のように多次元のリストを一元的に集計したいようなときに便利。いくつかやり方はあるようだが、sum()の意外な使い方も。
Pythonでflatten（多次元リストを一次元に平坦化）