Python備忘録　OpenStreetMapとgeocoder、さらにChatGPTの活用（２）

- 3月 30, 2023

　「Python備忘録　OpenStreetMapとgeocoder、さらにChatGPTの活用（１）」で作成した四国88箇所霊場の位置情報を利用して、PythonGeopandasほかの理解を深めていきたいと思います。

　まず、位置情報を読み込みます。

#Google ColabでGoogle Driveを使用するときのおまじない
from google.colab import drive
drive.mount('/content/drive')

!pip install geopandas
import geopandas as gpd
from shapely.geometry import Point

#前述のCodeで作成・保管した四国88箇所霊場のGeoPackage形式
#位置データを読み込み
gdf_88temples = gpd.read_file('/content/drive/MyDrive/Colab Notebooks/202303_Shikoku88Data/Shikoku88Temples_lnglat.gpkg')

gdf_88temples.head()

　出力はこちらです。

	Name	Address	geometry
0	霊山寺　徳島県	徳島県鳴門市大麻町　板東	POINT (134.51733 34.18415)
1	極楽寺　徳島県	徳島県鳴門市大麻町桧字ダンノ上	POINT (139.52867 35.30921)
2	金泉寺　徳島県	徳島県板野郡板野町大寺字亀山下	POINT (134.49020 34.15297)
3	大日寺　徳島県	徳島県板野郡板野町黒谷字居内	POINT (134.46259 34.03811)
4	地蔵寺　徳島県	徳島県板野郡板野町羅漢字林東	POINT (134.58525 34.07603)

　ここでこけたるわけには行きません。

　次に、各霊場の周囲5kmの範囲をbuffer()メソッドを利用して定義します。ちなみにChatGPTさんにbuffer()メソッドについて尋ねると、以下のような答えが返ってきました。

2 / 2

geopandas の　buffer()の使い方を教えて

geopandasのbuffer()メソッドは、地理空間オブジェクトに対してバッファーを作成するために使用されます。バッファーは、地理空間オブジェクトの周りに指定された距離の円または多角形を作成することができます。

以下のコードは、buffer()メソッドの使用例を示しています。

import geopandas as gpd
from shapely.geometry import Point

# ポイントオブジェクトを作成する
point = Point(0, 0)

# GeoDataFrameを作成する
gdf = gpd.GeoDataFrame(geometry=[point])

# バッファーを作成する
buffer_distance = 1.0  # バッファーの距離（メートル）
gdf_buffer = gdf.buffer(buffer_distance)

# GeoDataFrameをプロットする
ax = gdf.plot(color='blue')
gdf_buffer.plot(ax=ax, color='red', alpha=0.5)

この例では、Point()関数を使用して、(0, 0)の座標を持つポイントオブジェクトを作成します。次に、このポイントオブジェクトをGeoDataFrameに追加します。

buffer_distance変数で指定された距離のバッファーを作成するには、gdf.buffer()メソッドを使用します。このメソッドは、GeoDataFrameの各要素に対してバッファーを作成するため、結果は新しいGeoDataFrameオブジェクトとして返されます。この例では、バッファーの距離は1メートルです。

最後に、GeoDataFrame.plot()メソッドを使用して、元のポイントとバッファーをプロットします。プロットは、青色のポイントと、赤色の透過性のあるバッファーで表示されます。

バッファーの形状を設定するために、buffer()メソッドにいくつかのオプション引数を渡すことができます。たとえば、buffer(cap_style=3)のように、cap_style引数を使用して、多角形の端の形状を設定することができます。詳細については、geopandasのドキュメントを参照してください。

　ふむふむ。

　さらに、座標表示が緯度経度の場合のbuffer()メソッドの距離の単位を尋ねると「度なので、メートル単位を使うには適切な単位変換が重要です」と丁寧に答えてくれました。

　なので、実は世界測地系の緯度経度は距離や面積を計算するのが苦手で、そのままの座標系で面積を求めたりすると、わけのわからない単位・桁の数値を返してきます。これは、以前のPhtyon備忘録でも触れました。

　というわけで、ある程度の精度を求めて座標系変換を行います。どの座標系がよいかChatGPTさまに尋ねてみました。

JGD2011から距離や面積を高い精度で求めることができる座標系に変換するにはどうすればよいですか

JGD2011測地系は、日本国内で広く使用されている測地系の1つであり、測地系によって定義された座標系であるため、そのままでは距離や面積の計算には向いていません。より高い精度で距離や面積を求めるためには、JGD2011座標系に変換する必要があります。
JGD2011座標系は、主に以下の2つのEPSGコードで表されます。
EPSG:6669 - JGD2011 / UTM zone 54N
EPSG:6670 - JGD2011 / UTM zone 55N
UTM（ユニバーサル横メルカトル図法）は、地球を6度幅の帯に分割し、それぞれに座標系を割り当てたものです。UTM座標系は、距離や面積の計算に適した座標系であり、JGD2011座標系もUTM座標系の一種です。（以下、省略）

EPSG:6669 - JGD2011 / UTM zone 54N
EPSG:6670 - JGD2011 / UTM zone 55N　のzone54N, zone55Nはどういう意味ですか。

EPSG:6669とEPSG:6670は、UTM座標系であり、それぞれ日本のJGD2011測地系を基準としています。UTM座標系は、地球を60個の緯度帯に分割し、それぞれの緯度帯で直交座標系を定義する方法です。各緯度帯は6度の幅を持ち、AからZまでのアルファベットで表されます。
JGD2011 / UTM zone 54Nは、日本の東経129度から132度の範囲に相当し、緯度帯は54Nに相当します。同様に、JGD2011 / UTM zone 55Nは、日本の東経132度から135度の範囲に相当し、緯度帯は55Nに相当します。
つまり、EPSG:6669とEPSG:6670は、日本の座標系としてUTM座標系を採用しており、それぞれ東経129度から132度、132度から135度の範囲を対象としています。

