Merge pull request #2 from Mijadaj/conv

問題点の修正

_config.yml

@@ -0,0 +1,2 @@
+plugins:
+  - jekyll-sitemap

assets/style.css

@@ -30,7 +30,7 @@ li {
 .left {
     width: 50%;
     aspect-ratio: 4 / 3;
-    max-width: 800px;
+    max-width: 600px;
 }
 .right {
     width: 50%;
@@ -68,9 +68,6 @@ button, select {
     font-size: 16px;
     cursor: pointer;
 }
-select {
-    padding-right: 30px;
-}
 button:hover, select:hover {
     background-color: white;
 }

conv.html

@@ -20,10 +20,10 @@
         <a href="./" class="headerLink"><img src="./images/header.webp" alt="大札幌座標系" /></a>
     </header>
     </a>
-    <p>以下、道内各市街地の碁盤の目を基準とした測地基準系を <b>HGS</b> (Hokkaido Geodetic Systems) と呼ぶことにする。</p>
+    <p>以下，道内各市街地の碁盤の目を基準とした測地基準系を <b>HGS</b> (Hokkaido Geodetic Systems) と呼ぶことにする．</p>
     <h2>WGS から HGS への座標変換</h2>
     <p>
-        <b>WGS84</b> における地心直交座標を \((x, y, z)\), 緯度経度を \((\phi, \lambda)\), <b>HGS</b> における直交座標を \((X, Y, Z)\), 緯度経度を \((\Phi, \Lambda)\) で表す。最終的に求める条丁目は HGS の経緯度のみで決定し、楕円体高はすべて \(h = 0\) して考える。このとき、測地基準系 HGS を以下のように定義する。
+        <b>WGS84</b> における地心直交座標を \((x, y, z)\), 緯度経度を \((\phi, \lambda)\), <b>HGS</b> における直交座標を \((X, Y, Z)\), 緯度経度を \((\Phi, \Lambda)\) で表す．最終的に求める条丁目は HGS の経緯度のみで決定し，楕円体高はすべて \(h = 0\) して考える．このとき，測地基準系 HGS を以下のように定義する．
     </p>
     <p>\[
         \begin{array}{ll}
@@ -36,10 +36,10 @@ <h2>WGS から HGS への座標変換</h2>
     <p>\[
         \begin{eqnarray}\left\{ \begin{array}{ll}
         \Phi &= \sin^{-1} \dfrac{Z}{\sqrt{X^2 + Y^2 + Z^2}} \\
-        \Lambda &= \mathrm{sgn} (Y) \cos^{-1} \dfrac{X}{\sqrt{X^2 + Y^2}}
+        \Lambda &= \tan^{-1} \dfrac{Y}{X}
         \end{array}\right.\end{eqnarray}
     \]</p>
-    <p>ただし、\[
+    <p>ただし，\[
         \begin{array}{ll}
         \mathrm{O} (\phi_o , \lambda_o)&:&\text{HGS 経緯度原点の WGS 座標}\\
         \phi_o'&:&\text{地点Oにおける地心緯度}\\
@@ -49,11 +49,40 @@ <h2>WGS から HGS への座標変換</h2>
         R_x(\theta), R_y(\theta), R_z(\theta)&:& x軸, y軸, z軸\text{まわりの回転行列（反時計回りが正）}
         \end{array}
     \]</p>
-    <p>札幌座標系を例に考えれば、\(O (\phi_o, \lambda_o)\) は大通と創成川の交点の座標「北緯43度3分40.10秒 東経141度21分25.48秒」に対応し、\(\delta\) は真北に対する創成川の傾きである「西偏 10.5842138889 度」に対応する。</p>
-    <p>\((\Phi, \Lambda)\) から条丁目への変換は、あらかじめ定義された1区画あたりの緯度差\(\mathit{\Delta} \Phi\), 経度差 \(\mathit{\Delta} \Lambda\) にもとづいている。</p>
+    <p>札幌座標系を例に考えれば，\(O (\phi_o, \lambda_o)\) は大通と創成川の交点の座標「北緯43度3分40.10秒 東経141度21分25.48秒」に対応し，\(\delta\) は真北に対する創成川の傾き「西偏 10.5842138889 度」に対応する．</p>
+    <p>\((\Phi, \Lambda)\) から条丁目への変換には，あらかじめ定義された1区画あたりの緯度差\(\mathit{\Delta} \Phi\), 経度差 \(\mathit{\Delta} \Lambda\) を用いている．</p>
     <p>\[
-        \mathit{\Delta} \Phi = 
+        \begin{array}{cl}
+        \mathit{\Delta} \Phi &:= \dfrac{180}{\pi M(0)} l_\phi \\
+        \mathit{\Delta} \Lambda &:= \dfrac{180}{\pi N(0)} l_\lambda
+        \end{array}
     \]</p>
+    <p>ただし，\[
+        \begin{array}{cl}
+        M(\Phi) &:& \text{HGS} の緯度 \Phi における子午線曲率半径\\
+        N(\Phi) &:& \text{HGS} の緯度 \Phi における卯酉線曲率半径\\
+        l_\phi &:& 1区画の南北方向の長さ\\
+        l_\lambda &:& 1区画の東西方向の長さ\\
+        \end{array}
+    \]区画の1辺の長さは地球の半径に対して非常に小さいため，原点付近の緯度1度，経度1度あたりの長さを \(\pi M(0)/180, \pi N(0)/180 \) で近似し，\(l_\phi, l_\lambda\) をこれで除したものを1区画当たりの緯度差，経度差と定義した．なお各曲率半径の計算には Python の Numpy ライブラリを使用した．</p>
     <h2>HGS から WGS への座標変換</h2>
-    <p>上記の回転行列の逆行列 \(R_z(\lambda_o) R_y(-\phi_o') R_x(\delta) \) を用いる。</p>
+    <p>上記の回転行列の逆行列 \(R_z(\lambda_o) R_y(-\phi_o') R_x(\delta) \) を用いる．上述のように \(h=0\) であり経緯度さえ求まれば良いので，目的の地点を指す単位ベクトルを用いて計算する．</p>
+    <p>\[
+        \begin{array}{ll}
+        \begin{bmatrix}n_x \\ n_y \\ n_z\end{bmatrix}
+        &= R_x(-\delta) R_y(\phi_o') R_z(-\lambda_o) \begin{bmatrix}n_X \\ n_Y \\ n_Z\end{bmatrix} \\
+        &= R_x(-\delta) R_y(\phi_o') R_z(-\lambda_o) \begin{bmatrix}\cos \Phi \cos \Lambda \\ \cos \Phi \sin \Lambda\\ \sin \Phi \end{bmatrix}
+        \end{array}
+        \]
+    </p>
+    <p>\[
+        \begin{eqnarray}\left\{ \begin{array}{ll}
+        \phi &= \sin^{-1} \dfrac{n_z}{\sqrt{(1 - e^2)^2 (n_x^2 + n_y^2) + n_z^2}} \\
+        \lambda &= \tan^{-1} \dfrac{n_y}{n_x}
+        \end{array}\right.\end{eqnarray}
+    \]</p>
+    <footer>© 2024
+        <a href="https://github.com/Mijadaj" target="_blank">Міја</a>;
+        <a href="https://github.com/Mijadaj/Sapporo-Coordinate-System" target="_blank">GitHub</a>
+    </footer>
 </body>

google96ced66e8ce92e16.html

@@ -0,0 +1 @@
+google-site-verification: google96ced66e8ce92e16.html

index.html

@@ -25,33 +25,33 @@
         <option id="initialOption" value="">--座標系を変更--</option>
     </select>
     <button class="head" onclick="window.open('./conv.html')">計算式</button>
-    <p>大通を赤道、創成川（南6条～北11条間, 真北に対して西に10.6度の傾き）を本初子午線とする座標系をもとに地球上のあらゆる地点を札幌の住所表示に変換します。</p>
+    <p>大通を赤道，創成川（南6条～北11条間, 真北に対して西に10.6度の傾き）を本初子午線とする座標系をもとに地球上のあらゆる地点を札幌の住所表示に変換します．</p>
     <div class="container">
         <div class="left" id="map"></div>
         <div class="right">
             <table id="converter">
                 <caption>座標変換</caption>
                 <tr>
-                    <th>緯度, 経度（10進法）</th>
+                    <th>緯度, 経度 (10進法)</th>
                     <td><input type="text" id="latLng" placeholder="緯度, 経度"></td>
                 </tr>
             </tr>
             <th>条丁目</th>
             <td>
                 <div id="jochome">
                     <select class="jochome" id="selectJo"></select>
-                    <input type="text" class="jochome" id="joNumber"><span id="jo">条</span>
+                    <input type="number" class="jochome" id="joNumber"><span id="jo">条</span>
                     <select class="jochome" id="selectChome"></select>
-                    <input type="text" class="jochome" id="chomeNumber"><span>丁目</span>
+                    <input type="number" class="jochome" id="chomeNumber"><span>丁目</span>
                 </div>
             </td>
             <tr>
             </table>
             <ul id="description">
-                <li>地図上のピンを動かすとその位置の経緯度、および選択した座標系における条丁目が表示されます。ピン、経緯度、条丁目は連動しています。</li>
-                <li>赤い点は選択した座標系の原点です。</li>
-                <li>ピンの初期位置は各座標系のもととなった自治体の姉妹都市に設定してあります。</li>
-                <li>この座標変換における条丁目は楕円体座標に基づいています。平面直角座標ではありません。したがって南北方向の条（丁目）の数が大きくなるほど、東西方向の丁目（条）の間隔は狭くなります。札幌座標系の場合、北76919条（アイルランド沖 南西約1300kmの海域）付近ですべての丁目が1点に集まります。</li>
+                <li>地図上のピンを動かすとその位置の経緯度，および選択した座標系における条丁目が表示されます．ピン，経緯度，条丁目は連動しています．</li>
+                <li>赤い点は選択した座標系の原点です．</li>
+                <li>ピンの初期位置は各座標系のもととなった自治体の姉妹都市に設定してあります．</li>
+                <li>この座標変換における条丁目は楕円体座標に基づいています．平面直角座標ではありません．したがって南北方向の条（丁目）の数が大きくなるほど，東西方向の丁目（条）の間隔は狭くなります．札幌座標系の場合，北76919条（アイルランド沖 南西約1300kmの海域）付近ですべての丁目が1点に集まります．</li>
             </ul>
         </div>
     </div>

js/conv.js

@@ -2,7 +2,7 @@ const hgsSet = {
     Sapporo: {
         name: "大札幌 座標系",
         format: [["北", "南", "大通"], ["東", "西"]],
-        joOffset: ["latitudinal", 175, 110],
+        joOffset: ["latitudinal", -0.000990033565819702, 0.00157505340016771], //30.8631511197875 m/sec
         origin: [43.061138, 141.357079],
         declination: 10.5842138889,
         blockSizeLatLng: [0.001170046801872075, 0.0011661354271942782],
@@ -11,7 +11,7 @@ const hgsSet = {
     Obihiro: {
         name: "帯広 座標系",
         format: [["西", "東", "大通"], ["南", "北"]],
-        joOffset: ["longitudinal", 69, 69],
+        joOffset: ["longitudinal", -0.000618898932322025, 0.000618898932322025], //30.9689767839087 m/sec
         origin: [42.932044, 143.204010],
         declination: 5.59292777778,
         blockSizeLatLng: [0.001170168504264614, 0.0011660447761194029],
@@ -20,7 +20,7 @@ const hgsSet = {
     Nayoro: {
         name: "名寄 座標系",
         format: [["西", "東", "大通"], ["南", "北"]],
-        joOffset: ["longitudinal", 74, 74],
+        joOffset: ["longitudinal", -0.000663668630632269, 0.000663668630632269], //30.9726188745316 m/sec
         origin: [44.356392, 142.463787],
         declination: 0.0992,
         blockSizeLatLng: [0.001169909916936396, 0.0011659012740710033],
@@ -29,7 +29,7 @@ const hgsSet = {
     Nakashibetsu: {
         name: "中標津 座標系",
         format: [["東", "西", "大通"], ["南", "北"]],
-        joOffset: ["longitudinal", 45.5, 45.5],
+        joOffset: ["longitudinal", -0.000408735777291158, 0.000408735777291158], //30.9219050327609
         origin: [43.548025, 144.973040],
         declination: 42.1626611111,
         blockSizeLatLng: [0.0008177655011993894, 0.0008174683912222061],
@@ -39,6 +39,7 @@ const hgsSet = {
     template: {
         name: " 座標系",
         format: [[, ], [, ]],
+        joOffset: [],
         origin: [, ],
         declination: ,
         blockSizeLatLng: ,
@@ -73,7 +74,7 @@ const pi = acos(-1);
 const radians = (deg) => Decimal.div(pi, 180).mul(deg);
 const degrees = (rad) => Decimal.div(180, pi).mul(rad);
 
-let latLng = [];
+let latLng;
 let blockPosition = [];
 let jochome = [];
 let coordinateSystem = "Sapporo";
@@ -128,17 +129,62 @@ function hgs2wgs([LAT, LNG]) {
 
 //HGS経緯度から区画の位置を算出
 function hgs2blockPosition([LAT, LNG]) {
-    let latitudinalBlock = Decimal.abs(LAT).div(hgsSet[coordinateSystem].blockSizeLatLng[0]).floor().add(1);
-    let longitudinalBlock = Decimal.abs(LNG).div(hgsSet[coordinateSystem].blockSizeLatLng[1]).floor().add(1);
-    latitudinalBlock = Decimal(Math.sign(LAT)).mul(latitudinalBlock);
-    longitudinalBlock = Decimal(Math.sign(LNG)).mul(longitudinalBlock);
-    blockPosition = [latitudinalBlock, longitudinalBlock]
+    const offset = hgsSet[coordinateSystem].joOffset //大通の幅に対する補正
+    const latDegPerBlock = hgsSet[coordinateSystem].blockSizeLatLng[0];
+    const lngDegPerBlock = hgsSet[coordinateSystem].blockSizeLatLng[1];
+    const adjustJo = function(deg, degPerBlock) {
+        let pos;
+        if (deg <= offset[1]) {
+            pos = sub(deg, offset[1]).div(degPerBlock).ceil().sub(1);
+            return pos;
+        } else if (offset[2] <= deg) {
+            pos = sub(deg, offset[2]).div(degPerBlock).floor().add(1);
+            return pos;
+        } else {
+            return 0;
+        }
+    };
+    let latitudinalBlock, longitudinalBlock;
+    if (offset[0] == "latitudinal") {
+        latitudinalBlock = adjustJo(LAT, latDegPerBlock);
+        longitudinalBlock = Decimal.abs(LNG).div(lngDegPerBlock).floor().add(1);
+        longitudinalBlock = Decimal(Math.sign(LNG)).mul(longitudinalBlock);
+    } else {
+        latitudinalBlock = Decimal.abs(LAT).div(latDegPerBlock).floor().add(1);
+        latitudinalBlock = Decimal(Math.sign(LAT)).mul(latitudinalBlock);
+        longitudinalBlock = adjustJo(LNG, lngDegPerBlock);
+    };
+    blockPosition = [latitudinalBlock, longitudinalBlock];
     return blockPosition;
 };
 
 //区画の位置からHGS経緯度を算出
 function blockPosition2hgs([latitudinalBlock, longitudinalBlock]) {
-    const LAT = new Decimal(latitudinalBlock).mul(hgsSet[coordinateSystem].blockSizeLatLng[0]);
-    const LNG = new Decimal(longitudinalBlock).mul(hgsSet[coordinateSystem].blockSizeLatLng[1]);
+    const offset = hgsSet[coordinateSystem].joOffset; //大通の幅に対する補正
+    const latDegPerBlock = hgsSet[coordinateSystem].blockSizeLatLng[0];
+    const lngDegPerBlock = hgsSet[coordinateSystem].blockSizeLatLng[1];
+    const restoreDeg = function(jo, degPerBlock) {
+        let deg;
+        switch (Math.sign(jo)) {
+            case 1:
+                deg = sub(jo, 1).mul(degPerBlock).add(offset[2])
+                break
+            case -1:
+                deg = sub(jo, -1).mul(degPerBlock).add(offset[1])
+                break
+            case 0:
+                deg = 0
+                break
+        }
+        return deg;
+    };
+    let LAT, LNG;
+    if (offset[0] == "latitudinal") {
+        LAT = restoreDeg(latitudinalBlock, latDegPerBlock);
+        LNG = Decimal(longitudinalBlock).mul(lngDegPerBlock);
+    } else {
+        LAT = Decimal(latitudinalBlock).mul(latDegPerBlock);
+        LNG = restoreDeg(longitudinalBlock, lngDegPerBlock);
+    };
     return [LAT, LNG]
 };