🛤️

駅から一番遠い地点を探る

2023/01/10に公開

Python

ネタ

日本の陸地で、鉄道駅から一番遠い地点はどこか？

普通に考えれば南鳥島^[1]になってしまいますが、これを鉄道がある主要5島（本州・北海道・九州・四国・沖縄本島）の中において考えてみます。^[2]

道のりで判定できると最高ですが、難しいので^[3]直線距離とします。厳密な解ではなく、ある程度この辺、と分かれば満足です。

地図情報もPython数値計算も不得手なので、諸々うまい実装ではないと思います。

成果物

Streamlit

東京23区の例です。駅から遠いほど赤・白になっています。江東区若洲海浜公園 (東京ゲートブリッジ付近) が駅から一番遠い結果です。

距離変換画像

こちらはボロノイ図で、それぞれの駅の最寄り領域がわかります。若洲海浜公園の最寄りは緑の米印で示されていて、新木場駅です (約3.8km)。

ボロノイ図

実装

手法の流れ

全部の駅の座標（緯度・経度）を入手
日本の輪郭データを入手
遠い地点の判定には、島の輪郭内における、駅を基点とする距離変換
どの駅が最寄りなのかの判定には、駅を母点とするボロノイ図

環境

Python 3.10

requirements.txt

opencv-python-headless==4.7.0.68
streamlit==1.16.0
matplotlib==3.6.2
pyproj==3.4.1
shapely==2.0.0
more-itertools==9.0.0

データ参照元

鉄道駅LOD: https://uedayou.net/jrslod/
Linked Open Addresses Japan: https://uedayou.net/loa/
地図蔵 - 47都道府県のポリゴンデータ geojson: https://japonyol.net/editor/article/47-prefectures-geojson.html

手順

駅の緯度経度の取得

以下の手順で、ひたすら全国の路線を選んではダウンロードしていきました。

今回も、GeoJSON特化のライブラリは使わずにゴリゴリ書きました。あまり効率は良くなさそうなコードですが、以前の記事同様に駅の位置を集めていきます。

import json

def load_train_geojson(path: str) -> pd.DataFrame:
    with open(path, encoding="utf-8-sig") as file:
        json_obj = json.load(file)
        rows = [
            {
                "name": f["properties"]["name"],
                "lon": f["geometry"]["coordinates"][0],
                "lat": f["geometry"]["coordinates"][1],
            }
            for f in json_obj["features"]
            if f["geometry"]["type"] == "Point"
        ]
        return pd.DataFrame(data=rows)


df = load_train_geojson("geojson/近畿日本鉄道近鉄京都線.geojson")
print(df)

     name        lon       lat
0      京都  135.75692  34.98459
1      東寺  135.75259  34.97968
2      十条  135.75247  34.97443
3    上鳥羽口  135.75254  34.96540
...
24     平城  135.78396  34.70117
25  大和西大寺  135.78401  34.69315

この調子で全駅の緯度経度を集めます。欲しいデータとしては各路線ごとにデータが分かれていなくてよいので、pandas.concatで全部連結してしまいます。

from pathlib import Path

paths = Path("geojson").glob("*.json")
lonlat = pd.concat(load_train_geojson(str(p)) for p in paths)

もし「四国だけ」「関東だけ」といった地域限定の処理を行いたい場合は、この時点で使うgeojsonの絞り込みをしておくことになります。路線名だけで絞り込めることが多いですが、特にJRの長大な路線だと限界があります（例: 「中央本線」は東京都から愛知県まであり、関東だけ対象にしたいとしても山梨県以西が入ってしまいます）。駅単位の絞り込みを要するでしょう。本記事では詳細は省きますが、私の実装では路線単位の絞り込みと緯度経度による絞り込み（都道府県境の輪郭の内側かどうか）を併用しました [路線名絞り込みの参考] [緯度経度絞り込みの参考]。

駅の位置を平面直角座標系に変換

