GIS奮闘記

元GISエンジニアの技術紹介ブログ。主にPythonを使用。

スポンサーリンク

ArcGIS Runtime SDK for .NET を使ってみよう!

さて、本日は ArcGIS Runtime SDK for .NET を使ってみようと思います。今まで興味はあったもののなかなか触る機会のなかった SDK なのでとても楽しみです。

ArcGIS Runtime SDK for .NET とは

ArcGIS Runtime SDK for .NET は Windows 及び iOS、Android プラットフォーム上で動作するネイティブ GIS アプリケーションの開発キットであり、Windows デスクトップや Windows タブレット、iOS、Android などの多様なプラットフォーム向けに様々な GIS 機能を持ったアプリケーションを開発することができます。

ArcGIS Runtime SDK for .NET | ESRIジャパン

特徴

ArcGIS Runtime SDK for .NET はファイルジオデータベースを読み込んで SA アプリを開発することもできるのですが、どちらかというと Web のデータにアクセスすることを前提にしている SDK だと思います。GIS は Web 化の傾向にあるので、時代のニーズにあった SDK といえるのではないでしょうか。

Web に重きを置いているので、バリバリデータを編集するようなアプリの開発を検討しているという方は ArcGIS Runtime SDK for .NET よりも ArcGIS Pro SDK を使用した方がいいかと思います。興味のある方は以下のエントリーを読んでみてください。

www.gis-py.com

当ブログでは ArcGIS Pro SDK のサンプルなどは公開していないのですが、以下サイトでサンプルコードを公開しているので、ぜひ参考にしてみてください。

1−ArcGIS Pro SDK for .NET を使用した機能開発 ~マップと... | GeoNet, The Esri Community | GIS and Geospatial Professional Community

開発の前に必要なこと

ArcGIS Developer Subscription というサブスクリプションプログラムに登録する必要があります(無料)。製品詳細 | ESRIジャパン
また、SDK のインストールも行います。https://developers.arcgis.com/downloads/apis-and-sdks

今回やってみること

マップの作成とレイヤーの読み込み(始めてなのでまずは簡単なところから挑戦してみようと思います)

読み込むレイヤーは以下の日本地図です(ArcGIS Online で表示)。 f:id:sanvarie:20200119204052p:plain

実行環境

Windows10 64bit
Visual Studio 2017
ArcGIS Runtime SDK for .NET 100.7.0

サンプル

以下のサイトを参考にしてみました。 Add layers to a map | ArcGIS for Developers

using System;
using System.Collections.Generic;
using System.ComponentModel;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using System.Runtime.CompilerServices;
using Esri.ArcGISRuntime.Data;
using Esri.ArcGISRuntime.Geometry;
using Esri.ArcGISRuntime.Location;
using Esri.ArcGISRuntime.Mapping;
using Esri.ArcGISRuntime.Security;
using Esri.ArcGISRuntime.Symbology;
using Esri.ArcGISRuntime.Tasks;
using Esri.ArcGISRuntime.UI;

namespace RuntimeTest
{
    /// <summary>
    /// Provides map data to an application
    /// </summary>
    public class MapViewModel : INotifyPropertyChanged
    {
        public MapViewModel()
        {
            CreateNewMap();
        }

        private async void CreateNewMap()
        {
            // オープンストリートマップをベースマップにする
            Map newMap = new Map(Basemap.CreateOpenStreetMap());

            // フィーチャレイヤーを読込
            FeatureLayer featurelayer = new FeatureLayer(new Uri("https://services.arcgis.com/P3ePLMYs2RVChkJx/ArcGIS/rest/services/JPN_Boundaries_ECM/FeatureServer/0"));
            await featurelayer.LoadAsync();

            // フィーチャレイヤーをマップに追加
            newMap.OperationalLayers.Add(featurelayer);

            // 起動時の表示範囲を設定
            newMap.InitialViewpoint = new Viewpoint(featurelayer.FullExtent);

            Map = newMap;
        }

