GIS奮闘記

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

スポンサーリンク

Pythonでのモジュールとパッケージの作成

本日は「Pythonでのモジュールとパッケージの作成」です。簡単なスクリプトなら一つのファイルで十分なのですが、ある程度の規模になるとクラスを分けたくなりますよね。そこでパッケージとモジュールです。今回はパッケージとモジュールの作り方について書いてみます。

フォルダ構成はこんな感じにしました。
f:id:sanvarie:20160104105624p:plain
f:id:sanvarie:20160104105703p:plain
f:id:sanvarie:20160104105706p:plain

■main.py

# -*- coding: utf-8 -*-
from Lib.ReturnText import ReturnText

rt = ReturnText(u"あけまして")
print rt.content()

■ReturnText.py

# -*- coding: utf-8 -*-
class ReturnText:
    def __init__(self, text):
        self.text = text

    def content(self):
        return self.text + u"おめでとうございます。"


■__init__.py

#!/usr/bin/env python
# -*- coding: utf-8 -*-

「main.py」で「ReturnText.py」(モジュールファイル)をimportしています。「__init__.py」は空で構いませんが、必ず配置してください。

Pythonではフォルダに「__init__.py」があれば、そのフォルダを「パッケージ」として扱えるようになっています。

実行結果を見てみます。

f:id:sanvarie:20160104110510p:plain

正しい結果が得られたことが確認できました。簡単ではありますが、Pythonでのモジュールとパッケージの作成についてでした。

Pythonで緯度経度を平面直角座標(XY座標)に変換する

明けましておめでとうございます。新年最初の投稿は「Pythonで緯度経度を平面直角座標(XY座標)に変換する」です。

Pythonでジオコーディングをやってみる - GIS奮闘記で取得した緯度経度を平面直角座標に変換したいと思います。

pyprojのインストール

pyprojとは、「PythonでProjection(投影座標系、測地系)を扱うためのlibrary」です。pip install pyproj でインストール完了です。

サンプルコードは以下です。国立競技場の緯度経度を平面直角座標へ変換しました。

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

lat = 35.679933
lon = 139.714465

EPSG4612 = pyproj.Proj("+init=EPSG:4612")
EPSG2451 = pyproj.Proj("+init=EPSG:2451")
y,x = pyproj.transform(EPSG4612, EPSG2451, lon,lat)

print "x = %s" % x
print "y = %s" % y

EPSGコードに関しては以下サイトをご参照願います。
GISのための測地成果、測地系、楕円体、投影座標系、EPSGコードのまとめ - 自然環境保全のための周辺技術

結果を見てみます。
f:id:sanvarie:20160104084556p:plain

x = -35503.3049192
y = -10759.6736005

という座標を取得することができました。これを以下サイトで確認してみます。
平面直角座標 ⇔ 緯度・経度変換 with Google Map

f:id:sanvarie:20160104085321p:plain

変換結果が問題ないことがわかりました。簡単すぎてびっくりです。なんて日だ!

以上、本日は「Pythonで緯度経度を平面直角座標(XY座標)に変換する」でした。

投影法、座標系、測地系、それぞれの違いについて

今回はGISの基本的な概念である「投影法、座標系、測地法」とそれぞれの違いについて書いてみようと思います。GISに関わっているとこれらの用語には必ず出くわすと思います。

 

しかし、それぞれの違いについてよくわからない、という方も多いのではないでしょうか?例えば、業務では平面直角座標系(JGD2000)を使用しているが、平面直角座標系、JGD2000の違いを即座に答えられる方は少ないような気がします。

 

以下のようなサイトで確認しても難解な言葉の連続でよくわからないですね。←国土地理院さん。もう少し簡単にお願いします・・・

日本の測地系|国土地理院

 

GISを扱う上で必須なこれらの概念、つまり、「投影法、座標系、測地系」のそれぞれの関係性を表にしてみました。

 

この表だけあれば問題ないという方が大半かと思います。これ以降の説明はおまけくらいに思ってください。

f:id:sanvarie:20151231093709p:plain

 

測地系

測地系(そくちけい)は、地球上の位置を経緯度(経度・緯度)及び標高を用いる座標によって表すとき、前提とする条件のことです。座標系やら投影法やらありますが、まず測地系ありきということですね。

