【ArcMapとArcpy】オープンデータを使って自分だけの地図を作ろう！ ～阪神・淡路大震災編～

さて、本日は神戸市のオープンデータを利用してみようと思います。阪神・淡路大震災「1.17の記録」に阪神・淡路大震災の記録写真が公開されています。前回行ったようにこの写真にジオタグを追加して、その緯度経度にポイントを配置してみます。

データのダウンロード

上記サイトにおいて、画像データの一括ダウンロード＞CSVでダウンロードを行います。
こんな感じのCSVがダウンロードできました。緯度経度の項目があるにも関わらず、ここは空白になっています・・・仕方ないので以前行ったジオコーディングを使用したいと思います。ただ、場所が空白、もしくは、ジオコーディングできない（衣掛町5丁目2「こうべ地震災害対策広報」みたいな記述）写真は取込対象外とします。

ポイントとなる技術

主に以下のような技術を使用します。
①画像のダウンロード
②画像にジオタグ(geotag)の追加
③ジオコーディング
④ジオタグを追加した画像をポイントに変換

ご注意ください

本エントリー作成時には「②画像にジオタグ(geotag)の追加」が正しく動作したのですが、現状、正しく動作が判明しないことがわかりました。以下エントリーの方法で同じことが実現できますので、「②画像にジオタグ(geotag)の追加」に関してはこちらを参照してください。

www.gis-py.com

日本地図のダウンロード

こちらも以前行ったようにESRIのサイトからダウンロードを行ってください。
全国市区町村界データ | ESRIジャパン

サンプルコード

上述した「ポイントとなる技術」を盛り込んだサンプルコードです。パスなどはご自身の環境にあったものに書き換えてください。

# -*- coding: utf-8 -*-
import pandas as pd
import urllib2
import cookielib
import os
import pyexiv2
from PIL import Image
import geocoder
import arcpy

def main():
    #ダウンロードしたcsvの読込
    df = pd.read_csv(r"D:\python\geotag\20160120142702.csv",encoding="SHIFT-JIS")

    #各種設定
    folder = "D:\python\geotag"
    outFeatures = r"C:\ArcPySample\ArcPySample.gdb\geotag"
    badPhotosList = r"C:\ArcPySample\ArcPySample.gdb\empty"
    photoOption = "ONLY_GEOTAGGED"
    attachmentsOption = ""

    cj = cookielib.CookieJar()
    opener = urllib2.build_opener(urllib2.HTTPCookieProcessor(cj))

    #CSVの行数分ループ
    for key,rowD in df.iterrows():

        #リスエスト
        req = urllib2.Request(rowD.image)

        #出力先を指定
        with open(os.path.join(folder,rowD.filename), 'wb') as f:
            f.write(opener.open(req).read())

        #場所が無いデータは無視
        if type(rowD.place) != float:

            #緯度経度取得
            loc = get_coordinate(rowD.place)

            #井戸経度が取得できた場合
            if loc.lat is not None:
                #ジオタグを画像に追加
                set_gps_location(os.path.join("D:\python\geotag",rowD.filename),loc.lat,loc.lng)
            loc = None

    #ジオタグ付き写真をポイントに変換
    arcpy.GeoTaggedPhotosToPoints_management(folder, outFeatures, badPhotosList, photoOption, attachmentsOption)

def get_coordinate(location_name):
    try:
        #地名から座標を取得する
        ret = geocoder.google(location_name)
    except KeyError, e:
        print "KeyError(%s)" % str(e)
        return

    return ret

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
    min = int(t1)
    sec = round((t1 - min)* 60, 5)
    return (deg, min, sec, loc_value)

def set_gps_location(file_name, lat, lng):

    try:
        lat_deg = to_deg(lat, ["S", "N"])
        lng_deg = to_deg(lng, ["W", "E"])

        # 緯度、経度を10進法→60進法(度分秒)に変換
        exiv_lat = (pyexiv2.Rational(lat_deg[0]*60+lat_deg[1],60),pyexiv2.Rational(lat_deg[2]*100,6000), pyexiv2.Rational(0, 1))
        exiv_lng = (pyexiv2.Rational(lng_deg[0]*60+lng_deg[1],60),pyexiv2.Rational(lng_deg[2]*100,6000), pyexiv2.Rational(0, 1))
        metadata = pyexiv2.ImageMetadata(file_name)
        metadata.read()

        #メタデータを付与
        metadata["Exif.GPSInfo.GPSLatitude"] = exiv_lat
        metadata["Exif.GPSInfo.GPSLatitudeRef"] = lat_deg[3]
        metadata["Exif.GPSInfo.GPSLongitude"] = exiv_lng
        metadata["Exif.GPSInfo.GPSLongitudeRef"] = lng_deg[3]
        metadata["Exif.Image.GPSTag"] = 654
        metadata["Exif.GPSInfo.GPSMapDatum"] = "WGS-84"
        metadata["Exif.GPSInfo.GPSVersionID"] = '2 0 0 0'

        #メタデータを上書き
        metadata.write()
    except:
        print u"エラー。たまに出る。とりあえず飛ばす"

if __name__ == '__main__':
    main()

結果を確認してみます。それっぽい箇所にポイントが配置されましたが、全然関係ない場所にも配置されてしまっています。これは例えば「若松町」という場所が他の都道府県にも存在するという原因かと思います。完璧に行う場合はもう少し工夫が必要ですがとりあえず気にせず進みます（笑）

うまくいったっぽい箇所をGoogleMapと比較して確認してみます。「ポートターミナル周辺」でジオコーディング、GoogleMapで検索した結果です。

大体あってますかね。ちなみにこんな写真です。ジオタグもちゃんと追加されていますね。

ArcMapでポイントの個別属性から画像を参照することができます。

これだけで終わるのもあれなので、ここで一工夫してみようと思います。まずフィールドの追加を行います。

arcpy.AddField_management("geotag", "HTML_Path", "String")

こうなったらOKです。

その後、以下のようにフィールド演算を行います。

codeblock = """
def setHTMLPath(Path):
    return "<img src='"  + Path + "'" + "Width='300'>"
"""

expression = "setHTMLPath(!Path!)"
arcpy.CalculateField_management("geotag", "HTML_Path", expression,"PYTHON_9.3",codeblock)

こうなったらOKです。

ArcMapでHTMLポップアップ機能を使用します。

任意のポイントを選択するとポップアップが起動します。属性にもっと詳細な情報を入れたらより充実したマップになりますね。

今回わかったことは
①ダウンロードしたCSVの情報がかなり大雑把。そのため、ジオコーディングに一工夫が必要。
②ジオタグ追加→ポイント変換のコンボはかなり面白い！GPSの精度がさらに向上すれば、こういった案件も増えてきそうな予感。
③写真を見るとあらためて震災の怖さというものを感じた。

もっとオープンデータが充実すればさらに精度の高いマップを作成することができると思うので、今後そこを期待したいですね。

以上、「【ArcMapとArcpy】オープンデータを使って自分だけの地図を作ろう！　～阪神・淡路大震災編～ 」でした。