💧

「日本で最も丸い湖」を探して - GIS, 円形度, 最急降下法

sorami

2024/04/04に公開

gis

openstreetmap

tech

はじめに

温泉で有名な北海道の登別（のぼりべつ）には「のぼりべつクマ牧場」というものがあり、そこから隣の白老（しらおい）にある 倶多楽湖（くったらこ） という湖を望むことができます。

私はこれまで二度、クマ牧場を訪れましたが、一度目は霧で隠れており、二度目にその姿を目にすることができました。

のぼりべつクマ牧場から見た倶多楽湖 - 2023年4月

のぼりべつクマ牧場から見た倶多楽湖 - 2024年2月

この倶多楽湖は 日本で最も丸い と言われているそうです。

倶多楽湖は直径3キロメートル、周囲約8キロメートルの小さな丸いカルデラ湖。その丸さは日本で最も円に近いといいます。

倶多楽湖（くったらこ）／白老町 | たびらい

たしかに、地図で見ると、とても丸い感じがします。

OpenStreetMapでの倶多楽湖

しかし、本当にこれより丸い湖はないのでしょうか？

そもそも「丸い」とは何か

札幌には地元民から「まるいさん」という愛称で親しまれている百貨店がありますが、そもそも、「丸い」とはどう定義されるのでしょうか。

真円度

まず、英語版Wikipediaにある「Roundness」という記事を参照してみましょう。

https://en.wikipedia.org/wiki/Roundness

ここでは、ISOでの定義として「図形の”内接円”と”外接円”の半径の比率」というものが紹介されています。

「内接円」（inscribed circle）、つまり図形にすっぽり収まる円のうち最大のもの（その多角形の内部にあり全ての辺に接する円）、そして「外接円」（circumscribed circle）、すなわち図形をすっぽりと収めることができる円のうち最小のもの（その多角形の全ての頂点を通る円）、これらの大きさの比をもとにしています。

ISO Roundness

図: Claudio Rocchini, CC BY-SA 4.0 DEED, File:Iso roundness.svg - Wikipedia

上記の例では、左に示した「四角形」のroundnessは、 ${{\frac {1}{\sqrt {2}}}\simeq 0.7}$ 、右にある「八角形」のroundnessは ${{\frac {\sqrt {2+{\sqrt {2}}}}{2}}\simeq 0.92}$ ということになり、八角形の方が「丸い」ということになります。

真円のとき、このroundnessは $0$ となります。

この定義は、日本語では「真円度」と呼ばれているようです。JISB0621:1984 幾何偏差の定義及び表示を参照してみましょう:

4.3 真円度真円度とは，円形形体の幾何学的に正しい円（以下，幾何学的円という。）からの狂いの大きさをいう。

p.2

5.3 真円度真円度は，円形形体 (C) を二つの同心の幾何学的円で挟んだとき，同心二円の間隔が最小となる場合の，二円の半径の差 (f) で表し（図6），真円度＿mm又は真円度＿μmと表示する

p.4

円形度

一方で、別の尺度として「円形度」というのもあるそうです。以下のように定義されています:

円形度 = 4\pi \frac{\text{面積}}{\text{周囲長}^2}

円形度が1に近いほど、円らしいと言えます。

円の性質として「周長（図形の周りの長さ）が同じとき、その面積がもっとも大きくなる図形である」ということが知られています。

円の面積は、 $\pi r^2$ で求められます（ $r$ が半径）。円周は、 $2 \pi r$ です。これを先ほどの式に当てはめると、真円の円形度は $1$ になります。

真円度と円形度

では、真円度と円形度という二つの尺度、どちらの方が我々の目的（丸い湖を探したい）に適しているでしょうか。

次の記事では、いくつかの図形で、それぞれの真円度と円形度を算出しています:

「真円度」と「円形度」で図形を数値化した例

図: 【ImageJで円形度解析】見た目の「丸さ」を表すには、「真円度」よりも「円形度」？ - 60へぇ研究所より引用

真円度を見ると、中央と一番右が同等の値となっています。

円形度では、右に行くにつれ数値が下がる、つまり円らしくない、ということになります。

どちらが正しいということはなく、それぞれの特性がある、ということが見て取れます。

ただ、今回の用途においては、時に複雑な湖の形に対し、そもそも、真円度の算出に必要な内接円と外接円が存在しないことも多いでしょう。真円度は、円形にある程度近しい、複雑ではない図形に適しているもので、それがどれくらい真円から歪んでいるかを測るための尺度なのかなと思います。

一方、円形度は、面積と周長さえ分かれば、シンプルに算出できます。

円形度で湖を見る