　四国はおおむね東経132から135度ということで、EPSG:6670 - JGD2011 / UTM zone 55Nを選択して、座標変換を試みます。

#Google ColabでGoogle Driveを使用するときのおまじない
from google.colab import drive
drive.mount('/content/drive')

!pip install geopandas
import geopandas as gpd
from shapely.geometry import Point

のあとの、

#前述のCodeで作成・保管した四国88箇所霊場のGeoPackage形式
#位置データを読み込み, EPSG:6668(JGD2011測地系)
gdf_88temples = gpd.read_file('/content/drive/MyDrive/Colab Notebooks/202303_Shikoku88Data/Shikoku88Temples_lnglat.gpkg')

#JGD2011測地系座標(EPSG:6668)をUTM zone 55N(EPSG:6670…日本の東経132度から135度に適用)に変換
gdf_88temples_UTM55N = gdf_88temples.to_crs(epsg=6670)

gdf_88temples_UTM55N.head()

です。結果は以下の通り。

	Name	Address	geometry
0	霊山寺　徳島県	徳島県鳴門市大麻町　板東	POINT (324288.224 136923.912)
1	極楽寺　徳島県	徳島県鳴門市大麻町桧字ダンノ上	POINT (776502.671 289673.529)
2	金泉寺　徳島県	徳島県板野郡板野町大寺字亀山下	POINT (321904.527 133376.613)
3	大日寺　徳島県	徳島県板野郡板野町黒谷字居内	POINT (319788.536 120541.401)
4	地蔵寺　徳島県	徳島県板野郡板野町羅漢字林東	POINT (330974.880 125141.716)

POINTの数値が緯度経度から変換されています。

チェックのために図化してみます。

#buffer()メソッドを使って各霊場を中心にした5kmの範囲を定義
buffer_distance = 10_000 #1km=1,000m
gdf_88temples_UTM55N_buffer = gdf_88temples_UTM55N.buffer(buffer_distance)

gdf_88temples_UTM55N_buffer.plot()
gdf_88temples_UTM55N.plot()

結果はこちら。

　なんかグラフ真ん中上方と右上あたりの変なところに２グループほどの点がありますねぇ。元のデータにも同じような点が存在し、東経が139度というオーダーになっているので明らかな間違いですね。

　同名の寺院が四国外に存在するのでしょう。もとのデータに立ち返って、国土地理院HPより四国の東端(徳島県阿南市伊島→東経134度49分20秒≒134.82度未満)、西端（愛媛県伊方町佐田岬→東経132度00分45秒≒132度）を参考に、霊場の名称から求めた緯度経度が132度以上134.82以下ならOKで、その範囲を逸脱していたら番地なしの代表住所から求めた緯度経度を求めるというルールに変更します。以下、修正したコード（修正箇所は#→2023/04/01修正から#←2023/04/01修正までです。）

#四国遍路巡礼88箇所の寺社名、住所一覧表の読み込み
# Read the CSV file into a pandas DataFrame
df_temples = pd.read_csv(
    '/content/drive/MyDrive/Colab Notebooks/202303_Shikoku88Data/20230327_四国四県GIS情報/20230325_Shikoku88ShrineTable_v01.csv')

#霊場の名称、所在地、番地無し住所、緯度経度を格納する
#GeoDataFrameを定義
columns = ['Name', 'Address', 'geometry']
gdf_temples = gpd.GeoDataFrame(columns = columns)

#名称＋県名を用いてOpenStreetMapで経度緯度を検索
for index,row in df_temples.iterrows():
    ret = geocoder.osm(row["名称"] + '日本', timeout=5.0)

    # 所在地から番地を取り除き'番地無し住所'欄に入力
    pattern_address = r'\d{1,2}[-\d]{0,4}\s*[番号号]?(?![\d-])'
    address = row["所　在　地"]
    match = re.search(pattern_address, address)
        
    if match:
        row["番地なし住所"] = re.sub(pattern_address, '', address).strip()

    # '名称'から経度緯度を特定できない場合には'番地無し住所'から当該市町村
    # の代表地点緯度経度を出力する
    while ret.latlng == None:
        ret = geocoder.osm(row["番地なし住所"], timeout=5.0)
        if ret.latlng == None:
            # '番地なし住所'の右端の文字列を削除
            row["番地なし住所"] = row["番地なし住所"][:-1]
        else:
            break

#→2023/04/01_追加
    lat, lng = ret.latlng
    while ret.latlng == None or lng < 132 or lng > 134.82:
        print("**計算中**", row["名称"], row["番地なし住所"], "lng: ", lng, "lat: ", lat)
        ret = geocoder.osm(row["番地なし住所"], timeout=5.0)
        if ret.latlng == None:
            # '番地なし住所'の右端の文字列を削除
            row["番地なし住所"] = row["番地なし住所"][:-1]
        elif ret.latlng != None:
            #経度緯度に新たに取得した"番地なし住所"の数値を代入
            lat, lng = ret.latlng
        else:
            break
#←2023/04/01_追加

    gdf_temples.loc[len(gdf_temples)] = [row["名称"], row["番地なし住所"], Point(lng, lat)]
    print(row["名称"], ", ", row["番地なし住所"], "lng: ", lng, "lat: ", lat)