        private Map _map = new Map(Basemap.CreateStreetsVector());

        /// <summary>
        /// Gets or sets the map
        /// </summary>
        public Map Map
        {
            get { return _map; }
            set { _map = value; OnPropertyChanged(); }
        }

        /// <summary>
        /// Raises the <see cref="MapViewModel.PropertyChanged" /> event
        /// </summary>
        /// <param name="propertyName">The name of the property that has changed</param>
        protected void OnPropertyChanged([CallerMemberName] string propertyName = null)
        {
            var propertyChangedHandler = PropertyChanged;
            if (propertyChangedHandler != null)
                propertyChangedHandler(this, new PropertyChangedEventArgs(propertyName));
        }

        public event PropertyChangedEventHandler PropertyChanged;
    }
}

アプリが起動し、想定通りレイヤーも読み込まれました。ほとんど時間をかけずにここまでできました。本当に素晴らしいですね。自分でやったこととしては CreateNewMap() の作成とそれをコンストラクタで呼び出すことくらいです。あとはデフォルトのままで問題ありません。

f:id:sanvarie:20200119203905p:plain

今回初めて ArcGIS Runtime SDK for .NET を使ってみたのですが、単純にレイヤーを読み込んで少しカスタマイズするくらいだったら一日二日あればそれなりのものが作れてしまうのでは思うくらい簡単に扱える SDK だと思いました。さすが ArcGIS といったところでしょうか。次回以降はもっと機能を作りこもうと思います。本日は以上です。

ezdxf を使って Python で DXF を扱ってみよう!

さて、本日は ezdxf を使って Python で DXF を扱ってみようと思います。

DXF とは

以下エントリーでも少し紹介しているのですが、DXFは、Autodesk社が開発したファイル形式の1つで、CAD 系のデータですね。しかし、GIS の世界でも DXF はよく使われます。

www.gis-py.com

ezdxf とは

DXF を扱うための Python ライブラリです。データの読み書きを行うことができます。

pypi.org

サンプルデータ

ArcGIS Pro で以下のようなデータを用意しました。

1.浄水場名(アノテーション)
2.浄水場(ポイント)
3.河川(ライン)
4.神奈川県(ポリゴン)

f:id:sanvarie:20200118081144p:plain

f:id:sanvarie:20200118081240p:plain

これをジオプロで「ExportCAD.dxf」という名称のDXFに出力します。

f:id:sanvarie:20200118081650p:plain

できあがった DXF を ArcGIS Pro で読み込むとこのようになります。 f:id:sanvarie:20200118082247p:plain

ちなみにですが、「2.浄水場(ポイント) 」「3.河川(ライン) 」は国土数値情報からダウンロードしてきました。

国土数値情報 上水道関連施設データの詳細

国土数値情報 河川データの詳細

ダウンロードの自動化について過去に紹介しているので、興味のある方はぜひ読んでみてください。 www.gis-py.com

また、「1.浄水場名(アノテーション) 」ですが、「2.浄水場(ポイント) 」から作成しました。作成方法は以下のエントリーで紹介しているので、ぜひ読んでみてください。

www.gis-py.com

実行環境

Windows10 64bit
Python3.6.6
ezdxf 0.10.2
Jupyter 6.0.1

サンプルコード

ezdxf の使い方を紹介します。

レイヤー名称一覧を抽出

import ezdxf

doc = ezdxf.readfile(r"D:\python\data\CAD\ExportCAD.dxf")
for layer in doc.layers: 
    print(layer.dxf.name)

このようになります。「0」と「Defpoints」というレイヤーが抽出されていますが、これは フィーチャクラスを DXF に変換する際にできたもので、特に気にする必要はありません。

f:id:sanvarie:20200118084615p:plain

特定のレイヤのデータを抽出する

「河川」レイヤーのデータを抽出します。

import ezdxf