今回も最終的に画像処理的なアプローチをするので、緯度経度を平面直角座標系に変換します。

説明は以前の記事に譲ります。

前回の記事で示したget_transformerメソッドを利用してpyproj.Transformerオブジェクトを得ます。そしてDataFrameの各緯度経度を書き換えていきます。

transformer = get_transformer("京都府", "")

def map_row(row: pd.Series):
    lat, lon = transformer.transform(row["lat"], row["lon"])
    return {"name": row["name"], "lon": lon, "lat": lat}
	
utm = lonlat.apply(map_row, axis=1, result_type="expand")

     name           lon            lat
0      京都 -22192.291199 -112621.103493
1      東寺 -22588.953196 -113164.797225
2      十条 -22601.351727 -113747.149542
3    上鳥羽口 -22597.439803 -114748.856836
...
24     平城 -19791.364085 -144065.671865
25  大和西大寺 -19788.692602 -144955.293993

島・都道府県の輪郭の緯度経度を得る

再掲ですが以下2つのサイトからgeojsonを頂きました。都道府県単位の処理であれば後者のみで賄えます。市町村単位などもっと細かい処理をしたければ前者が使えます。

Linked Open Addresses Japan: https://uedayou.net/loa/
地図蔵 - 47都道府県のポリゴンデータ geojson: https://japonyol.net/editor/article/47-prefectures-geojson.html

駅座標の工程と異なるのは、欲しいデータがgeojsonから直接は得られない点です^[4]。例えば九州について処理したい場合、「九州島」の輪郭データは無いので、以下工程が必要になります。

福岡・佐賀・長崎・熊本・大分・宮崎・鹿児島の7県の輪郭を結合する
九州島以外の離島（対馬、屋久島等）を除外する [任意]

輪郭の結合のため、shapelyを利用しました。

以下は九州の輪郭を得る例です。まず九州7県の輪郭をshapeオブジェクトとして集めるコードです。

import shapely

with open("geojson/prefectures.geojson", encoding="utf-8-sig") as file:
    gj = json.load(file)
    
shapes = tuple(
    shapely.geometry.shape(f["geometry"])
    for f in gj["features"]
    if f["properties"]["name"] in {"福岡県", "佐賀県", "長崎県", "熊本県", "大分県", "宮崎県", "鹿児島県"}
)

shapeの種類は .geom_type で調べることができ、MultiPolygonまたはPolygonになっています。多角形同士の連結はshapely.unary_unionメソッドでできます。

merged_shape = shapely.unary_union(polygons)

unionしても離島等はくっつかないため、多角形はまだ1つになりません。もし1つに絞りたければ、最大のものだけ残すといった処理を書きます。

largest_contour = max(merged_shape.geoms, key=lambda poly: poly.area)

最後に輪郭データも駅同様にpyproj.Transformerで平面直角座標系に変換しておきます。

transformer = get_transformer("熊本県", "")
area_contour = np.array([transformer.transform(lat, lon)[::-1] for lon, lat in largest_contour])

座標の正規化

駅や輪郭の座標が (-22597.439803 -114748.856836) のような画像処理的に扱いづらい状態なので、正規化します。前回の記事同様です。

駅と輪郭とで同じ倍率で縮めたいので、いったん駅座標行列と輪郭座標行列（双方Nx2）を縦にくっつけてから正規化し、また切り離すということをしました。（これもセオリーとは逸れてる気がする...）