#GeoPackege形式で保存
gdf_temples.to_file('/content/drive/MyDrive/Colab Notebooks/202303_Shikoku88Data/Shikoku88Temples_lnglat.gpkg', driver = 'GPKG')

計算過程を含めて、出力です。**計算中**と表示されている部分が、今回の再計算対象です。徳島県の極楽寺、安楽寺、十楽寺、熊谷寺、法輪寺、常楽寺、国分寺、井戸寺、鶴林寺、薬王寺が再計算の対象になりました。不思議と徳島県に集まっているのはただの偶然でしょうか？

霊山寺　徳島県 ,  徳島県 鳴門市大麻町　板東 lng:  134.5173311 lat:  34.1841514
**計算中** 極楽寺　徳島県 徳島県 鳴門市大麻町桧字ダンノ上 lng:  139.5286725 lat:  35.3092118
**計算中** 極楽寺　徳島県 徳島県 鳴門市大麻町桧字ダンノ lng:  139.5286725 lat:  35.3092118
**計算中** 極楽寺　徳島県 徳島県 鳴門市大麻町桧字ダン lng:  139.5286725 lat:  35.3092118
**計算中** 極楽寺　徳島県 徳島県 鳴門市大麻町桧字ダ lng:  139.5286725 lat:  35.3092118
**計算中** 極楽寺　徳島県 徳島県 鳴門市大麻町桧字 lng:  139.5286725 lat:  35.3092118
**計算中** 極楽寺　徳島県 徳島県 鳴門市大麻町桧 lng:  139.5286725 lat:  35.3092118
極楽寺　徳島県 ,  徳島県 鳴門市大麻町桧 lng:  134.4974121 lat:  34.1669586
金泉寺　徳島県 ,  徳島県 板野郡板野町大寺字亀山下 lng:  134.4902033 lat:  34.1529667
大日寺　徳島県 ,  徳島県 板野郡板野町黒谷字居内 lng:  134.46258518184646 lat:  34.038105599999994
地蔵寺　徳島県 ,  徳島県 板野郡板野町羅漢字林東 lng:  134.58524816906916 lat:  34.07603225
**計算中** 安楽寺　徳島県 徳島県 板野郡上板町引野字寺ノ西北 lng:  136.033438 lat:  35.4478005
**計算中** 安楽寺　徳島県 徳島県 板野郡上板町引野字寺ノ西 lng:  136.033438 lat:  35.4478005
**計算中** 安楽寺　徳島県 徳島県 板野郡上板町引野字寺ノ lng:  136.033438 lat:  35.4478005
**計算中** 安楽寺　徳島県 徳島県 板野郡上板町引野字寺 lng:  136.033438 lat:  35.4478005
**計算中** 安楽寺　徳島県 徳島県 板野郡上板町引野字 lng:  136.033438 lat:  35.4478005
**計算中** 安楽寺　徳島県 徳島県 板野郡上板町引野 lng:  136.033438 lat:  35.4478005
安楽寺　徳島県 ,  徳島県 板野郡上板町引野 lng:  134.3903269 lat:  34.1231987
**計算中** 十楽寺　徳島県 徳島県 阿波市土成町　高尾字法教田 lng:  135.8763788 lat:  35.0638987
**計算中** 十楽寺　徳島県 徳島県 阿波市土成町　高尾字法教 lng:  135.8763788 lat:  35.0638987
**計算中** 十楽寺　徳島県 徳島県 阿波市土成町　高尾字法 lng:  135.8763788 lat:  35.0638987
**計算中** 十楽寺　徳島県 徳島県 阿波市土成町　高尾字 lng:  135.8763788 lat:  35.0638987
**計算中** 十楽寺　徳島県 徳島県 阿波市土成町　高尾 lng:  135.8763788 lat:  35.0638987
十楽寺　徳島県 ,  徳島県 阿波市土成町　高尾 lng:  134.3725236 lat:  34.1270227
**計算中** 熊谷寺　徳島県 徳島県 阿波市土成町　土成字前田 lng:  139.3834687194688 lat:  36.146302750000004
**計算中** 熊谷寺　徳島県 徳島県 阿波市土成町　土成字前 lng:  139.3834687194688 lat:  36.146302750000004
**計算中** 熊谷寺　徳島県 徳島県 阿波市土成町　土成字 lng:  139.3834687194688 lat:  36.146302750000004
**計算中** 熊谷寺　徳島県 徳島県 阿波市土成町　土成 lng:  139.3834687194688 lat:  36.146302750000004
熊谷寺　徳島県 ,  徳島県 阿波市土成町　土成 lng:  134.3337537 lat:  34.1169467
**計算中** 法輪寺　徳島県 徳島県 阿波市土成町　土成字田中 lng:  135.677215 lat:  35.010001
**計算中** 法輪寺　徳島県 徳島県 阿波市土成町　土成字田 lng:  135.677215 lat:  35.010001
**計算中** 法輪寺　徳島県 徳島県 阿波市土成町　土成字 lng:  135.677215 lat:  35.010001
**計算中** 法輪寺　徳島県 徳島県 阿波市土成町　土成 lng:  135.677215 lat:  35.010001
法輪寺　徳島県 ,  徳島県 阿波市土成町　土成 lng:  134.3337537 lat:  34.1169467
切幡寺　徳島県 ,  徳島県 阿波市市場町　切幡字観音 lng:  133.2110846 lat:  34.4006552
藤井寺　徳島県 ,  徳島県 吉野川市鴨島町飯尾 lng:  133.73728775270988 lat:  34.57373475
焼山寺　徳島県 ,  徳島県 名西郡神山町下分字地中 lng:  134.310411 lat:  33.985057
大日寺　徳島県 ,  徳島県 徳島市一宮町西丁 lng:  134.46258518184646 lat:  34.038105599999994
**計算中** 常楽寺　徳島県 徳島県 徳島市国府町延命 lng:  139.5407424 lat:  35.3482067
常楽寺　徳島県 ,  徳島県 徳島市国府町延命 lng:  134.4807317 lat:  34.0496833
**計算中** 国分寺　徳島県 徳島県 徳島市国府町矢野－ lng:  135.78861594346432 lat:  35.48001445
国分寺　徳島県 ,  徳島県 徳島市国府町矢野－ lng:  134.4746384 lat:  34.0578445
観音寺　徳島県 ,  徳島県 徳島市国府町字観音寺 lng:  134.54505145000002 lat:  34.0612512
**計算中** 井戸寺　徳島県 徳島県 徳島市国府町井戸－ lng:  135.21070011427173 lat:  34.70657555
井戸寺　徳島県 ,  徳島県 徳島市国府町井戸－ lng:  134.4843104 lat:  34.0831464
恩山寺　徳島県 ,  徳島県 小松島市田野町恩山寺谷 lng:  134.5778614 lat:  33.9860981
立江寺　徳島県 ,  徳島県 小松島市立江町 lng:  134.6071878 lat:  33.9653211
**計算中** 鶴林寺　徳島県 徳島県 勝浦郡勝浦町鶴山 lng:  139.4629112 lat:  35.5012402
**計算中** 鶴林寺　徳島県 徳島県 勝浦郡勝浦町鶴 lng:  139.4629112 lat:  35.5012402
**計算中** 鶴林寺　徳島県 徳島県 勝浦郡勝浦町 lng:  139.4629112 lat:  35.5012402
鶴林寺　徳島県 ,  徳島県 勝浦郡勝浦町 lng:  134.5049636 lat:  33.9257826
太竜寺　徳島県 ,  徳島県 阿南市加茂町龍山 lng:  134.55370197224073 lat:  34.0720771
平等寺　徳島県 ,  徳島県 阿南市新野町秋山 lng:  134.58250779182066 lat:  33.851728949999995
**計算中** 薬王寺　徳島県 徳島県 海部郡日和佐町奥河内寺前－ lng:  139.5342141 lat:  35.5603633
**計算中** 薬王寺　徳島県 徳島県 海部郡日和佐町奥河内寺前 lng:  139.5342141 lat:  35.5603633
**計算中** 薬王寺　徳島県 徳島県 海部郡日和佐町奥河内寺 lng:  139.5342141 lat:  35.5603633
**計算中** 薬王寺　徳島県 徳島県 海部郡日和佐町奥河内 lng:  139.5342141 lat:  35.5603633
**計算中** 薬王寺　徳島県 徳島県 海部郡日和佐町奥河 lng:  139.5342141 lat:  35.5603633
**計算中** 薬王寺　徳島県 徳島県 海部郡日和佐町奥 lng:  139.5342141 lat:  35.5603633
**計算中** 薬王寺　徳島県 徳島県 海部郡日和佐町 lng:  139.5342141 lat:  35.5603633
**計算中** 薬王寺　徳島県 徳島県 海部郡日和佐 lng:  139.5342141 lat:  35.5603633
薬王寺　徳島県 ,  徳島県 海部郡日和佐 lng:  134.5307127 lat:  33.7291961
最御崎寺　高知県 ,  高知県 室戸市室戸岬町－ lng:  134.1759439 lat:  33.2490846
津照寺　高知県 ,  高知県 室戸市室津 lng:  133.53528266042073 lat:  33.565381900000006
金剛頂寺　高知県 ,  高知県 室戸市 lng:  134.1522806 lat:  33.2898523
神峰寺　高知県 ,  高知県 安芸郡安田町唐の浜 lng:  134.008901 lat:  33.429891
大日寺　高知県 ,  高知県 香美郡野市町母代寺 lng:  133.70523012135828 lat:  33.57769055
国分寺　高知県 ,  高知県 南国市国分 lng:  133.0242384 lat:  34.0252747
善楽寺　高知県 ,  高知県 高知市一宮 lng:  133.573539 lat:  33.596092
竹林寺　高知県 ,  高知県 高知市五台山 lng:  133.585875 lat:  33.54279
禅師峰寺　高知県 ,  高知県 南国市十市 lng:  133.624401 lat:  33.532119
雪蹊寺　高知県 ,  高知県 高知市長浜－ lng:  133.54175613260276 lat:  33.5437543
種間寺　高知県 ,  高知県 吾川郡春野町秋山 lng:  134.008901 lat:  33.429891
清滝寺　高知県 ,  高知県 土佐市高岡町丁 lng:  133.4508805076673 lat:  33.42593935
青竜寺　高知県 ,  高知県 土佐市宇佐町　竜 lng:  133.4508805076673 lat:  33.42593935
岩本寺　高知県 ,  高知県 高岡郡四万十町茂串－ lng:  133.13466725867693 lat:  33.2081737
金剛福寺　高知県 ,  高知県 土佐清水市足摺岬－ lng:  133.0100828 lat:  32.728457
延光寺　高知県 ,  高知県 宿毛市平田町中山 lng:  132.7754017 lat:  32.9601995
観自在寺　愛媛県 ,  愛媛県 南宇和郡愛南町御荘平城－ lng:  132.5617799 lat:  32.9610057
龍光寺　愛媛県 ,  愛媛県 宇和島市三間町戸雁 lng:  132.601247 lat:  33.295845
佛木寺　愛媛県 ,  愛媛県 宇和島市三 lng:  132.573381 lat:  33.235656
明石寺　愛媛県 ,  愛媛県 西予市宇和町明石 lng:  132.7397885 lat:  33.8916797
大宝寺　愛媛県 ,  愛媛県 上浮穴郡久万高原町菅生二番耕地 lng:  132.74221817531122 lat:  33.84169675
岩屋寺　愛媛県 ,  愛媛県 上浮穴郡久万高原町七鳥 lng:  132.9862943 lat:  33.6578348
浄瑠璃寺　愛媛県 ,  愛媛県 松山市浄瑠璃町 lng:  132.815242 lat:  33.749282
八坂寺　愛媛県 ,  愛媛県 松山市浄瑠璃町 lng:  132.815242 lat:  33.749282
西林寺　愛媛県 ,  愛媛県 松山市高井町 lng:  132.816453 lat:  33.796264
浄土寺　愛媛県 ,  愛媛県 松山市鷹子町 lng:  132.9927153009765 lat:  34.0468397
繁多寺　愛媛県 ,  愛媛県 松山市畑寺町 lng:  132.42194126001942 lat:  33.46302005
石手寺　愛媛県 ,  愛媛県 松山市石手－－ lng:  132.7963926 lat:  33.8467695
太山寺　愛媛県 ,  愛媛県 松山市太山寺町 lng:  132.7203772 lat:  33.8839201
円明寺　愛媛県 ,  愛媛県 松山市和気 lng:  132.739489 lat:  33.8911074
延命寺　愛媛県 ,  愛媛県 今治市阿方甲 lng:  132.9639973 lat:  34.0668412
南光坊　愛媛県 ,  愛媛県 今治市別宮町－ lng:  132.9953058 lat:  34.0682265
泰山寺　愛媛県 ,  愛媛県 今治市小泉－－ lng:  132.7203772 lat:  33.8839201
栄福寺　愛媛県 ,  愛媛県 今治市玉川町八幡 lng:  132.9782314 lat:  34.0297786
仙遊寺　愛媛県 ,  愛媛県 今治市玉川町別所 lng:  132.9773496 lat:  34.0132837
国分寺　愛媛県 ,  愛媛県 今治市国分－－ lng:  133.0242384 lat:  34.0252747
横峰寺　愛媛県 ,  愛媛県 西条市小松町石鎚 lng:  133.1110909 lat:  33.8376181
香園寺　愛媛県 ,  愛媛県 西条市小松町南川 lng:  133.1034109 lat:  33.893603
宝寿寺　愛媛県 ,  愛媛県 西条市小松町新屋敷 lng:  133.1151159 lat:  33.8973011
吉祥寺　愛媛県 ,  愛媛県 西条市氷見乙 lng:  133.1291511 lat:  33.8957732
前神寺　愛媛県 ,  愛媛県 西条市州之内甲 lng:  133.1605615 lat:  33.8903607
三角寺　愛媛県 ,  愛媛県 四国中央市金田町三角寺 lng:  133.5866363 lat:  33.9673725
雲辺寺　徳島県 ,  徳島県 三好市池田町白地－ lng:  133.7229222 lat:  34.0351293
大興寺　香川県 ,  香川県 三豊市山本町辻 lng:  133.71914297374104 lat:  34.1022869
神恵院　香川県 ,  香川県 観音寺市八幡町－－ lng:  133.64730950356878 lat:  34.13398665
観音寺　香川県 ,  香川県 観音寺市八幡町－－ lng:  133.6554569 lat:  34.124273
本山寺　香川県 ,  香川県 三豊市豊中町本山甲 lng:  133.6937985133103 lat:  34.139606549999996
弥谷寺　香川県 ,  香川県 三豊市三野町大見乙 lng:  133.7241876 lat:  34.230281
曼荼羅寺　香川県 ,  香川県 善通寺市吉原町－ lng:  133.7503177 lat:  34.2233323
出釈迦寺　香川県 ,  香川県 善通寺市吉原町 lng:  133.7502513 lat:  34.2191392
甲山寺　香川県 ,  香川県 善通寺市弘田町－ lng:  133.7657237 lat:  34.2330662
善通寺　香川県 ,  香川県 善通寺市善通寺町－－ lng:  133.78922 lat:  34.2299924
金倉寺　香川県 ,  香川県 善通寺市金蔵寺町 lng:  133.7809675 lat:  34.2495723
道隆寺　香川県 ,  香川県 仲多度郡多度津町北鴨－－ lng:  133.76289294017013 lat:  34.27669405
郷照寺　香川県 ,  香川県 綾歌郡宇多津町 lng:  133.8252787 lat:  34.3107957
高照院　香川県 ,  香川県 坂出市西庄町－ lng:  133.8586386 lat:  34.3186022
国分寺　香川県 ,  香川県 高松市国分寺町国分 lng:  133.9441627 lat:  34.3030901
白峯寺　香川県 ,  香川県 坂出市青海町 lng:  133.9270371 lat:  34.3344627
根来寺　香川県 ,  香川県 高松市中山町 lng:  133.966296 lat:  34.3443761
一宮寺　香川県 ,  香川県 高松市一宮町 lng:  134.0265283971522 lat:  34.2863763
屋島寺　香川県 ,  香川県 高松市屋島東町 lng:  134.10106303395688 lat:  34.3574992
八栗寺　香川県 ,  香川県 高松市牟礼町牟礼 lng:  134.133009 lat:  34.3493325
志度寺　香川県 ,  香川県 さぬき市志度 lng:  134.1796409 lat:  34.3242502
長尾寺　香川県 ,  香川県 さぬき市長尾西 lng:  134.1717088 lat:  34.2666753
大窪寺　香川県 ,  香川県 さぬき市多和兼割 lng:  134.2068728 lat:  34.1913639