doc = ezdxf.readfile(r"D:\python\data\CAD\ExportCAD.dxf")
msp = doc.modelspace()

rivers = msp.query('*[layer=="河川"]')
for river in rivers:
    print(river)

こんな感じで結果の結果になりました。

f:id:sanvarie:20200118100421p:plain

色々調べてみたのですが、一つ一つのアイテムの属性にアクセスすることはできませんでした。何かやり方があるのかもしれませんが、今回はあきらめます。

DXF を作成する

ラインを一本作図してそれを DXF として出力します。

import ezdxf

doc = ezdxf.new()
msp = doc.modelspace()
doc.layers.new(name='MyLines', dxfattribs={'linetype': 'DASHED', 'color': 7})
msp.add_line((0, 0), (10, 0), dxfattribs={'layer': 'MyLines'})
doc.saveas(r"D:\python\data\CAD\ExportCAD2.dxf")

想定通りラインが作図された状態で DXF を作成することができました。 f:id:sanvarie:20200118101319p:plain

まとめ

使ってみた感想としては、簡易的に DXF を扱う分には十分だけど、細かな操作を行いたい場合はちょっと厳しいかもしれないと思いました。DXF は専用のソフトで を扱うのが一番かと思います。ただ、今回初めて Python で DXF を扱ってみたのでそれ自体はとてもいい経験になったと思います。興味のある方はぜひ ezdxf を使ってみてください。本日は以上です。

ArcGIS Pro でレイヤーにフィルターをかける方法

さて、本日は ArcGIS Pro の機能紹介になります。レイヤーを参照している際、特定のフィーチャのみを表示したいということはないでしょうか?例えば、日本地図を表示している際に東京都のフィーチャのみを表示したいといった感じです。ご存知の方も多いかと思うのですが、意外とご存じない方もいらっしゃるので、ここで紹介したいと思います。

対象データ

「日本地図」というレイヤーを使用します。こちらは以下のエントリーでも紹介しているESRIジャパンさんの「全国市区町村界データ」をフィーチャクラスに変換したものです。

www.gis-py.com

f:id:sanvarie:20200115205554p:plain

実行環境

Windows10 64bit
ArcGIS Pro 2.3.3

操作方法

1.「日本地図」レイヤーを右クリックして「プロパティ」を押下
f:id:sanvarie:20200115205000p:plain

2.フィルター設定 > 「新しいフィルター設定」ボタンを押下
f:id:sanvarie:20200115205620p:plain

3.「KEN が 東京都 と等しい」 と設定し、OK。
f:id:sanvarie:20200115205649p:plain

東京都のポリゴンだけが表示されました。
f:id:sanvarie:20200115205358p:plain

まとめ

レイヤーのフィルター設定はこのように簡単に行うことができます。状況に応じて、必要なフィーチャの表示/非表示を切り替えることによって、より効率的に作業を進めることができるようになるかと思います。

本日は以上です。

Python で SQLite を使ってみよう!

明けましておめでとうございます。少し更新が止まっていましたが、また再開しようと思います。

去年は私にとっては色々変化の大きい年でした。自分の勉強不足なども痛感し、今年はより一層努力しなければと思います。

さて、本日は Python で SQLite を使ってみようと思います。GIS とはあまり関係ないと思われるかもしれませんが、GIS データを扱うにあたってデータベースは避けては通れない道かと思います。今回だけではなく、今後もデータベースについて色々紹介できればと思います。

SQLite とは

SQLiteは、パブリックドメインの軽量な関係データベース管理システム のことですね(SQLite - Wikipedia)。簡単に言うと、SQLSERVER や Oracle などと異なり、ファイルベースで簡単に扱えるデータベース、といった感じでしょうか。

使用するライブラリ

Python 標準ライブラリである sqlite3 を使用します。

SQLite のGUI 管理ツール

SQLite のGUI 管理ツールですが、DB Browser for SQLite を使用します。

f:id:sanvarie:20200113100352p:plain