それではこの円形度を手がかりに、日本にある湖の丸さを調べてみましょう。

処理は全てPythonで行いました。地理空間情報を扱うライブラリGeoPandasを用いています。以下のJupyter Notebookで、手順を確認することができます:

roundest-lake/japan.ipynb at main · sorami/roundest-lake

データ

国土数値情報の「湖沼データ」を利用しました:

「数値地図25000（地名・公共施設）」に含まれる約600の湖、日本の湖沼アトラス（国土地理院技術資料D1-No,299）に含まれる58の自然湖沼、「国土数値情報・統一フォーマット（湖沼、湖岸線）」に含まれる約300湖沼について、その形状を数値地図25000（空間データ基盤）より取得し、属性情報として湖沼名を数値地図25000（地名・公共施設）の注記、国土数値情報統一フォーマット（湖沼、湖岸線）の湖沼名、「日本の湖沼アトラス」の湖沼名から付与した。さらに「日本の湖沼アトラス」に記載された自然湖沼については、最大水深・水面標高を属性として付与した。

国土数値情報「湖沼データ」をQGISで表示した図

このデータセットには556の湖沼が含まれていました。そのうち107（20%弱！）が、北海道にあるものでした。

円形度の算出

元データの座標系は「JGD2000 / （B, L）」（「B, L」 = 「地理座標系（緯度経度）」）、すなわちEPSG 4612に該当します。

これを、それぞれの面積と周囲長を算出するために、地理座標系から投影座標系へ投影変換します。今回は簡単のため全てに対し「JGD2000 / UTM zone 53N」（EPSG 3099）を用いました。

面積と周囲長は、GeoPandasの内部で用いられているライブラリShapelyが算出してくれます。

面積の最も大きく、周囲長の最も長い湖沼は、滋賀県の琵琶湖で、データ上では、面積669km², 周囲長260kmです。

面積の最も小さく、周囲長の最も短い湖沼は、青森県の青池（十二湖）で、面積429m²、周囲長92mでした。面積は琵琶湖の150万分の1です。

あとはこれらの値を、先述した円形度の式に当てはめるだけです:

gdf["円形度"] = 4 * math.pi * gdf["面積"] / (gdf["周囲長"] ** 2)

円形度による「丸い湖」

以下が、円形度の値が最も大きい湖トップ5です:

0.954554 - 明月湖（山形県）
0.945450 - 倶多楽湖（北海道）
0.920083 - 東雲湖（東小沼）（北海道）
0.907386 - 四尾連湖（山梨県）
0.902002 - 巣鷹湖（長野県）

... なんということでしょう、冒頭で述べた我らが倶多楽湖が、1位ではないではありませんか！

この明月湖が、日本で一番丸い湖なのでしょうか。福島県と山形県の県境に位置する吾妻山にあるものだそうです（参考: 明月湖 -日本湖沼めぐり-）。

OpenStreetMapでの明月湖

見てみると、確かに、丸いですが、少しひしゃげているようにも思えます。

再び、倶多楽湖を見てみましょう。

OpenStreetMapでの倶多楽湖

なんとなく、倶多楽湖のほうが、「丸い」気がする！これは贔屓からだけではないように思えます。

そもそも、明月湖は、面積978m², 周囲長113m。倶多楽湖は、面積5km², 周囲長8kmです。

縮尺が違うため、ポリゴンの”細かさ”も違います。そのため、面積と周囲長だけをもとに算出される「円形度」という値が、単純に比較できるものではなく、小さな湖の方が有利になっているのではないでしょうか ...

「国」のかたち

ここでは話を変えて、先行事例として、国（国土）の「四角さ」「丸さ」が検討されたものを見ていきます。

国の四角さ（Barry, 2016）

David Barryによる "The rectangularness of countries" という考察では、国の「四角さ（Rectangularness）」を検討しています:

https://pappubahry.com/misc/rectangles/

「トルコはとても四角い」と聞いたことから調べてみたそうです（結果としてトルコは15位、1位はエジプト、日本は187位）。

ここでは四角さを「”国土と同じ面積の四角形”との重なりの最大割合」と定義しています。

しかし、「同じ面積の四角形」の中から、国土との重なりが”最大”になるものはどう特定できるでしょうか。

Barryは、シンプルな最急降下法により、「重なりの割合」、つまり「四角さのスコア」が最大になる四角形を探索しています。

まず、可能なすべての「同じ面積の四角形」を、以下4つのパラメーターで表します:

中心点の緯度
中心点の経度
回転度
アスペクト比