　数値だけ見てもわかりにくいので図化して確認します。

　グラフを横並びにして、適度な間隔を開けて見やすくしました。やり方は全部ChatGPT先生にお聞きしました。

#buffer()メソッドを使って各霊場を中心にした5kmの範囲を定義
buffer_distance = 10_000 #1km=1,000m
gdf_88temples_UTM55N_buffer = gdf_88temples_UTM55N.buffer(buffer_distance)

#2つのグラフを横に並べるfigureを作成
fig, axs = plt.subplots(1, 2, figsize=(8,4))

#左側のグラフを作成
gdf_88temples_UTM55N_buffer.plot(ax=axs[0], edgecolor='black')
axs[0].set_title('Buffered Data after CSR Transform')
gdf_88temples.plot(ax=axs[1], edgecolor='black')
axs[1].set_title('Original Data before CSR Transform')

plt.subplots_adjust(wspace=0.5)

　全ての点が一定の場所にまとまったように見えます。経度も132～135程度に収まっていますので、まずはこれでデータができたとします。

　次に、四国の地図上に四国巡礼88霊場の位置をプロットしてみようと思います。

　四国の形状は四国四県それぞれのポリゴンデータを合成した国土交通省の国土数値情報（https://nlftp.mlit.go.jp/ksj/）の行政区域のデータを活用します。国土数値情報の行政区域のページに行き、徳島県、香川県、愛媛県、高知県のデータをShapeFile形式でダウンロードして、Zipを解凍したうえで、GoogleDriveの任意のフォルダーにアップロードしておきます。