世界測地系

日本測地系より新しい基準。日本測地系は使用せずにこちらに移行しましょう、という傾向がみられる。「WGS1984(WGS84系)」はアメリカの基準でGPS関係などに使用されているとか。なので、とりあえず「JGD2000」を使用しておけば問題ありません。

※最近は東日本大震災による影響でJGD2011なるものが出てきました。すみませんが、本ブログでは割愛させていただきます。

日本測地系

世界測地系より古い基準。ただ、田舎ではまだまだ使用しているところが多い気がする。「日本測地系」=「TOKYO」で問題ないと思います。

 

投影法

投影法とは3次元の立体を2次元の平面上に表現する方法のことをいいます。要するに投影法とは「地理座標系」「UTM座標系」「平面直角座標系」の総称、また、投影座標系は「UTM座標系」「平面直角座標系」の総称ということで問題ないと思います。おそらく「平面直角座標系」を使用しているという方が多いのではないでしょうか。以下にそれぞれの説明を記載しますが、こういったものがあるのか、くらいの認識で大丈夫だと思います。

 

地理座標系

経度は赤道を基点(0°)として南北の、緯度は旧グリニッジ天文台跡を基点として東西の、それぞれ角度を表した数値です。

 

UTM座標系

東西方向に水平に保たれた円筒の中に地球をはめ込んで、地球の中心と地表面上の点とを結んだ直線が円筒と交わる点に印を付けて、円筒を展開したもの。投影誤差を小さくするため地球を60の帯(ゾーン)に分けて、それぞれの平面に投影しています。

 

平面直角座標系

日本国内を土地の状況に合わせて19に分割します。1つの区画は、原点から東西に130km、計260kmの幅になります。それらの分割区画は、第1系~ 第19系と番号付けされ、地図を扱う上では目的とする地点がどの系に含まれるかを知らなければなりません。1/2500地図のような比較的狭い範囲の地図 に適した座標系です。

 

上記の表を見ればなんとなく、投影法、座標系、測地系、それぞれの関係性が見えてきたのではないでしょうか。すべてを一気に覚える必要はなく、まずは世界測地系と日本測地系の違いから覚えてみるのがいいかもしれません。

 

 タイトルとは直接関係ないのですが、投影法、座標系、測地系以外にもGISに関して知りたいことがたくさんあるという方は、以下のエントリーをぜひ読んでみてください。

www.gis-py.com

 

以上、「投影法、座標系、測地系、それぞれの違いについて」でした!

 

GDALでGeoTiffを作ってみよう!

さて、本日は「GDALでGeoTiffを作ってみよう!」です。 GDALとはラスターデータのメタデータの検索や、データフォーマットの変換などが可能になる便利なソフトです。

何を言っているかわからないって?!とりあえずGDALっていうソフトを使うとGeoTiffが作れるくらいの認識で大丈夫です(他にも色々できますが、今回はGeoTiffについてだけをとりあげます)。

 今度はGeoTiffって何?!となる方も多いかと思います。簡単に言うと座標つきの画像ですね。航空写真なんかがそれにあたります。航空写真については以下を参照してください。

 

www.gis-py.com

 

 

 このような地図に航空写真をはめ込みたいと思った時に航空写真自体に座標をもっていないと、地図にうまくはめ込めないですよね。

f:id:sanvarie:20151228111544p:plain

 座標を持っているからこそ、このように地図と画像がばちっとはまります。

 f:id:sanvarie:20151228111746p:plain

 今回はこのような座標つきの画像(GeoTiff)を作ってみたいと思います。

以下サイトから「gdalwin32exe160.zip」をダウンロードして任意のフォルダに格納してください。

 http://download.osgeo.org/gdal/win32/1.6/

 コマンドプロンプトで以下のようなコマンドを叩いてください。

gdal_translate -of "GTiff"  -a_ullr 左上X座標 左上Y座標 右下X座標 右下Y座標 input.tif output.tif

 

こんな感じで作成完了です。

f:id:sanvarie:20151228112442p:plain

ちょっと画像がわかりにくいのですが、原点を確認するとちゃんと作成されたっぽいですね!(今回はSIS MapModellerで確認しています。フリーソフトでも確認できます。)