# 緯度経度座標列(UTM座標)から平面画像上の座標に正規化し変換
def normalize(
    station_df: pd.DataFrame,
    island_contours: tuple[npt.NDArray[np.float64], ...],
    image_size: int
) -> tuple[pd.DataFrame, tuple[npt.NDArray, ...]]:
    station_coordinates = station_df[["lon", "lat"]].values
    combi = np.vstack((station_coordinates, *island_contours))

    x, y = combi.T
    xmin, xmax, _, _ = cv2.minMaxLoc(x)
    ymin, ymax, _, _ = cv2.minMaxLoc(y)
    if (xmax - xmin) > (ymax - ymin):
        scale = xmax - xmin
    else:
        scale = ymax - ymin

    normx = ((x - xmin) / scale) * image_size
    normy = ((y - ymin) / scale) * image_size
    combi_norm = np.vstack((normx, normy)).T

    ret_station_coordinates, all_island_contours = np.vsplit(combi_norm, [len(station_df)])
    ret_station_df = station_df.copy()
    ret_station_df["lon"] = ret_station_coordinates[:, 0]
    ret_station_df["lat"] = ret_station_coordinates[:, 1]

    sections = list(itertools.accumulate(len(c) for c in island_contours))[:-1]
    ret_island_contours = np.vsplit(all_island_contours, sections)

    return ret_station_df, ret_island_contours

結果が問題ないかどうか描画してみる一例として、matplotlibを使う方法を置いておきます。st.pyplotでプロットを埋め込めるので、Streamlit開発の最中でも使いやすいです。注意として、アスペクト比が自動調整されるので地図としては歪みます。必要なら ax.set_aspect("equal") のようにして調整します。

import matplotlib.pyplot as plt

fig, ax = plt.subplots()
for data in (station_df["lon", "lat"].values, island_contours):
    for x, y in data: 
        ax.scatter(x, y, s=1)
ax.grid()
st.pyplot(fig)

天草や甑島などが消えているのは、最大contourだけに絞った影響です。

距離変換

座標が出そろったので、駅からの距離の計算に入っていきます。

今回は、そもそも座標を正規化している時点で当然ですが厳密解は求めていないので、簡単さを考えOpenCVのcv2.distanceTransform（距離変換）を使うことにしました。距離変換については以下などをご参照ください。

つまり、黒一色の画像に対し、駅の位置に白い点を打っていき、距離変換をすれば、駅から遠い箇所が現れるという作戦です^[5]。

距離変換ののち、対象の陸地部分だけをマスク処理で切り抜くようにします。

# 距離変換の結果を得る
def run_distance_transform(
    station_df: pd.DataFrame,
    island_contours: tuple[npt.NDArray[np.float64], ...],
    image_size: int,
) -> tuple[npt.NDArray[np.float64], FarthestPoint]:
    img = np.full((image_size, image_size, 1), 255, dtype=np.uint8)
    for p in station_df[["lon", "lat"]].values:
        cv2.rectangle(img, p.astype(int), (p + 1).astype(int), (0, 0, 0), thickness=2)
    dist = cv2.distanceTransform(img, cv2.DIST_L2, 5)

    mask = np.zeros_like(dist, dtype=np.uint8)
    for one_island in island_contours:
        cv2.fillPoly(
            mask,
            [np.array([f.astype(int) for f in one_island])],
            (255, 255, 255),
        )
    dist = dist * (mask / 255)
    cv2.normalize(dist, dist, 0, 1.0, cv2.NORM_MINMAX)

ラスタ画像はy軸が下向きなので、南北が逆転した結果になります。表示直前に反転させればよいです。

dist = cv2.flip(dist, 0)
st.image(dist, caption="Distance Transform")

距離変換画像 (グレースケール)

ここでdistは0.0～1.0のfloat画像で、自動的に0～255のグレースケールに変換して表示してくれています。いまいちグレースケールだと見た目がはっきりしないので、疑似カラーを使います。cv2.applyColorMapを使います。

様々な色合いがフラグ1つで選べるので、お好みのものを選択します。以下など参考になります。

dist_u8 = cv2.applyColorMap((dist * 255).astype(np.uint8), cv2.COLORMAP_HOT)
dist_u8 = cv2.copyMakeBorder(dist_u8, 50, 50, 50, 50, cv2.BORDER_CONSTANT, (0, 0, 0))