　次に、四県それぞれのShapeFileを呼び出して、その１行目(Index:0)を取り出し、新たに定義した'shikoku_gdf'というGeoDataFrameに追加していきます。また、四国地図は今後使えるかもしれないので、'all_shikoku.shp’という名称で保管しておきます。コードは以下の通り。

#四国地方全体のshapefileを作成するコード

# 四国地方のpolygonデータを格納するGeoDataFrame'shikoku_gdf'を生成
shikoku_gdf = gpd.GeoDataFrame(columns = ["N03_001", "N03_002", "N03_003", "N03_004", "N03_007", "geometry"])

# Read the shapefile into a geopandas GeoDataFrame from 徳島県
temp_gdf = gpd.read_file('/content/drive/MyDrive/Colab Notebooks/202303_Shikoku88Data/20230327_四国四県GIS情報/N03-20220101_36_GML/N03-22_36_220101.shp')
shikoku_gdf.loc[len(shikoku_gdf)] = temp_gdf.loc[0]

# Read the shapefile into a geopandas GeoDataFrame from 香川県
temp_gdf = gpd.read_file('/content/drive/MyDrive/Colab Notebooks/202303_Shikoku88Data/20230327_四国四県GIS情報/N03-20220101_37_GML/N03-22_37_220101.shp')
shikoku_gdf.loc[len(shikoku_gdf)] = temp_gdf.loc[0]