そして、現在のパラメーターを少し変更して、スコアへの影響を見る。もしスコアが上昇していたら、変更後の値を現在のパラメーターに設定する。これを、スコアが収束する（それ以上にならなくなる）か、規定の回数に達するまで繰り返します。

初期値として、中心地の経緯度には「国土ポリゴンの重心」、回転度は「0」、アスペクト比は「バウンディングボックス（ポリゴンを含む最小の四角）での比率」を設定したそうです。

以下に、本人による解説ページで紹介されている、コンゴ共和国でのパラメーター探索の様子を示します:

コンゴ共和国でのパラメーター探索の様子

また、本人による解説とR言語でのコードで、詳細が述べられています。

大部分のケースでは、このシンプルな最適化手法でうまく行ったそうですが、いくつかの問題も生じたそうです。

例えば、沢山の島がある場合には、重心点が陸地から遠くなりすぎて、重なりがゼロになります。そのときには、それらを含む一番大きな四角形を手動で初期値に設定したそうです。

最急降下法では200イテレーションを行ったそうですが、うまく収束しないケース（43カ国）では、アスペクト比の更新が遅かったようです。細長い国の場合に、島があったりすると、アスペクト比の初期値が小さくなるので、最適な値までイテレーション内に到達できなかったとのこと。収束しなかった場合は、更新値を少し大きくして再度実行したそうです。

元データはNatural Earthで、これをシンプルに正距円筒図法（x=経度, y=緯度）へ投影変換して用いています（理想的には国ごとに、適した投影法が適用できると良い）。

国の丸さ（Ciruelos, 2016）

上記のBarryによる「四角い国」の考察を踏まえて、Gonzalo Ciruelosによる "What is the roundest country?" という論説では、「国の丸さ」に着目しています:

まず、「丸さ」を定義する必要があります。

円形度についての節で述べたように、円の性質として「周長（図形の周りの長さ）が同じとき、その面積がもっとも大きくなる図形である」ということが知られています。

しかし「国」の場合は、周囲（国境）が細かくギザギザとしていることが多いため、想定以上に長い周長となり、この円の性質を利用するのには向いていないでしょう。

そのためCiruelosは、任意の国（2次元空間上の図形）を $C \subset \mathbb{R}^2$ としたとき、以下のように丸さを定義しています:

\text{roundness}(C) = \max\limits_{x \in \mathbb{R}^2, r \in \mathbb{R}_{\gt 0}} \frac{ \text{area}(C \cap D(p, r))} {\max\{\text{area}(D(p, r)), \text{area}(C)\}}

（元記事にある $x$ を、わかりやすさのため、 $p$ と書き換えて記述しています）

ここで $D(p, r)$ は、中心点 $p$ 、半径 $r$ の円（Disk）を表します。

分子 $\text{area}(C \cap D(p, r))$ は、国 $C$ と任意の円 $D$ の重なった面積を表します。

分母 $\text{area}(D(p, r)), \text{area}(C)\}$ は、国 $C$ か任意の円 $D$ 、大きい方の面積を表します。

このスコアは、0から1の間の値をとり、 $C$ が円の時に1となります。

丸さのスコアを得るためには、これらの比率が最大になる円のパラメーター $x$ と $r$ を探す必要があります。

ここでCiruelosは、Barryが四角さを求める際にやったのと同様に、最急降下法を用いてパラメーターを求めています。

算出のためのPythonコードは、GitHub gistで公開されています:

バウンディングボックス（bbox）は、その図形を囲む最小の長方形です。それをもとに、まず初期値として、以下の二つを設定しています:

bboxの内に収まる円のうち、面積が最大のもの（半径 = bboxの短辺 / 2）
bboxを収める円のうち、面積が最小のもの（半径 = bboxの対角線 / 2）

初期値の円: 内側（直径 = バウンディングボックスの短辺）

初期値の円: 外側（直径 = バウンディングボックスの対角線）

コードの該当部分:

def circle_inside(xmin, ymin, xmax, ymax):
    center = [(xmin + xmax) / 2, (ymin + ymax) / 2]
    radius = min((xmax - xmin) / 2, (ymax - ymin) / 2)
    return center, radius

def circle_outside(xmin, ymin, xmax, ymax):
    center = [(xmin + xmax) / 2, (ymin + ymax) / 2]
    radius = math.sqrt((xmax - xmin)**2.0 + (ymax - ymin)**2.0) / 2
    return center, radius

その上で、円を上下左右に移動させ、また半径を大きく・小さくすることで、重なりが大きくなるよう、パラメーターを最適化していきます。