関連エントリー

以下のエントリーで Python で SQLSERVER の操作について紹介していますので、興味がありましたらぜひ読んでみてください。

www.gis-py.com

環境

Windows10 64bit
Python 3.6.6
DB Browser for SQLite 3.11.2

sqlite3 の使用方法

sqlite3 を使用して以下のような操作をするための方法を書いてみました。

  1. データベースの作成/接続
  2. テーブルの作成
  3. 挿入
  4. 更新
  5. 削除
  6. 検索

データベースの作成/接続

対象のSQLiteのファイル(sample_data.db)がなかったら作成して、あったら接続します。

# -*- coding: utf-8 -*-
import sqlite3

# sample_data.db がなかったら作成、あったら接続
db_path = "D:\data\SQLite\sample_data.db"
conn = sqlite3.connect(db_path)
conn.close()

実行前

f:id:sanvarie:20200113102206p:plain

実行後

f:id:sanvarie:20200113102229p:plain

sample_data.db が作成されたことがわかります。sample_data.db を DB Browser for SQLite で読み込んだ結果が以下になります。

f:id:sanvarie:20200113102352p:plain

テーブルの作成

「SAMPLE」というテーブルを作成するサンプルです。コネクションを取得して、作成した SQL を実行するだけですね。

# -*- coding: utf-8 -*-
import sqlite3

db_path = "D:\data\SQLite\sample_data.db"
conn = sqlite3.connect(db_path)

def create_table():
    sql = '''CREATE TABLE IF NOT EXISTS SAMPLE
                 (ID    INTEGER PRIMARY KEY,
                  TEST1 TEXT,
                  TEST2 TEXT,
                  TEST3 TEXT)'''

    conn.execute(sql)
    conn.commit()
    conn.close()

create_table()

以下のような結果になりました。カラムもしっかり作成されていることがわかります。 f:id:sanvarie:20200113103215p:plain

挿入

テーブルの作成と基本的には一緒で SQL をINSERT 文にするだけですね。

# -*- coding: utf-8 -*-
import sqlite3

db_path = "D:\data\SQLite\sample_data.db"
conn = sqlite3.connect(db_path)

def insert_record():
    sql = "INSERT INTO SAMPLE VALUES (1, 'あ', 'い', 'う')"
    conn.execute(sql)

    sql = "INSERT INTO SAMPLE VALUES (2, 'てすと1', 'てすと2', 'てすと3')"
    conn.execute(sql)

    sql = "INSERT INTO SAMPLE VALUES (3, '12345', '22222', '55555')"
    conn.execute(sql)

    conn.commit()
    conn.close()

insert_record()

想定通りの結果になりました。 f:id:sanvarie:20200113125244p:plain

更新

こちらも SQL を UPDATE 文に変えただけですね。

# -*- coding: utf-8 -*-
import sqlite3

db_path = "D:\data\SQLite\sample_data.db"
conn = sqlite3.connect(db_path)

def update_record():
    sql = "UPDATE SAMPLE SET TEST1 = 'か' WHERE ID = 1"
    conn.execute(sql)
    conn.commit()
    conn.close()

update_record()

しっかり UPDATE されていることがわかります。 f:id:sanvarie:20200113125551p:plain

削除

こちらも SQL を DELETE 文に変えただけですね。

# -*- coding: utf-8 -*-
import sqlite3

db_path = "D:\data\SQLite\sample_data.db"
conn = sqlite3.connect(db_path)

def delete_record():
    sql = "DELETE FROM SAMPLE WHERE ID = 3"
    conn.execute(sql)
    conn.commit()
    conn.close()

delete_record()

ID = 3 のレコードが削除されていることがわかります。 f:id:sanvarie:20200113125731p:plain

検索

最後に検索ですが、こちらは他のとは少し異なりますのでご注意ください。 SQL を SELECT 文に変えて、そのあと execute した結果をループします。

# -*- coding: utf-8 -*-
import sqlite3