# Read the shapefile into a geopandas GeoDataFrame from 愛媛県
temp_gdf = gpd.read_file('/content/drive/MyDrive/Colab Notebooks/202303_Shikoku88Data/20230327_四国四県GIS情報/N03-20220101_38_GML/N03-22_38_220101.shp')
shikoku_gdf.loc[len(shikoku_gdf)] = temp_gdf.loc[0]

# Read the shapefile into a geopandas GeoDataFrame from 高知県
temp_gdf = gpd.read_file('/content/drive/MyDrive/Colab Notebooks/202303_Shikoku88Data/20230327_四国四県GIS情報/N03-20220101_39_GML/N03-22_39_220101.shp')
shikoku_gdf.loc[len(shikoku_gdf)] = temp_gdf.loc[0]

shikoku_gdf.to_file('/content/drive/MyDrive/Colab Notebooks/202303_Shikoku88Data/all_shikoku.shp', encoding='utf-8')

shikoku_gdf.head()

　'shokoku_gdf'を出力すると以下のようになります。

	N03_001	N03_002	N03_003	N03_004	N03_007	geometry
0	徳島県	NaN	NaN	NaN	NaN	MULTIPOLYGON (((134.31889 33.53861, 134.31889 ...
1	香川県	NaN	NaN	NaN	NaN	MULTIPOLYGON (((133.63603 34.10358, 133.63601 ...
2	愛媛県	NaN	NaN	NaN	NaN	MULTIPOLYGON (((132.45428 32.88539, 132.45406 ...
3	高知県	NaN	NaN	NaN	NaN	MULTIPOLYGON (((132.54156 32.70275, 132.54152 ...

　少しまどろっこしいですが、せっかく'all_shikoku.shp'を作成したので、それを読み込んで、図化します。

shikoku_gdf = gpd.read_file('/content/drive/MyDrive/Colab Notebooks/202303_Shikoku88Data/all_shikoku.shp')

shikoku_gdf.plot()

　そうすると（１）に掲載した四国地図が図化できます。ちなみに、横軸・縦軸は経度・緯度になっていることがわかります。

　次は四国地図に四国88箇所霊場をプロットします。

　四国地図も四国88箇所霊場の位置情報も同じJGD2011世界測地系の経度緯度で設定されているので、そのままプロットしてOKです。

#四国地図上に四国88箇所霊場の位置をプロット
fig, ax = plt.subplots(figsize=(10,10))

#四国地図をプロット
shikoku_gdf.plot(ax=ax)

#四国88箇所霊場の位置をプロット
gdf_88temples.plot(ax=ax, edgecolor='black', color='red', alpha=1.0, label='四国88箇所霊場')

ax.legend()

　なぜか瀬戸内海の海上に２点プロットされています。よくわからないので、各プロットの横に霊場の名前を書きます。matplotlibで日本語を表示するためには、

　import japanize_matplotlib

があらかじめ必要になりますので、お忘れ無きように。。。

字が小さくて見づらいかもしれませんが、瀬戸内海の海上に存在するプロットは切幡寺と藤井寺でした。この２つの霊場については、アナログチックにGeoDataFrameの緯度経度を直接的にGoogleMapで計測した経度緯度に修正します。

GeoDataFrameの特定の行に入力されている位置情報を修正する方法をChatGPTさまにお聞きしてみました。ChatGPTさまは「もう少しSimpleな方法はないですか？」と聞くと、ほぼ毎回さらにシンプルな方法を提案してくれます。

Cityの名称を指定して、geometryを修正する方法を教えて

PythonのGeopandasにおいて、Cityの名称を指定して、その行のgeometryを修正するには、まず"loc"メソッドを使用して、Cityの名称が一致する行を取得し、その行のインデックスを取得します。次に、取得したインデックスを使用して、"at"メソッドを使用して"geometry"列に新しいgeometryオブジェクトを代入する方法があります。以下に例を示します。

import geopandas as gpd
from shapely.geometry import Point

# Create a GeoDataFrame with two points
df = gpd.GeoDataFrame({'City': ['New York', 'Los Angeles'],
                       'Coordinates': [Point(-74.0060, 40.7128), Point(-118.2437, 34.0522)]})

# Print the original geometry of the second point
print(df.at[1, 'Coordinates'])

# Modify the geometry of the second point by City name
city_name = 'Los Angeles'
new_geometry = Point(-118.2437, 37.7749)
df.at[df.loc[df['City'] == city_name].index[0], 'Coordinates'] = new_geometry

# Print the new geometry of the second point
print(df.at[df.loc[df['City'] == city_name].index[0], 'Coordinates'])

この例では、2つのPOINTからなるGeoDataFrameを作成し、"Los Angeles"というCity名を持つ行のgeometryを修正します。まず、2つ目の点の元のgeometryを印刷し、その後、"loc"メソッドを使用して、"City"列で"Los Angeles"と一致する行を取得し、その行のインデックスを取得します。次に、取得したインデックスを使用して、"at"メソッドを使用して、"Coordinates"列に新しいPointオブジェクトを設定します。最後に、修正されたgeometryを印刷して、変更が正しく適用されたことを確認します。

もう少しシンプルな方法はないですか

はい、よりシンプルな方法があります。Pandasのデータフレームと同様に、GeopandasのGeoDataFrameでも、条件を指定して行を選択し、列を指定して値を変更することができます。

以下は、City名が"Los Angeles"の行のgeometryを修正する例です。

import geopandas as gpd
from shapely.geometry import Point

# Create a GeoDataFrame with two points
df = gpd.GeoDataFrame({'City': ['New York', 'Los Angeles'],
                       'Coordinates': [Point(-74.0060, 40.7128), Point(-118.2437, 34.0522)]})

# Print the original geometry of the second point
print(df.loc[df['City'] == 'Los Angeles', 'Coordinates'])

# Modify the geometry of the second point by City name
df.loc[df['City'] == 'Los Angeles', 'Coordinates'] = Point(-118.2437, 37.7749)

# Print the new geometry of the second point
print(df.loc[df['City'] == 'Los Angeles', 'Coordinates'])

この例では、"Los Angeles"のCity名を持つ行のgeometryを修正します。まず、条件式df['City'] == 'Los Angeles'を使用して、City名が"Los Angeles"の行を選択します。そして、"Coordinates"列に新しいPointオブジェクトを設定して、geometryを修正します。最後に、修正されたgeometryを印刷して、変更が正しく適用されたことを確認します。

この方法は、条件式による行の選択と、列指定による値の変更を一度に行うことができるため、より簡潔かつわかりやすい方法です。

ちょっと手間がかかりますが、素晴らしい回答です。藤井寺と切幡寺の経度緯度を力業で修正します。

#藤井寺と切幡寺の位置情報をGoogleMapから取得した緯度経度に修正
#藤井寺...34.05152663252319, 134.34871241761365
#切幡寺...34.1075959887686, 134.3040334198125

gdf_88temples.loc[gdf_88temples['Name'] == '藤井寺　徳島県', 'geometry'] = Point(134.34871241761365, 34.05152663252319)
gdf_88temples.loc[gdf_88temples['Name'] == '切幡寺　徳島県', 'geometry'] = Point(134.3040334198125, 34.1075959887686)

再度、図化してみます。コードは先ほど掲載しているものと同じです。

#四国地図上に四国88箇所霊場の位置をプロット
fig, ax = plt.subplots(figsize=(10,10))

#四国地図をプロット
shikoku_gdf.plot(ax=ax)

#四国88箇所霊場の位置をプロット
gdf_88temples.plot(ax=ax, edgecolor='black', color='red', alpha=1.0, label='四国88箇所霊場')

# 各点にラベルを表示
for i, point in gdf_88temples.iterrows():
    ax.text(point.geometry.x+0.02, point.geometry.y, point['Name'], fontsize=8)

ax.legend()

だいたい理想とするところに近づいたような気がします。改めて見ると、四国８８か所霊場は四国全体に比較的均等に分散している訳ではないということがよくわかりました。徳島県の西部、高知県北部あたりは本当に霊場がありません。

　これでとりあえず満足して、それぞれの霊場について半径10kmの範囲を定義したいと思います。距離や面積を扱う場合には経度緯度表示（JGD2011測地系座標、EPSG:6668）は不適なので、直交座標系（UTM zone 55N EPSG:6670）に変換します。

#四国巡礼88箇所霊場のGeoDataFrameのReferenceCoordinateSystemを
#JGD2011測地系座標(EPSG:6668)をUTM zone 55N(EPSG:6670…日本の東経132度から135度に適用)
#に変換
gdf_88temples_UTM55N = gdf_88temples.to_crs(epsg=6670)

gdf_88temples_UTM55N.head()

#buffer()メソッドを使って各霊場を中心にした5kmの範囲を定義
buffer_distance = 10_000 #1km=1,000m
gdf_88temples_UTM55N_buffer = gdf_88temples_UTM55N.buffer(buffer_distance)

#2つのグラフを横に並べるfigureを作成
fig, axs = plt.subplots(1, 2, figsize=(8,4))

#左側のグラフを作成
gdf_88temples_UTM55N_buffer.plot(ax=axs[0], edgecolor='black')
axs[0].set_title('Buffered Data after CSR Transform')
gdf_88temples.plot(ax=axs[1], edgecolor='black')
axs[1].set_title('Original Data before CSR Transform')

plt.subplots_adjust(wspace=0.5)

　図化するとこのようになりました。

　まあまあ四国の形通りになっています。四国地図の座標系はJGD2011測地系座標、EPSG:6668なので、Buffer()メソッドで定義した半径10㎞の図形をJGD2011測地系座標に再度変換します。

#buffer()したGeoDataFrameno座標系をJGD2011測地系座標(EPSG:6668)に戻す
gdf_88temples_EPSG6668_buffer = gdf_88temples_UTM55N_buffer.to_crs(epsg=6668)

gdf_88temples_EPSG6668_buffer.head()

0    POLYGON ((134.62561 34.18100, 134.62471 34.172...
1    POLYGON ((134.60567 34.16382, 134.60477 34.155...
2    POLYGON ((134.59844 34.14984, 134.59755 34.141...
3    POLYGON ((134.57068 34.03501, 134.56979 34.026...
4    POLYGON ((134.69338 34.07283, 134.69247 34.064...
dtype: geometry

　出力は経度緯度表示になっているので、よろしいようです。

　では、これを四国の地図と合成して表示します。

#四国地図上に四国88箇所霊場の位置をプロット
fig, ax = plt.subplots(figsize=(10,10))

#四国地図をプロット
shikoku_gdf.plot(ax=ax)

#四国88箇所霊場の位置をプロット
gdf_88temples.plot(ax=ax, edgecolor='black', color='red', alpha=1.0, label='四国88箇所霊場')

# 各点にラベルを表示
for i, point in gdf_88temples.iterrows():
    ax.text(point.geometry.x+0.02, point.geometry.y, point['Name'], fontsize=8)

#Buffer()したデータを表示
gdf_88temples_EPSG6668_buffer.plot(ax=ax, edgecolor='pink', linestyle='solid', linewidth=1, color='gray', alpha=0.2)

#判例を表示
#ax.legend(['四国88箇所霊場'])
ax.legend()

　出力です。

　最後に苦労して作成した四国巡礼88箇所霊場の経度緯度情報をshapeファイル形式で保存しておくことにします。

gdf_88temples.to_file('/content/drive/MyDrive/Colab Notebooks/202303_Shikoku88Data/shikoku_88temples.shp', encoding='utf-8')

gdf_88temples.head()

　引き続き四国巡礼88箇所霊場の位置でGeoPandasの世界を楽しんでいきたいと思います。

	Name	Address	geometry
0	霊山寺　徳島県	徳島県鳴門市大麻町　板東	POINT (134.51733 34.18415)
1	極楽寺　徳島県	徳島県鳴門市大麻町桧	POINT (134.49741 34.16696)
2	金泉寺　徳島県	徳島県板野郡板野町大寺字亀山下	POINT (134.49020 34.15297)
3	大日寺　徳島県	徳島県板野郡板野町黒谷字居内	POINT (134.46259 34.03811)
4	地蔵寺　徳島県	徳島県板野郡板野町羅漢字林東	POINT (134.58525 34.07603)

このブログを検索

Shikuu Kitashirakawa - Kitchen garden, Running/Marathon, Python etc.

Python備忘録　OpenStreetMapとgeocoder、さらにChatGPTの活用（２）

コメント

コメントを投稿

このブログの人気の投稿

【完結】ランニング、お食事　２０２２年５月～２０２２年１２月

ランニング、グルメ、ドライブ　２０２３年４月〜

ランニング、グルメ、クライム　２０２３年７月〜

Python備忘録 OpenStreetMapとgeocoder、さらにChatGPTの活用（２）

コメント

コメントを投稿

このブログの人気の投稿

【完結】ランニング、お食事 ２０２２年５月～２０２２年１２月

ランニング、グルメ、ドライブ ２０２３年４月〜

ランニング、グルメ、クライム ２０２３年７月〜

Python備忘録　OpenStreetMapとgeocoder、さらにChatGPTの活用（２）

【完結】ランニング、お食事　２０２２年５月～２０２２年１２月

ランニング、グルメ、ドライブ　２０２３年４月〜

ランニング、グルメ、クライム　２０２３年７月〜