元データはNatural Earthで、WGS84（経緯度）の形式です。これを正距円筒図法にしています。しかし、多くの国が歪んでしまうため、そこからさらに、国ごとの方位図法へ変換しています。方位図法の中心としては、その国のどこか（具体的には、その国のボーダーからランダムな2点を選び、その中間点）を設定しています。

また、このgistへのJan Freybergによるコメントも少々興味深いものでした。コメントでは、国の「丸さ」と「経度」（首都の経度の絶対値）が相関していると指摘されています。おそらく、太平洋にある島国（領土がすべて島から構成される国）が影響しているのだろうと述べられていました:

Scatterplot - Roundness of Country vs. Longitude of Capital

改めて、丸い湖を探す

さて、日本の湖に戻ります。

先の試行では、面積と周囲長から解析的に求める円形度を用いましたが、これは縮尺によって単純な比較ができなさそうであるという問題がありました。

では、Ciruelosによる例を参考に、湖の丸さを、数値的に見てみましょう。湖と円との重なり割合を求めます。

最急降下法による探索

以下のJupyter Notebookで、処理を見ることができます:

roundest-lake/japan_gd.ipynb at main · sorami/roundest-lake

基本的には、Ciruelosと同様の考え方で進めています。

以下のように、円が動き、大きさが変わり、スコアが上昇していきます。

北海道のシュンクシタカラ湖:

1970年代に、人工衛星によってその存在が確認された湖で、2007年現在、日本国内最後に発見された湖とされる。しかし、地元の人々はそれ以前より湖が存在することを知っており、1920年の「国土地理院5万分の1地図」にも「シュンクシタカラ沼」と表記がある[1]。

シュンクシタカラ湖 - Wikipedia

シュンクシタカラ湖 - パラメーター最適化の様子

北海道のコムケ湖:

コムケ湖 - パラメーター最適化の様子

三日月湖: 初期値の調整

最初、円の初期値として、中心点は「バウンディングボックスの中心」に、半径は「湖の面積」と同じになるよう設定してみました。

しかし、多くのケースではこれでも上手くいくのですが、例えば北海道の「しのつ湖」では、スコアがゼロになってしまいます。三日月湖のような凹図形では、中心点が領域外に配置され、また半径が小さすぎて円が全く重ならないためです:

しのつ湖 - スコアが0になっている

これは、初期値の半径を、湖の面積から求めるのではなく、「バウンディングボックスを囲む円」（半径 = バウンディングボックスの対角線 / 2）とすることで、局所で収束せずに、それなりの値になりました:

def get_initial_circle_params(lake):
    """
    探索の初期値

    中心: バウンディングボックスの中心
    半径: バウンディングボックスの対角線の長さの1/2
    """
    x_min, y_min, x_max, y_max = lake.bounds
    center = [(x_min + x_max) / 2, (y_min + y_max) / 2]

    radius = math.sqrt((x_max - x_min) ** 2 + (y_max - y_min) ** 2) / 2

    return center, radius

しのつ湖 - パラメーター最適化の様子

ちなみに、しのつ湖は冬には凍って、その上でワカサギ釣りができます。私は、ここで人生初のワカサギ釣りをやりました。朝一番に4人で行って、3時間かけて、4匹しか釣れませんでした。

しのつ湖でのワカサギ釣りの様子

ダム湖: 更新幅の調整

円を「どれくらい」移動させたり、大きく・小さくさせるか。

当初は「半径の1%」をベースに、その値を足したり引いたりして、移動と拡大縮小を行っていました:

def neighbors(center, radius, delta_perc=0.01):
    delta = radius * delta_perc
    delta_options = [-delta, 0, delta]

    for cx_delta in delta_options:
        for cy_delta in delta_options:
            for r_delta in delta_options:
                yield ((center[0] + cx_delta, center[1] + cy_delta), radius + r_delta)

しかしこれでは、局所最適解に陥ることがありました。以下は、沖縄にある「伊集の湖」の例です:

伊集の湖 - 局所最適解になっている

このようなダム湖は、細々としており、円がうまく当てはまりません。そして、端の方に大きな領域がある場合、更新幅が小さすぎると、そこまで探索が辿り着かずに収束してしまいます。

そのため、大きな更新幅を用意（具体的には「半径の0.5%」に加えて「半径の10%」も考慮）することで、局所最適へ陥りづらいようにしました:

def neighbors(center, radius, delta_perc=0.005):
    delta = radius * delta_perc
    delta_options = [-delta * 20, -delta, 0, delta, delta * 20]

    for cx_delta in delta_options:
        for cy_delta in delta_options:
            for r_delta in delta_options:
                yield ((center[0] + cx_delta, center[1] + cy_delta), radius + r_delta)