db_path = "D:\data\SQLite\sample_data.db"
conn = sqlite3.connect(db_path)

def select_record():
    sql = "SELECT * FROM SAMPLE WHERE ID = 2"
    for row in conn.execute(sql):
        print(row)

select_record()

このように SELECT した結果を取得することができます。 f:id:sanvarie:20200113130200p:plain

Python での SQLite の操作について一通り紹介しましたが、すごく簡単にできることがわかりました。私は業務で SQLite はあまり使ったことがないのですが、おそらく世間的にはたくさん使われているデータベースだと思います。やはりファイルベースで扱えるというメリットが人気の秘密でしょうか。

本日は以上です。

GeoPandas、Fiona、Shapely を 使ってファイルジオデータベース内のデータを Shape 出力する方法

さて、本日は GeoPandas、Fiona、Shapely を 使ってファイルジオデータベース内のデータを Shape 出力してみようと思います。どれも GIS 界では有名でとても優秀なライブラリですね。これらについては過去のエントリーで紹介しているので、詳しくはそちらを読んでみてください。

www.gis-py.com

www.gis-py.com

www.gis-py.com

実行環境

実はGeoPandas、Fiona、Shapely 以外にも Pandas を使います。ただ、Pandas は Shape 出力には関係ないところで使用します。詳細は後述します。

Windows 10
Python 3.6.6
ArcGIS Pro 2.4
GeoPandas 0.6.1
Fiona 1.8.4
Shapely 1.6.4
Pandas 0.25.2

出力するデータ

このようなデータをフィーチャクラス単位で出力してみようと思います。フィーチャデータセットも含めてみました。また、アノテーションが入っていますが、アノテーションは出力対象外になるように制御をかけようと思います。

f:id:sanvarie:20191030174559p:plain

一部ではありますが、データはこのように日本地図のポリゴン、各県の外郭線のライン、横須賀市内に作成したポイントなどがあります。

f:id:sanvarie:20191030175020p:plain

f:id:sanvarie:20191030175047p:plain

属性には日本語データが入っています。 f:id:sanvarie:20191030175507p:plain

f:id:sanvarie:20191030175532p:plain

今回直面した問題

最初はこんなにたくさんのライブラリを使うつもりはなく、単純に GDAL とかで出力しようかと思ったのですが、日本語データが含まれているとうまく出力できず、また、GeoPandas や Fiona でも同じように試したのですが、日本語データがネックになり断念しました。もしかしたら何かやり方があるのかもしれないのですが、とりあえず今回は別のアプローチをとりました。ちなみに日本語データが含まれていなければ以下のような感じにすれば大丈夫かと思います(GeoPandas と Fiona を使用した例)。

import fiona
import geopandas as gpd
import os

gdb = r"D:\python\data\Sample.gdb"
layers = fiona.listlayers(gdb)
anno_flg = 0

for layer in layers:

    gdf = gpd.read_file(gdb, driver="FileGDB", layer=layer)

    for column in gdf.columns:

        if column == "AnnotationClassID":# アノテーションは対象外
            anno_flg = 1

        if column == "Shape_Length":
            gdf.drop('Shape_Length', axis=1, inplace=True)

        if column == "Shape_Area":
            gdf.drop('Shape_Area', axis=1, inplace=True)

    if anno_flg == 1:
        anno_flg = 0
        continue

    out = os.path.join(r"D:\python\data\shape", layer + ".shp")
    gdf.to_file(out)

直面した問題を回避するために

出力するファイルを以下の二つに分けました。
1.Shape・・・図形データを出力
2.CSV・・・属性データを出力

このようにファイルを分けて出力すれば日本語データが含まれていても問題ないかと思います。双方にキーとなるカラムも含めていますので、データ出力後はそのキーを使って図形と属性を紐づけるだけですね。Pandas は CSV 出力の処理で使用します。

サンプルコード