f:id:sanvarie:20151228112853p:plain

SIS MapModellerって何?という方は以下のエントリーを読んでみてください。

www.gis-py.com

簡単ではありますが、「GDALでGeoTiffを作ってみよう!」でした。

 

【Pythonで画像処理】Pillowを使用してシンボル一覧を作成する

本日は「【Pythonで画像処理】Pillowを使用してシンボル一覧を作成する」です。Pillowとは、Pythonの画像処理ライブラリで、Python Imaging Library (PIL)のforkプロジェクトです。PILは開発が停滞しPython2.7までの対応にとどまっていますが、PillowはPython3.3に対応しています。これからはPILではなくPillowを使用した方がいいと思います。

さて、GISといえば様々なマップ表現を行うことができますが、その中でもシンボルというのはかなり重要な要素です。通常、マップにポイントを配置した場合、このようになります。

f:id:sanvarie:20151225134248p:plain

ただのポイントですね。しかし、このポイントが一体何なのかよくわかりませんね。シンボルファイルを使用してこれを一目でどの目標物か識別できるようにすることができます。

f:id:sanvarie:20151225134341p:plain

ポイントを蛇のマークにしてみました。例えば、コンビニのシンボルはこれ、ガソリンスタンドのシンボルはこれ、といったように定義しておくと直感的にマップがわかりやすくなりますね。実際に業務を行っていると、何十、何百ものシンボルを使用することになりますが、その一覧を作成するのはなかなか難儀な話です(一つ一つが画像になっている場合などは特に)。そこで、対象のシンボルファイルを一覧にした画像を作成してみます。

※今回、元ファイルの拡張子はemfにしますが、通常の画像ファイルなら何でも大丈夫です。

ディレクトリ構成は以下のようなものを前提としています。

f:id:sanvarie:20151225135023p:plain

各フォルダにこのような画像を配置しています。

f:id:sanvarie:20151225164410p:plain

まずは以下サイトでPillowをダウロードしてインストールを行います。
https://pypi.python.org/pypi/Pillow/2.7.0

ソースは以下です。

# -*- coding: utf-8 -*-
from PIL import Image,ImageDraw,ImageFont
import glob

#対象ディレクトリ
dir = 'D:\python\Symbol'

canvas = Image.new('RGB', (1600, 1600), (255, 255, 255))

x = 0
y = 0

#対象ディレクトリをループ
for name in glob.glob('D:\python\Symbol/*'):
    #対象画像でループ
    for n in glob.glob(name+'/*.emf'):
        #キャンバスの高さの限界に達したらスライド
        #対象画像(300×300),テキスト画像(400,100)
        if y != 0 and y % 1600 == 0:
            y = 0
            x += 400

        #対象画像をサムネイルにする
        a_jpg = Image.open(n, 'r')
        canvas.paste(a_jpg, (x, y))
        y += 300

        #対象画像挿し込むテキストの設定
        imgTxt = Image.new('RGB', (400,100),(255, 255, 255))
        draw = ImageDraw.Draw(imgTxt)

        #フォントの設定
        font = ImageFont.truetype('C:\Windows\Fonts\msgothic.ttc', 48)

        #ディレクトリ名からファイル名を取得する
        f = n[n.rfind("\\"):]
        f = f[1:]

        #対象画像のファイル名をサムネイルに挿し込む
        draw.text((0, 0), f.decode('cp932'), font = font, fill='#000000')
        canvas.paste(imgTxt, (x, y))
        y += 100

canvas.save('D:\python\Symbol\sample.jpg', 'JPEG', quality=100, optimize=True)

結果を見てみます。

f:id:sanvarie:20151225141949p:plain

シンボルファイルが一覧になりました。また、おまけでファイル名も画像として挿し込みました。こんな簡単に画像処理ができるなんて、なんて日だ!

以上、本日は「【Pythonで画像処理】Pillowを使用してシンボル一覧を作成する」でした。

Arcpyでフィーチャクラスの座標系を一括で変換する

さて、本日は「Arcpyでフィーチャクラスの座標系を一括で変換する」です。久しぶりにArcGISネタですね。なんとArcMap上ではフィーチャクラスごとにしか座標系の定義をすることができないので、フィーチャクラスが大量に存在している場合、座標系の変換はとても面倒です。そこで、今回は対象のフィーチャクラスの座標系を一括変換するスクリプトを書いてみました。前提条件としては以下です。

