Python × GIS の基礎（その1） #Python

なぜPythonを使う？

GISソフトだけで大量の空間データを処理するには時間がかかり、研究を進めていく上でコード処理の必要を感じていました。またQGISに存在するプラグインではネットワーク解析や図形単純化を十分に行えず, Pythonの豊富なライブラリは有効な手段となります。

今回はヘルシンキ大学が提供しているAutomating GIS processesを勉強の素材として用います。
7週構成になっており、GeoPnadasやmatplotlibを用いた可視化について特に充実して学ぶことができます。また各週に受講者向けの宿題(練習問題)が付けられており, 本記事では主にその解答例を載せます。
その1：Week1-Week2
その2：Week3-Week4
その3：Week5-Week6(Week7は課題なし)
その4(作成中)：Final Assignment

下記は週毎の教材の内容の一部です。

週	内容
Week1	ポイント, ライン, ポリゴンの描画
Week2	Geopandasの導入, CRS設定, 距離計算
Week3	空間結合, 近傍点分析
Week4	空間結合2, データ階層化
Week5	静的・動的な地図の描画,
Week6	OSMの分析, ネットワーク解析
Week7	QGISプラグイン

Week 1

1-1 ポイント, ライン, ポリゴンを作る関数の作成

問題の条件として関数が正常に機能するよういくつかの要件が求められています。

from shapely.geometry import Point, LineString, Polygon

#ポイント
def create_point_geom(x, y):
    point = Point(x, y)
    return point

#ライン リストを引数に取ります
def create_line_geom(points):
    assert type(points)=="list", "Input should be a list!"
    assert len(points)>=2, "LineString object requires at least two Points!"
    line = LineString([points[0], points[1]])
    return line

#ポリゴン リストを引数に取ります
def create_poly_geom(coords):
    assert type(coords) is list, "Input should be a list!"
    assert len(coords)>=3, "Polygon object requires at least three Points!"
    for i in coords:
        assert type(i) is tuple, "All list values should be coordinate tuples!"
    poly = Polygon(coords)
    return poly

1-2 重心, ポリゴン面積, 距離の取得

assert文は本質的ではないので一部省略しています。

#重心　引数はポイントorラインorポリゴン
def get_centroid(gem):
    assert type(gem) == 'shapely',  "Input should be a Shapely geometry!"
    assert gem.geom_type in ['Point', 'LineString', 'Polygon'], "Input should be a Shapely geometry!"
    centroid = gem.centroid
    return centroid
#面積 
def get_area(poly):
    return poly.area
#距離
def get_length(geom):
    if geom.geom_type == 'LineString':
        return geom.length
    elif geom.geom_type == 'Polygon':
        return geom.exterior.length

1-3 テキストデータの取得と表示

ヘルシンキの移動データを用いた簡単なデータ加工の練習問題です。
from_x, from_y は出発点, to_x, to_y は到着点を示しており, これらのデータから総移動距離を求めるというのが目標です。

import pandas as pd

#データ読み込み 
data = pd.read_table("data/travelTimes_2015_Helsinki.txt", sep=";",)
data = data[['from_x','from_y', 'to_x', 'to_y', 'total_route_time',]]
#座標データの作成
orig_points = []
dest_points = []
from shapely.geometry import Point
for index, row in data.iterrows():
    orig = Point(row['from_x'], row['from_y'])
    dest = Point(row['to_x'], row['to_y'])
    orig_points.append(orig)
    dest_points.append(dest)
#移動ラインの作成
from shapely.geometry import LineString
lines = []
for orig, dest in zip(orig_points, dest_points):
    line = LineString([orig, dest])
    lines.append(line)
#総移動距離の取得
total_length = 0
for line in lines:
    total_length += line.length

Week2

Pandasの拡張形であるGeopandasというライブラリを学びます。Pandas同様非常に使いやすく、空間操作も簡単に行えます。

2-1 座標データのGeoDataFrameへの変換

import geopandas as gpd
import matplotlib.pyplot as plt
from shapely.geometry import Polygon

#使う座標データ
longitudes = [29.99671173095703, 31.58196258544922, 27.738052368164062, 26.50013542175293, 26.652359008789062, 25.921663284301758, 22.90027618408203, 23.257217407226562,
           23.335693359375, 22.87444305419922, 23.08465003967285, 22.565473556518555, 21.452774047851562, 21.66388702392578, 21.065969467163086, 21.67659568786621,
           21.496871948242188, 22.339998245239258, 22.288192749023438, 24.539581298828125, 25.444232940673828, 25.303749084472656, 24.669166564941406, 24.689163208007812,
           24.174999237060547, 23.68471908569336, 24.000761032104492, 23.57332992553711, 23.76513671875, 23.430830001831055, 23.6597900390625, 20.580928802490234, 21.320831298828125,
           22.398330688476562, 23.97638702392578, 24.934917449951172, 25.7611083984375, 25.95930290222168, 26.476804733276367, 27.91069221496582, 29.1027774810791, 29.29846954345703,
           28.4355525970459, 28.817358016967773, 28.459857940673828, 30.028610229492188, 29.075136184692383, 30.13492774963379, 29.818885803222656, 29.640830993652344, 30.57735824584961,
           29.99671173095703]