import os
import geopandas as gpd
import fiona
from shapely.geometry import Point, LineString, Polygon, MultiPolygon, MultiLineString, MultiPoint
import pandas as pd

def make_shape_file():
    gdb = r"D:\python\data\ToShapeFile.gdb" # 出力対象のGDB
    layers = fiona.listlayers(gdb)

    for layer in layers:
        new_data = gpd.GeoDataFrame()
        new_data['geometry'] = None

        gdf = gpd.read_file(gdb, driver="FileGDB", layer=layer)

        if 'AnnotationClassID' in gdf[0:1].columns: # アノテーションは処理対象外
            continue

        geometry_type = gdf[0:1]['geometry'].geom_type[0]

        for index, row in gdf.iterrows():

            if geometry_type == "MultiPolygon":
                new_data.loc[index, 'geometry'] = MultiPolygon(row.geometry)

            elif geometry_type == "Polygon": # 今回は不使用
                new_data.loc[index, 'geometry'] = Polygon(row.geometry)

            elif geometry_type == "LineString": # 今回は不使用
                new_data.loc[index, 'geometry'] = LineString(row.geometry)

            elif geometry_type == "MultiLineString":
                new_data.loc[index, 'geometry'] = MultiLineString(row.geometry)

            elif geometry_type == "Point":
                new_data.loc[index, 'geometry'] = Point(row.geometry)

            elif geometry_type == "MultiPoint": # 今回は不使用
                new_data.loc[index, 'geometry'] = MultiPoint(row.geometry)

        # Shapeファイル出力
        out = os.path.join(r"D:\python\data\shape", layer + ".shp")
        new_data.to_file(out)

        # CSVファイル出力
        df = pd.DataFrame(gdf)
        df.to_csv(os.path.join(r"D:\python\data\csv", layer + ".csv"))

make_shape_file()

結果

Shape も CSV もきちんと出力されましたね。

f:id:sanvarie:20191030181847p:plain

f:id:sanvarie:20191030181816p:plain

Shape を ArcGIS Pro で読み込むと図形データのみ出力されていることがわかります。

f:id:sanvarie:20191030182008p:plain

CSV はこのようになっています。左端の列がキーです。

f:id:sanvarie:20191030182422p:plain

色々苦労はしたものの、プログラム自体はけっこうあっさりしてますね。本当に便利なライブラリ達です。次は GDAL で同じことを実現させてみたいですね。やり方がごっそりと変わると思いますが面白そうです。

今回は出力するファイルを Shape と CSV に分ける方法をとったのですが、もし Shape ファイルを出力するだけで問題ないやり方をご存知の方がいましたらコメントいただけるととても助かります。よろしくお願いします。本日は以上です。

piexif を使って画像に付与されているジオタグからポイントを作成して Shape ファイルとして出力する方法

さて、本日は前回のエントリーに続きジオタグのお話です。前回は画像にジオタグを追加する方法を紹介しましたが(以下エントリー参照)、今回は画像に付与されているジオタグを抽出する方法を紹介しようと思います。また、その緯度経度を使用してポイントを作成しShape ファイルとして出力してみようと思います。

www.gis-py.com

使用するライブラリ

piexif
GDAL

piexif で画像に付与されているジオタグを抽出して、その情報を GDAL を使ってポイントを作成し、Shapeファイルとして抽出してみようと思います。なお、ポイント作成とShape 出力の部分は以下エントリーをベースに書いてみようと思います。

www.gis-py.com

使用する画像

こんなサイトがあったので、使ってみました。

www.geoimgr.com

世界各地の風景写真を掲載しているサイトですね。これは京都の桂川の写真らしいです。

f:id:sanvarie:20191026144148p:plain

緯度経度も設定されていることがわかります。

f:id:sanvarie:20191027085913p:plain

こんな感じにデータを格納してみました。これらのジオタグを抽出して、ポイントとして Shape ファイルを作成してみようと思います。

f:id:sanvarie:20191026142150p:plain