・対象のGDBを全て任意のフォルダに格納する

これだけです。とりあえずサンプルとしてこんな感じのフィーチャクラスを作成してみました。

f:id:sanvarie:20151224151941p:plain

座標系はすべて「平面直角座標系第11系(JGD2000)」としてフィーチャクラスを作成しました。

f:id:sanvarie:20151224152727p:plain

この座標系を「平面直角座標系第9系(JGD2000)」に変換したいと思います。

# -*- coding: utf-8 -*-
import arcpy
import os

#対象のフォルダを設定
folder = u"C:\ArcPySample\座標系変換"

#平面直角座標系第9系(JGD2000)にする
projection ="PROJCS['JGD_2000_Japan_Zone_9',GEOGCS['GCS_JGD_2000',DATUM['D_JGD_2000',SPHEROID['GRS_1980',6378137.0,298.257222101]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Transverse_Mercator'],PARAMETER['False_Easting',0.0],PARAMETER['False_Northing',0.0],PARAMETER['Central_Meridian',139.8333333333333],PARAMETER['Scale_Factor',0.9999],PARAMETER['Latitude_Of_Origin',36.0],UNIT['Meter',1.0]]"

#対象フォルダ
dir = os.listdir(folder)

#gdbフォルダをリストに格納
gdb = []

#対象フォルダ内のgdbフォルダをリストに格納
for d in dir:
    if d.find('.gdb') > -1:
        gdb.append(d)

#mdx情報を取得(これを入れるとArcMap上での実行が必要なので、削除してもOK)
mxd = arcpy.mapping.MapDocument(r"CURRENT")

#gdbの数分ループ
for g in gdb:
    ws = os.path.join(folder , g)

    # ワークスペースの設定
    arcpy.env.workspace = ws

    # フィーチャクラスのリストを取得
    fcList = arcpy.ListFeatureClasses()

    # すべてのフィーチャクラスに対してループで実行
    for fc in fcList:
        #投影法の定義
        arcpy.DefineProjection_management(fc,projection)

        #いちいちマップに読み込んでしまうので一旦削除。ってか他にいい処理の仕方があるはず。
        for df in arcpy.mapping.ListDataFrames(mxd):
            for lyr in arcpy.mapping.ListLayers(mxd, "", df):
                arcpy.mapping.RemoveLayer(df, lyr)

結果を見てみます。

f:id:sanvarie:20151224155153p:plain

「平面直角座標系第9系(JGD2000)」に変換することができました。
簡単ではありますが、本日は「Arcpyでフィーチャクラスの座標系を一括で変換する」でした。

PythonでRSSを使用してアンテナサイトを作る

本日はアンテナサイトを作ってみます。アンテナサイトを作るにはRSSを利用します。RSS - Wikipedia
こんなものを用意してくれているなんて便利な世の中ですね。アンテナサイトといえば2ちゃんねるのまとめを対象としたものが多い気がしますが、本ブログはGISをメインとしてますので、GIS関連ニュースのアンテナサイトを作ってみようと思います。まったく需要がないと思うのですが(笑)

まず、RSS解析には「feedparser」というライブラリを使用しますので、pip install feedparser をお願いします。それと「pandas」もインストールをお願いします(←今回の例では使用しましたが、これを使わなくても全然実装可能です)。

また、本ブログでは以下のようなことは考慮していませんので、実際にサイトを立ち上げたりする場合はご注意願います。
・ページ読込時にRSSを取得してる
→本来はRSS取得用のプログラムを定期的に実行してデータを格納しておき、ページ読込時にそのデータを取得した方がいいと思います。そうしないと読込が重くなってしまう可能性があるので。

まずは、サーバーを立ち上げる必要があります。

■cgiserver.py

import CGIHTTPServer
CGIHTTPServer.test()

これだけで開発用サーバーの立ち上げができます。便利ですね!
詳細は以下サイトをご参照ください。
Pythonで学ぶWebアプリ開発のABC みんなのPython Webアプリ編 HTML版(無料) | TRIVIAL TECHNOLOGIES 4 @ats のイクメン日記

