日本一マクドナルドから遠い場所　をより精確に求めてみた

# 経緯と目的

j_catfish氏の[日本一マクドナルドから遠い場所](https://qiita.com/j_catfish/items/0d36a21b4744aa3e4870)の面白い着想に触発され、地図大好き人間として自分も書きたくなった。

目的としては以下の通り。

* Google Maps APIは地球を回転楕円体ではなく球として距離の計算をするので、最大0.3%程度の誤差が出る
* この誤差がもしかしたら結果に影響を与えるかもしれない、せたな町と知床の距離の検証を行う

# データの用意

今回用意したデータは以下の2つ。

* [基盤地樹情報ダウンロードサービス](https://fgd.gsi.go.jp/download/menu.php)から海岸線データ
    * せたな町の一部(2次メッシュ番号で633926、633936、633946の3箇所)
    * 知床半島の先端部(2次メッシュ番号で664532, 664542の2箇所)
* [日本マクドナルドの店舗検索ページ](https://map.mcdonalds.co.jp/)から、店舗の位置情報JSONを開発者ツールを使って保存
    * せたな町近辺のマクドナルド店舗のJSON(店舗数68)
    * 知床半島近辺のマクドナルド店舗のJSON(店舗数26)

計算を高速化するため、範囲を絞った。


ちなみに、具体的にどの店舗なのか不明なのだが、元記事の時点ではマクドナルドの総店舗数が2887だったのが、2018年5月2日時点では1店舗減って2886になっていた。
(追記: 2887店に戻ってました)

# 計算

測地線の距離計算には、誤差15nmと言われる[GeographicLib](https://geographiclib.sourceforge.io/)をJ/Link経由でMathematicaから呼び出す、拙作のパッケージ[mathematica-geodesic](https://github.com/330k/mathematica-geodesic/)を使用した。

## せたな町

まずは海岸線データのXMLファイルを読み込み、`posList`要素に記述された緯度経度データを抜き出す。
海岸線データには岩礁なども含まれるため、始点と終点の座標が同じになっているものは除去した。

なお海岸線データの`Curve`要素は、海岸線の管理者(国か県か)といった行政の都合によってバラバラに分割されていて非常に使いにくい。
そのため完全には除去できていないのだが、今回は最遠点を求めるのに支障がなければよしとする。

~~~Wolfram
xmlfiles1 = FileNames["data00/FG-GML-63*.xml"];
xmldata1 = Import /@ xmlfiles1;

curves1 = 
  Flatten[Cases[#, XMLElement[{_, "Curve"}, _, _], Infinity] & /@ 
    xmldata1, 1];

coastpoints1 = 
  Flatten[Select[
    Partition[
       Flatten[Cases[#, 
         XMLElement[{_, "posList"}, _, 
           a_] :> (ToExpression /@ StringSplit[a, Whitespace]), 
         Infinity]], 2] & /@ curves1, First@# != Last@# &], 1];
~~~

抜き出した緯度経度データの個数は5990。

~~~Mathematica
In[5]:= Length[coastpoints1]
Out[5]= 5990
~~~

次にせたな町周辺のマクドナルド店舗のJSONから、緯度経度データを抜き出す。

~~~Mathematica
mcdonalds1 = {"latitude", "longitude"} /. Import["mcdonalds_setana.json"];
~~~

こちらのデータの個数は68。

~~~Mathematica
In[7]:= Length[mcdonalds1]

Out[7]= 68
~~~

海岸線の緯度経度データとマクドナルド店舗の緯度経度データのすべての組み合わせの距離を求める。
一般化した外積`Outer`を使うと一発で楽ちん。

~~~Mathematica
<< Geodesic`
dists1 = Min /@ Outer[GeoDistanceGL, coastpoints1, mcdonalds1, 1];
~~~

さて・・・結果は？

~~~Mathematica
In[17]:= Max[dists1]

Out[17]= 93357.2

In[24]:= coastpoints1[[Ordering[dists1, -1]]]

Out[24]= {{42.3029, 139.767}}
~~~

緯度経度はFullFormで精確な値を出して、せたな町のマクドナルド最遠点の位置は

* 緯度:  42.302894444°
* 経度: 139.766917778°
* 距離: 93.357 km

尾花岬の北にも2つほど突端部があるが、そこではなく尾花岬。

![image.png](https://qiita-image-store.s3.amazonaws.com/0/197772/8727dbc9-ea08-b281-008c-88ad98bf1bb5.png)

## 知床

知床半島突端のデータに対しても同様に計算。

~~~Mathematica
In[10]:= xmlfiles2 = FileNames["data00/FG-GML-66*.xml"];

In[11]:= xmldata2 = Import /@ xmlfiles2;

In[12]:= coastpoints2 = 
 Flatten[Select[
   Partition[
      Flatten[Cases[#, 
        XMLElement[{_, "posList"}, _, 
          a_] :> (ToExpression /@ StringSplit[a, Whitespace]), 
        Infinity]], 2] & /@ curves2, First@# != Last@# &], 1]

In[13]:= Length[coastpoints2]

Out[13]= 6375

In[14]:= mcdonalds2 = {"latitude", "longitude"} /. 
   Import["mcdonalds_shiretoko.json"];

In[15]:= Length[mcdonalds2]

Out[15]= 26

In[18]:= dists2 = 
  Min /@ Outer[GeoDistanceGL, coastpoints2, mcdonalds2, 1];

In[19]:= Max[dists2]

Out[19]= 92740.

In[25]:= coastpoints2[[Ordering[dists2, -1]]]

Out[25]= {{44.3389, 145.345}}
~~~

ということで、こちらは以下のようになった。

* 緯度:  44.338895556°
* 経度: 145.344701944°
* 距離: 92.740 km

せたな町のほうが遠い。その差は617m。

ちなみに位置は、灯台の東側。
![image.png](https://qiita-image-store.s3.amazonaws.com/0/197772/419c1503-d3c4-7bef-0177-7e30680b141d.png)


# 結論

日本一マクドナルドから遠い場所は、やっぱり***北海道せたな町の尾花岬***！

![image.png](https://qiita-image-store.s3.amazonaws.com/0/197772/0f8fb719-6fc7-53b7-7aec-f44021e90df4.png)

航空写真を見るととても徒歩では到達できそうにない。

![image.png](https://qiita-image-store.s3.amazonaws.com/0/197772/3e62bccd-a2f5-5e7d-5c81-8595aa3b62e7.png)