実行環境

Windows 10
Python 3.6.6
piexif 1.1.3
GDAL 2.1.1
ArcGIS Pro 2.4

サンプルコード

import os
from glob import glob
import piexif
import osgeo.ogr as ogr

folder = r"D:\data\picture"
shape_file = r"D:\data\shape\geotag.shp"

def make_shape_for_point(info):
    """ポイントShapeファイルを作成"""

    # shapefileドライバ
    driver = ogr.GetDriverByName("ESRI Shapefile")

    # Shape出力先を設定
    data_source = driver.CreateDataSource(shape_file)

    # レイヤ作成
    layer = data_source.CreateLayer("geotag",  geom_type= ogr.wkbPoint)

    # フィールド追加
    field_name = ogr.FieldDefn("file_name", ogr.OFTString)
    field_name.SetWidth(24)
    layer.CreateField(field_name)

    layer.CreateField(ogr.FieldDefn("longitude", ogr.OFTReal))
    layer.CreateField(ogr.FieldDefn("latitude", ogr.OFTReal))

    # Shape作成
    for l in info:
      # フィーチャ作成
      feature = ogr.Feature(layer.GetLayerDefn())

      # 属性を設定
      feature.SetField("file_name", l[2])
      feature.SetField("longitude", l[1])
      feature.SetField("latitude", l[0])

      # 作成するポイントに緯度経度をセット
      wkt = "POINT(%f %f)" %  (float(l[1]) , float(l[0]))

      # ポイント作成
      point = ogr.CreateGeometryFromWkt(wkt)
      feature.SetGeometry(point)
      layer.CreateFeature(feature)
      feature = None

    data_source = None

def conv_deg(v):
    """緯度、経度を60進法(度分秒)→10進法に変換"""

    d = float(v[0][0]) / float(v[0][1])
    m = float(v[1][0]) / float(v[1][1])
    s = float(v[2][0]) / float(v[2][1])
    return d + (m / 60.0) + (s / 3600.0)

def read_geotag():
    """画像を読み込んで緯度経度を取得"""

    info = []

    # 指定するフォルダ内のjpgを取得
    for jpg in glob(folder + "\*.jpg"):

        # 画像を読み込む
        exif_dict = piexif.load(jpg)

        if (len(exif_dict["GPS"]) > 0):

            # 緯度、経度を60進法(度分秒)→10進法に変換
            lat = conv_deg(exif_dict["GPS"][2])
            lon = conv_deg(exif_dict["GPS"][4])

            # ファイル名を取得
            file_name = os.path.basename(jpg)

            info.append([lat, lon, file_name])

    return info

make_shape_for_point(read_geotag())

結果を確認するとこのように Shape ファイルが出力されたことがわかります。

f:id:sanvarie:20191026143502p:plain

出力された Shape ファイルを ArcGIS Pro で読み込んでみると、それっぽい場所にポイントが作成されていることがわかります。

f:id:sanvarie:20191026143133p:plain

京都の桂川のデータを確認してみると、場所は問題なさそうです。また、属性もしっかり格納されていますね。

f:id:sanvarie:20191026144259p:plain

このように piexif を使用すれば前回のようにジオタグを画像に追加したり、今回のように画像に付与されているジオタグを抽出することが簡単にできます。また、 GDAL を使用すれば、その抽出した情報からポイントを作成し、Shape ファイルとして出力することも簡単にできます。

二回連続でジオタグに関して書いてみましたが、いかがでしたでしょうか?もし他に何か知りたいことなどがありましたらコメントしていただければと思います。本日は以上です。

piexif を使って画像にジオタグを追加する方法

さて、本日は 久しぶりに Python を使ってみようと思います。piexif というライブラリがあるのですが、これを使って画像にジオタグを追加する方法を紹介します。ライブラリの詳細は以下のリンクを参照してください。

piexif.readthedocs.io

今回のエントリーに至った経緯