dist_u8 = cv2.flip(dist_u8, 0)
st.image(dist_u8, channels="BGR", caption="Distance Transform")

距離変換画像 (疑似カラー)

この例では明らかに襟裳岬が一番白いので、最も駅から遠いと分かりました。計算でそれを求めるには、distの値が最も大きいインデックスを求めます。上の図の×印はこのmax_locの位置を示しています。

_, _, _, max_loc = cv2.minMaxLoc(dist)

ボロノイ分割

上の結果から、例えば北海道では襟裳岬が最も遠いと分かりました。ではどの駅が最寄りなのでしょうか？

愚直に計算して十分間に合いますが、見た目が面白そうなのでボロノイ図を作ってみることにしました。

駅の位置を母点としてボロノイ図を作ります。OpenCVではcv2.Subdiv2Dが担当します。

def voronoi(station_df: pd.DataFrame, image_size: int):
    points = station_df[["lon", "lat"]].values
    subdiv = cv2.Subdiv2D((0, 0, image_size + 1, image_size + 1))
    for p in points:
        subdiv.insert((p[0], p[1]))
    facets, centers = subdiv.getVoronoiFacetList([])
        
    img = np.zeros((image_size + 100, image_size + 100, 3), np.uint8)
    for p in centers:
        cv2.drawMarker(img, (p + 50).astype(int), (0, 0, 255), thickness=3)
    cv2.polylines(img, [(f + 50).astype(int) for f in facets], True, (0, 255, 255), thickness=2)

輪郭の描画コードは省略しますが、こんな感じになります。

緑の※印が襟裳岬に直線距離で最短の駅で、新吉野駅です。計算で求めるには、（もっと良い方法があるかもしれませんが、）各ボロノイ領域 facets についてcv2.pointPolygonTestを使って襟裳岬座標が含まれるかをテストしていくと分かります。

import more_itertools

# 距離変換画像への cv2.minMaxLoc で求められた、画像内での距離最大の座標
farthest = (1162, 1770)

_, nearest_station_location = more_itertools.first_true(
    zip(facets, centers), pred=lambda c: cv2.pointPolygonTest(c[0], farthest, measureDist=False) >= 0
)

最寄りの駅との距離

以上で、最も駅から遠い地点と、そこの最寄り駅がわかったとして、どのくらい離れているか直線距離を求めてみます。

こちらなど参考になります。

メートル単位で求められます。この例では約100kmですね。

station_pos = (143.24036, 41.92318)  # 新吉野駅
farthest_pos = (143.609, 42.778)  # 襟裳岬

g = Geod(ellps="WGS84")
_, _, distance_2d = g.inv(*station_pos, *farthest_pos)
print(distance_2d)  # 99692.91134424254

Google Mapsの距離測定ツールで検算してみました。誤差300mくらいで、素人がやったにしてはまずまずの結果です。

結果を見る

Streamlit再掲: https://shimat-station-voronoi-main-czmu73.streamlit.app/

全国(主要五島)

先に本命の結果から。一番駅から遠いのは、例示していた 襟裳岬 でした。前述の通り、最寄りは新吉野駅で約100kmあります。

襟裳岬が首位に躍り出たのは、日高本線が2021年4月で大部分廃止になって以降で、それまでは最寄りの様似駅まで約34kmでした。

以下では各島・地域ごとに見ていきます。

北海道

北海道の結果はご承知の通りなので、ここでは日高本線(鵡川～様似)の廃止前はどうだったのかを見てみました。

今度は知床岬が首位になりました。最寄りの知床斜里駅まで約72km。ちなみにこうなった場合、全国首位は別の場所になります。

到達困難具合を含めれば、今も昔も知床岬が圧倒的であり続けていますね。

本州

本州島でのトップは、能登半島の先端（珠洲岬付近?）のようです。最寄りは穴水駅で約50km。