RSS対象サイト

以下がRSS対象サイトです。うーん、マニアックですね(笑)

What's New 更新情報 --- GISA 地理情報システム学会
国土地理院ホームページ:新着情報
ESRIジャパン
日本海洋データセンター 新着情報
NPO法人 全国GIS技術研究会

RSS.py

# -*- coding: utf-8 -*-
import feedparser
import time
import pandas as pd

def main():

    #RSSデータ取得
    pds = getRss()

    html_body = u"""<html>
                    <head><meta http-equiv="content-type" content="text/html;charset=utf-8">
                    <link href="../css/sample.css" rel="stylesheet" type="text/css"/>
                    <title>GISアンテナ</title>
                    </head>
                    <body>
                    <div id="pagebody">
                    	<!-- ヘッダ -->
                    	<div id="header"><h1>GISアンテナ</h1></div>
                    	<!-- メインメニュー -->
                    	<ul id="menu">
                    		<li id="menu01"><a href="https://www.google.co.jp/maps" target="_blank">GoogleMap</a></li>
                    		<li id="menu02"><a href="https://www.google.co.jp/intl/ja/earth/" target="_blank">GoogleEarth</a></li>
                    		<li id="menu03"><a href="http://user.numazu-ct.ac.jp/~tsato/webmap/map/gmap.html?data=history" target="_blank">地理院地図</a></li>
                    		<li id="menu04"><a href="http://www.mapion.co.jp/route/" target="_blank">距離測</a></li>
                    		<li id="menu05"><a href="http://www.its-mo.com/" target="_blank">いつもNAVI</a></li>
                    	</ul>
                        <div id="news"><h2>最新記事一覧</h2>
                        <div id="scr"><ul id>
                        """
    cnt = 0
    #最新記事一覧を作成
    for i,row in pds.sort("time",ascending=False).iterrows():
        #表示件数を絞る
        if cnt == 10:
            break
        if row.flg == "0":
            html_body += "<li id =" + "inline" + ">" + row.time + "   " + "<a href=" + row.link + " target" + "=_blank>" + row.title + "</a></li>"
        cnt += 1
    html_body += "</div></div></ul><div id=" + "content" + "><ul id=" + "list" + ">"

    for i,row in pds.iterrows():
        if row.flg == "1":
            html_body += u"<h2><a class=" + "feed-link href=" + row.link + " target" + "=_blank>"  + row.title + "</a></h2>"

        elif row.flg == "0":
            html_body += u"<li id= " + "inline" + "><a href=" + row.link + " target" + "=_blank>"  + row.title + "</a></li>"

    html_body += u"""</ul></div><div id="footer">
                     <address>Copyright (c) GISアンテナ All Rights Reserved.</address></div></div></body></html>
                     """
    print 'Content-type: text/html\n'
    print html_body.encode('utf-8')

def getRss():
    lists =[]

    #対象のRSS
    lists.append("https://www.gisa-japan.org/news/rss.xml")
    lists.append("http://www.gsi.go.jp/index.rdf")
    lists.append("http://www.esrij.com/feed/")
    lists.append("http://www.jodc.go.jp/jodcweb/cgi-bin/public/newsrss1.cgi?lang=ja")
    lists.append("http://www.npo-zgis.or.jp/rss.xml")

    rssList = []
    for i,l in enumerate(lists):
        rss = feedparser.parse(l)
        rssList.append([rss.feed.title,rss.feed.link,"",i,"1"])
        #RSS情報をリストに格納
        for entry in rss.entries:
            #RSS情報をリストに格納
            rssList.append([entry.title,entry.link,time.strftime('%Y/%m/%d %X',entry.updated_parsed),i,"0"])

    #RSS情報をデータフレームに格納
    pds = pd.DataFrame(rssList)
    pds.columns = ["title","link","time","category","flg"]
    return pds

if __name__ == '__main__':
    main()

http://localhost:8000/cgi-bin/RSS.py にアクセスしたときにhtmlを生成するプログラムです。
一応CSSものせます。

■sample.css

