SREホールディングス株式会社でデータサイエンティストをしている岡林です。
当社ではAI査定をはじめとする不動産系のプロダクト、AIの開発を行っています。
不動産情報を扱うにあたって、切っても切り離せないのがGIS（地理情報システム）です。
今回は国土交通省の国土数値情報で公開されているデータを使用して、PythonでGISを扱う手順についてご紹介します。
前提
- 使用するライブラリとバージョン
- Python 3.11.0
- geopandas 0.12.2
- fiona 1.8.22
- shapely 2.0.1
- folium 0.14.0
- numpy 1.24.2
※fionaはgeopandasのインストール時に一緒にインストールされますが、現時点の最新バージョン（1.9.1）では、geopandas.read_file()の速度が遅くなるようなので、バージョンを指定してインストールしています。
- 動作環境
- M1 MacBook Pro
- メモリ：16GB
扱う内容
- シェープファイルのロード
- 地図へのプロット
- 距離計算
扱わない内容
- 空間結合
- shapelyや座標参照系の詳細
- 細かい見た目の調整
データのダウンロード
サンプルデータとして地価公示を使用します。
国土数値情報 -> 国土（水・土地） -> 地価公示（ポイント）から現時点の最新版である令和4年のデータ（L01-22_GML.zip）をダウンロードします。
解凍すると以下のようなディレクトリ構造になります。
GIS_test.ipynbにスクリプトを書いていきます。
.
├── L01-22_GML
│ ├── KS-META-L01-22.xml
│ ├── L01-22.dbf
│ ├── L01-22.geojson
│ ├── L01-22.prj
│ ├── L01-22.shp
│ ├── L01-22.shx
│ └── L01-22.xml
└── GIS_test.ipynb
データのロード
geopandasを使用することで、シェープファイル（.shp）を直接ロードしてpandas.DataFrameと同じように扱うことができます。
import geopandas as gpd
gdf = gpd.read_file('./L01-22_GML/L01-22.shp')
列数が140と多いので、使用する列のみに絞ります。併せて列名も変更しておきます。
gdf = gdf[['L01_006', 'L01_024', 'geometry']].rename(columns={'L01_006': 'kouzi_price', 'L01_024': 'address'})
gdf
- L01_006: kouzi_price が公示価格（円/㎡）
- L01_024: address が基準地の所在
- geometry が座標情報
となっています。
geometryの詳細については、shapelyの公式ドキュメントが分かりやすいです。
データの可視化
次にgeometryがどのようなデータなのか可視化してみます。foliumというライブラリを使用して、OpenStreetMap上に描写します。
import folium
fmap = folium.Map(
location=(35.6809591, 139.7673068), # 東京駅の座標
zoom_start=14
)
folium.GeoJson(gdf.to_json()).add_to(fmap)
fmap
folium.Map()でベースとなる地図の座標（中心点）とズームレベルを指定します。
folium.GeoJson(gdf.to_json()).add_to(fmap)でプロットするデータを指定します。
GeoDataFrameはDataFrameと同じように扱うことができるので、例えば地価公示が1,000万円/㎡以上の高額な地点だけプロットしたい場合は以下のように実現できます。
fmap = folium.Map(
location=(35.6809591, 139.7673068), # 東京駅の座標
zoom_start=5
)
folium.GeoJson(gdf[gdf['kouzi_price'] > 10000000].to_json()).add_to(fmap)
fmap
foliumで作成した地図はhtmlで保存することもできます。
fmap.save('test_map.html')
htmlで保存した地図は、ブラウザでも動かして見ることができます。
距離を計算する
2地点間の距離計算として、以下の2パターンのやり方を試してみます。
-
geopandasのdistanceメソッドで計算する - 経度、緯度から近似計算で求める
距離を計算する対象としてポイントを1つ用意しておきます。
from shapely import Point
p = Point(139.7576692, 35.6802117) # 皇居の座標
p = gpd.GeoDataFrame(geometry=[p], crs=6668)
crs=6668は座標参照系 (CRS) をEPSG:6668に指定しています。
地図上にプロットして確認します。
fmap = folium.Map(
location=(35.6809591, 139.7673068), # 東京駅の座標
zoom_start=14
)
folium.GeoJson(p.to_json()).add_to(fmap)
fmap
皇居外苑にピンがプロットされているのが確認できます。
1. geopandasのdistanceメソッドで計算する
distanceメソッドでの距離計算は以下のように行います。
gdf.to_crs(6677).distance(p.to_crs(6677).geometry[0])
0 8.296499e+05
1 8.285137e+05
2 8.291691e+05
3 8.280606e+05
4 8.284030e+05
...
25988 1.566160e+06
25989 1.568438e+06
25990 1.569213e+06
25991 1.567645e+06
25992 1.568706e+06
Length: 25993, dtype: float64
distanceメソッドで距離を計算するためには、座標参照系を地理座標系から平面直角座標系に変換する必要があります。to_crs(6677)でEPSG:6677に変換しています。また、distanceメソッドで配列を指定すると、1対1の行単位での距離計算をしてしまうため、geometry[0]で1つのgeometryを指定しています。（詳細は公式ドキュメント）
2. 経度、緯度から近似計算で求める
経度、緯度から距離を算出する近似式はいくつもありますが、比較的精度が高く計算式もシンプルな球面三角法を用いた以下の式を使用します。
import numpy as np
def calc_distance(lat1, lng1, lat2, lng2):
R = 6371. # 地球の平均半径
# 度数法からラジアンに変換
lat1 = np.deg2rad(lat1)
lng1 = np.deg2rad(lng1)
lat2 = np.deg2rad(lat2)
lng2 = np.deg2rad(lng2)
d = R * np.arccos(
np.cos(lat1) * np.cos(lat2) * np.cos(lng1 - lng2)
+ np.sin(lat1) * np.sin(lat2)
)
return d
calc_distance関数で距離を計算します。calc_distance関数は配列を入力として受け取れるので、lat1, lng1にはgdfの緯度、経度の配列をそれぞれ指定します。
lat1 = gdf.centroid.y
lng1 = gdf.centroid.x
lat2 = p.geometry.y[0]
lng2 = p.geometry.x[0]
calc_distance(lat1, lng1, lat2, lng2)
0 830.885665
1 829.743913
2 830.399170
3 829.288630
4 829.637172
...
25988 1557.191266
25989 1559.487333
25990 1560.244389
25991 1558.674162
25992 1559.743072
Length: 25993, dtype: float64
これで経度、緯度から距離を計算することができました。
近似計算の場合は座標参照系の変換が必要なく、どの直交座標系を使用するかを迷う必要がないため、より使用しやすいかと思います。
おまけ
座標参照系の補足
GeoDataFrame.crsでGeoDataFrameに設定されているcrsを確認することができます。
# 地価公示データのcrs
gdf.crs
<Geographic 2D CRS: EPSG:4612>
Name: JGD2000
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: Japan - onshore and offshore.
- bounds: (122.38, 17.09, 157.65, 46.05)
Datum: Japanese Geodetic Datum 2000
- Ellipsoid: GRS 1980
- Prime Meridian: Greenwich
# 皇居にプロットしたポイントのcrs
p.crs
<Geographic 2D CRS: EPSG:6668>
Name: JGD2011
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: Japan - onshore and offshore.
- bounds: (122.38, 17.09, 157.65, 46.05)
Datum: Japanese Geodetic Datum 2011
- Ellipsoid: GRS 1980
- Prime Meridian: Greenwich
JGD2000とJGD2011は異なる座標参照系ですが、実用上は誤差が小さいため同列に扱っています。また、可視化で仕様したOpenStreetMapはWGS84という別の座標参照系ですが、こちらも実用上の誤差が小さいため、同列に扱っています。
distanceメソッドで距離計算するために変換したEPSG:6677は、東京周辺の距離を正確に表現できる平面直角座標系です。（詳細は国土地理院の説明ページ）
距離計算の精度
国土地理院の距離計算と比較してみます。
座標の設定
- 東京（皇居）: 139.75767, 35.68021
- 東京（近傍点）: 139.76262, 35.67752
- 沖縄: 124.14933, 24.34700
- 北海道: 141.67194, 45.42195
東京（皇居）からの距離（単位はkm）
| 地点 | 国土地理院 | distanceメソッド | 誤差率 | 近似計算 | 誤差率 |
|---|---|---|---|---|---|
| 東京（近傍点） | 0.538410 | 0.538485 | 0.0139% | 0.538058 | -0.0654% |
| 沖縄 | 1,956.375739 | 1,976.928000 | 1.0505% | 1,956.610837 | 0.0120% |
| 北海道 | 1,093.769116 | 1,093.749000 | -0.0018% | 1,095.147157 | 0.1260% |
distanceメソッドによる距離計算は沖縄との距離で誤差が1%を超えていますが、概ね1%には収まりそうです。
計算速度の違い
distanceメソッドとnumpyを使用した近似計算では、データ数が多い場合の計算速度に大きな差があります。25,993行の地価公示データに対して、1地点との距離計算を10,000回ループした時の時間は以下の通りです。（M1 MacBook Proで計測）
- distanceメソッド: 67s
- 近似計算: 6s
10倍程度の差があるため、データ数が多い場合には特に近似計算を使用した方が良さそうです。
実行したスクリプトは以下の通り。
distanceメソッド
gdf_ = gdf.to_crs(6677)
p_ = p.to_crs(6677)
for _ in range(10000):
d = gdf_.distance(p_.geometry[0])
近似計算
lat1 = gdf.geometry.y
lng1 = gdf.geometry.x
lat2 = p.geometry.y[0]
lng2 = p.geometry.x[0]
for _ in range(10000):
d = calc_distance(lat1, lng1, lat2, lng2)
最後に
GIS（地理情報システム）の取り扱いと距離計算について紹介しました。geopandasによってPythonでも比較的簡単にGISを扱えます。一方でGIS特有の座標参照系に対する理解や、計算によるズレが生じる点について注意する必要があります。１つ１つプロットしたり計算して、想定通りの結果になっているか、誤差が許容範囲内かを確認することが重要です。
Discussion