ここも廃線があり、のと鉄道能登線が2005年に廃止になるまでは能登半島の先端まで鉄道がありました。

次点が尾瀬の付近、その次が紀伊半島の真ん中（奈良県十津川村あたり）でしょうか。これらには廃線は無いようです。島根県付近の赤さも目立ちますね（三江線が2018年で廃止）。

四国

四国のトップは佐田岬で、最寄り八幡浜駅から約41km。まあ、すごい尖りっぷりだしそうだろうなあという順当な結果でした。

四国は割と廃線が少ないようで、その面からも順位の変動は昔からほぼなさそうです。

ただし佐田岬については、豊後水道を挟んで対岸、大分県の幸崎駅までが23km程度だったりします。それを加味すると、高知・愛媛または高知・徳島の県境一帯のほうが距離が上回ります。

九州

九州は、主要五島で唯一、内陸が一番遠い結果です。(32.51207, 131.25962) の付近で、住所は「宮崎県東臼杵郡椎葉村松尾」です。Google Mapsによると付近のスポットとしては「椎葉村大いちょう展望台」「仙人の棚田」などがあります。

最寄りは、処理結果としては太平洋側のJR日向市駅で約35.9kmと出ました。ただし熊本県側の内陸に延びる鉄道とも等距離のようで、くま川鉄道の湯前駅、南阿蘇鉄道の高森駅も同じく約36kmです。

九州は結構廃線が多いですが、順位を塗り替えそうな廃線はなさそうです (参考)。次点で大隅半島の白さも目立ち、これは国鉄大隅線の廃線によるところです。

沖縄本島

とても分かりやすい図になりました。北の先端付近、奥港あたりが最も遠く、最寄りはゆいレールのてだこ浦西駅で約87km。なんと全国2位の遠さです。沖縄は狭いと見くびってはいけませんね。

戦前は嘉手納まで鉄道があったそうで、その場合は10kmほど近くなります。

東京23区

本州のくくりが広かったので絞ってみました。はじめの方でも例示しました東京23区です。

遠かったのは江東区若洲海浜公園 (東京ゲートブリッジ付近)で、最寄りは新木場駅です (約3.8km)。

島全体の図とは違った趣があります。くまなく駅が敷き詰められていて、いかに鉄道会社または行政が交通の穴を埋めに行っているかを感じられました。特に日暮里・舎人ライナーのボロノイ領域がとても均等なのが個人的にツボですね。

私の手抜きの問題で領域外(23区外や埼玉県など)の駅を含めて処理できていないので、白く見えても実は駅に近いことがあります。例えば左上の練馬区大泉中央公園付近が白いですが、北の埼玉県側の和光市駅の方が近いので、本当はもっと黒くなるはずです。

大阪市

最後に大阪市です。こちらも駅がくまなくありそうなので。

沖合の埋立地、夢洲（ゆめしま）の先っぽが一番遠く、最寄りコスモスクエア駅から約3.7kmと出ました。

2025年の万博会場となり、駅も設置予定です。まもなく首位から降りるはずです。そうなると以後は西淀川区の中島PA付近（淀川河口の北）か、大正区か、というところですね。

地形データの問題なのか仕様なのか、水路があまり反映されておらず、陸地が広いように見えています。川の扱いか海の扱いかという差なのかもしれません。

大阪もまた、くまなく駅を敷き詰めに行っている感じがとてもよくわかります。今里筋線とか、地下鉄のルートがとても狙いに行っているように見えてきます。

ということで、技術的に目新しくはないですがまた昔から知りたかったことが満たされて満足しました。

脚注

(本州の)東京から約1900km ↩︎
東京のお台場や関西空港駅等、主要5島に近接して駅がある島を含みます。 ↩︎
技術的にはできそうですが、主に金銭的に。 ↩︎
別のオープンデータでは用意されているかもしれませんが。 ↩︎
見た目上の問題で、以下コードでは白黒逆にしました。 ↩︎