latitudes = [63.748023986816406, 62.90789794921875, 60.511383056640625, 60.44499588012695, 60.646385192871094, 60.243743896484375, 59.806800842285156, 59.91944122314453,
           60.02395248413086, 60.14555358886719, 60.3452033996582, 60.211936950683594, 60.56249237060547, 61.54027557373047, 62.59798049926758, 63.02013397216797,
           63.20353698730469, 63.27652359008789, 63.525691986083984, 64.79915618896484, 64.9533920288086, 65.51513671875, 65.65470886230469, 65.89610290527344, 65.79151916503906,
           66.26332092285156, 66.80228424072266, 67.1570053100586, 67.4168701171875, 67.47978210449219, 67.94589233398438, 69.060302734375, 69.32611083984375, 68.71110534667969,
           68.83248901367188, 68.580810546875, 68.98916625976562, 69.68568420410156, 69.9363784790039, 70.08860778808594, 69.70597076416016, 69.48533630371094, 68.90263366699219,
           68.84700012207031, 68.53485107421875, 67.69471740722656, 66.90360260009766, 65.70887756347656, 65.6533203125, 64.92096710205078, 64.22373962402344, 63.748023986816406]
#ポリゴン化
coordpairs = list(zip(longitudes, latitudes))
poly = Polygon(coordpairs)
#GeoDataFrameの作成
geo = gpd.GeoDataFrame(index=[0], columns=['geometry'])
geo['geometry'] = poly
#図示
import matplotlib.pyplot as plt
geo.plot()
#shpファイルの保存
fp = 'polygon.shp'
geo.to_file(fp)

2-2 ｃｓｖに格納された座標データのプロット

憎らしいCRSが登場します。apply関数はpandasとの相性が良く, 短いコードで加工ができます。

#ライブラリ
import pandas as pd
from shapely.geometry import Point, LineString, Polygon
from pyproj import CRS
import matplotlib.pyplot as plt
#csv読み込み
data = pd.read_csv('data/some_posts.csv')
#ポイントデータ作成
data = pd.read_csv('data/some_posts.csv')
make_point = lambda row:Point(row['lat'],row['lon'])
data['geometry'] = data.apply(make_point, axis=1)

#GeoDataFrameへの変換 geometryと座標系を指定する必要があります.
geo = gpd.GeoDataFrame(data, geometry='geometry',crs=CRS.from_epsg(4326).to_wkt())

#プロット
geo.plot()

2-3 移動距離の計算

SNSの投稿の座標データから各ユーザーの移動距離を計算します。
座標系の変換が登場しますが、以下の点に注意が必要です。CRSはどこでも厄介。。。

・座標系の定義(ex. data.crs=CRS.from_espg(4276) )だけではGeoDataFrameのデータは変換しないため, 座標系の変換(ex. data=data.to_crs(epsg=4276) ) が必要
・GeoDataFrameの座標系が定義されていない場合は, 座標系変換の前に座標系定義を行う


import geopandas as gpd
import pandas as pd
from pyproj import CRS
from shapely.geometry import Point, LineString, Polygon
#データ読み込み
data = gpd.read_file('Kruger_posts.shp')
#座標系変換
data = data.to_crs(epsg=32735)
#idごとに分類して移動ラインを作成(投稿が一つのみの場合移動していないため, None値を入れる)
grouped = data.groupby('userid')
movements = gpd.GeoDataFrame(columns=['userid', 'geometry'])
for key, group in grouped:
    group = group.sort_values('timestamp')
    if len(group['geometry'])>=2:
        line = (LineString(list(group['geometry'])))
    else:
        line=None
    movements.at[count, 'userid'] = key
    movements.at[count, 'geometry'] = line
movements.crs = CRS.from_epsg(32735)
#移動距離の計算
def cal_distance(x):
    if x['geometry'] is None:
        return None
    else:
        return x['geometry'].length
movements['distance'] = movements.apply(cal_distance, axis=1)
#移動距離の平均, 最大値, 最小値
print(mean(movements['distance'])) #138871.14194459998
print(max(movements['distance'].dropna())) #8457917.497356484
print(min(movements['distance'])) #0.0