実は以下のエントリーで過去にも同じことをしていたのですが、最近そのコードが動かないことが判明しました(当時使用していたライブラリ(pyexiv2)の仕様が変わった?!もしくは何か依存ライブラリが足りてないのか調べても原因はわかりませんでした・・・)。今回はせっかくなので、別のライブラリ(piexif) を使ってみようと思います。

www.gis-py.com

www.gis-py.com

ジオタグとは

ジオタグとは、写真や動画、あるいは SNS の投稿といったさまざまなメディアに追加することができる、位置情報 (緯度・経度) を示すメタデータのことです。

使用するライブラリ

piexif
Pillow

実行環境

Windows 10
ArcGIS Pro 2.4
Python 3.6.6
piexif 1.1.3
Pillow 6.1.0

使用するデータ

以下の画像を使用します。

f:id:sanvarie:20191026105621p:plain

使用する緯度経度

緯度: 35.678601
経度: 139.740286

ここがどこかといいますと永田町駅ですね。

f:id:sanvarie:20191026091618p:plain

サンプルコード

1ファイルだけを対象としたサンプルコードです。

from fractions import Fraction
import piexif
from PIL import Image

in_file = r"D:\data\picture\arcgis.jpg"         # 使用する画像
out_file = r"D:\data\picture\arcgis_geotag.jpg" # 出力する画像
lat = 35.678601
lon = 139.740286

def to_deg(value, loc):
    """緯度、経度を10進法→60進法(度分秒)に変換"""
    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 change_to_rational(number):
    """有理数に変換"""
    f = Fraction(str(number))
    return (f.numerator, f.denominator)

def atattch_geotag(file_name, lat, lng):
    """ジオタグに追加する情報を作成"""
    lat_deg = to_deg(lat, ["S", "N"])
    lng_deg = to_deg(lng, ["W", "E"])

    # 緯度、経度を10進法→60進法(度分秒)に変換
    exiv_lat = (change_to_rational(lat_deg[0]), change_to_rational(lat_deg[1]), change_to_rational(lat_deg[2]))
    exiv_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.GPSLatitudeRef: lat_deg[3],
        piexif.GPSIFD.GPSLatitude: exiv_lat,
        piexif.GPSIFD.GPSLongitudeRef: lng_deg[3],
        piexif.GPSIFD.GPSLongitude: exiv_lng,
    }

    exif_dict = {"GPS": gps_ifd}
    exif_bytes = piexif.dump(exif_dict)

    # ジオタグ付きの画像を作成
    create_picture_geotag(exif_bytes)

def create_picture_geotag(exif_bytes):
    """ジオタグを貼り付けた画像を作成"""
    im = Image.open(in_file)
    im.save(out_file, exif = exif_bytes)

atattch_geotag(in_file, lat, lon)

結果を確認してみると、出力された画像にジオタグが追加されていることがわかります。

f:id:sanvarie:20191026103815p:plain

追加されたジオタグが正しいか確認

ArcGIS Pro を使って追加されたジオタグが正しいか確認してみます。「ジオタグ付き写真 → ポイント」というジオプロがあるので、ArcPy を使って実行してみようと思います。

import arcpy
arcpy.GeoTaggedPhotosToPoints_management(r"D:\data\picture", r"D:\Sample\Sample.gdb\python_GeoTaggedPhotosToPoints", "", "ONLY_GEOTAGGED", "ADD_ATTACHMENTS")

設定した緯度経度(永田町)にポイントが作成されました。

f:id:sanvarie:20191026104804p:plain

ポップアップを構成して属性ウインドウに写真を表示してみました。
f:id:sanvarie:20191026105055p:plain

この piexif というライブラリですが、とても簡単に画像にジオタグを追加できる優秀なライブラリということがわかりました。以下のエントリーで画像に付与されているジオタグを抽出する方法を紹介しているので、興味のある方はぜひ読んでみてください。

www.gis-py.com

本日は以上です。