/*============================================
全般的なスタイル
============================================*/
* {
	margin:0; padding:0; 	/*全要素のマージン・パディングをリセット*/
	line-height:1.5;	/*全要素の行の高さを1.5倍にする*/
	color:#333333;		/*文字色*/

} 
body {
	background-color:#999999;	/*ページ全体の背景色*/
	text-align:center;		/*IE6以下でセンタリングするための対策*/
}
div#pagebody {
	width:920px; margin:0 auto;	/*内容全体をセンタリング*/
	text-align:left;	        /*テキストの配置を左揃えにする*/
	background-color:#ffffff;	/*内容全体の背景色*/
        
}

/*============================================
ヘッダ
============================================*/
div#header {
	height:77px;	
	background-color:#cccccc;		        /*ヘッダ部分の背景色*/
}
h1 {
	padding:20px 0px 0px 30px;		        /*見出し内容の位置調整*/
	font-family:Arial, Helvetica, sans-serif;	/*フォントの種類*/
}


/*============================================
メインメニュー
============================================*/
ul#menu {
	height:42px; background-color:#eeeeee; font-weight:bold;
}
li#menu01,li#menu02,li#menu03,li#menu04,li#menu05 {
	float:left;			/*リスト項目を横に並べる*/
	display:inline;			/*リスト項目をインライン表示にする*/
	list-style-type:none;		/*リストマーカー無しにする*/
}

ul#menu a {
	display:block;				        /*リンクをブロック表示にする*/
	height:42px; padding-top:4px; text-align:center;
	font-family:Arial, Helvetica, sans-serif;	/*フォントの種類*/
}
/*ボタン01~05にはそれぞれ異なる背景画像を指定する*/
li#menu01 {
	width:164px; height:42px;	                
}
li#menu02 {
	width:156px; height:42px;	
}
li#menu03 {
	width:156px; height:42px;	
}
li#menu04 {
	width:156px; height:42px;	
}
li#menu05 {
	width:164px; height:42px;	
}
ul#menu a {
	display:block;				/*リンクをブロック表示にする*/
	height:42px; padding-top:4px; text-align:center;
	text-decoration:none; 			/*リンクの下線を無くす*/
	font-family:Arial, Helvetica, sans-serif;	/*フォントの種類*/
}


/*============================================
コンテンツ
============================================*/
div#content {
	width:860px; margin:10px 20px 10px 0px;	/*幅の指定と位置調整*/
	float:right;				/*2カラム全体を右寄せにする*/
}

li#inline {
width:800px; margin:10px 20px 10px 0px;	/*幅の指定と位置調整*/
float:right;			        /*2カラム全体を右寄せにする*/
display: inline;
margin: 0 10px;
}

div#scr {
  width: 860px;
  max-height: 200px;
  height: auto !important; /* IE6 max-height対応 */
  height: 100px; /* IE6 max-height対応 */
  overflow-x: hidden; /* 横スクロール非表示 */
  overflow-y: scroll; /* 縦スクロール */
}


div#news {
	width:860px; margin:10px 20px 10px 0px;	/*幅の指定と位置調整*/
	float:right;				/*2カラム全体を右寄せにする*/
}



/*見出し・段落・水平線のスタイル指定*/
h2 {
	font-size:100%; margin-bottom:10px; padding-left:25px;
	font-size:95%; border-bottom:solid 1px #cccccc;
	background-image:url("images/icon.gif");
	background-repeat:no-repeat; background-position:left center;
}

/*============================================
フッタ
============================================*/
div#footer {
	height:42px; text-align:center;
	clear:both;					/*回り込みを解除する*/
	background-color:#cccccc;			/*フッタ部分の背景色*/
}
address {
	font-style:normal;			 /*フォントスタイルを標準にする*/
	font-size:small;			 /*フォントサイズを小さくする*/
	padding:5px 0px 5px 0px;		 /*要素内容の位置調整*/
}

では結果を見てみましょう。

f:id:sanvarie:20151216105139p:plain

f:id:sanvarie:20151216105147p:plain

シンプルですが、アンテナサイトっぽいものができました!でも、地味ですね。もっとデザインセンスを磨かねば。。。
おまけでGoogleMapなどへのリンクも作っておきました。今回はローカルでの作成ですが、もう少しサイトを充実させたうえで公開してみたいですね。

以上、本日は「PythonRSSを使用してアンテナサイトを作る」でした。