2018年1月15日月曜日

メルカトル図法の投影式からWebメルカトルの投影式を導く

地球を球体とみなしたときのメルカトル図法の投影式は次の通り。

\begin{align} x &= R(\lambda - \lambda_0) \\ y &= R \ln \left[\tan \left(\frac{\pi}{4} + \frac{\varphi}{2} \right)\right] \end{align}

\(y\)座標の導出は地図における数学が詳しい。\(x\)座標は引き伸ばしながら平面に展開するだけで、\(\lambda_0\)は地図の中央の経度。

Webメルカトルでは\(\lambda_0 = 0\)なので代入する。

\[x = R\lambda\]

このとき\(-\pi \le \lambda \le \pi\)より\(-R \pi \le x \le R \pi\)。

また、Webメルカトルは北極と南極周辺を無視することで正方形に投影するため、\(y\)座標のとる値の範囲は\(x\)座標と等しく\(-R \pi \le y \le R \pi\)。\(x \ge 0\), \(y \ge 0\)にしたいので、ともに\(R \pi\)を加え、\(y\)軸は下方向が正だから符号を反転すると、

\begin{align} x &= R (\lambda + \pi) \\ y &= R \left(\pi - \ln \left[\tan \left(\frac{\pi}{4} + \frac{\varphi}{2} \right)\right]\right) \end{align}

となり、\(0 \le x \le 2R \pi\), \(0 \le y \le 2R \pi\)。

Webメルカトルでは右下隅の座標が\((256, 256)\)になるので、\(2R \pi\)で割り、\(256\)を掛ける。

\begin{align} x &= \frac{128}{\pi} (\lambda + \pi) \\ y &= \frac{128}{\pi} \left(\pi - \ln \left[\tan \left(\frac{\pi}{4} + \frac{\varphi}{2} \right)\right]\right) \end{align}

\(0 \le x \le 256\), \(0 \le y \le 256\)となって、Webメルカトルの投影式が導かれた。

なお、Showing pixel and tile coordinates | Google Maps JavaScript API | Google Developersproject()や、google.maps.Projection.fromLatLngToPoint()で使われている投影式は、形が違うが変形すれば同じものになる。

\[ \tan(\alpha + \beta) = \frac{\sin(\alpha + \beta)}{\cos(\alpha + \beta)} = \frac{\sin \alpha \cos \beta + \cos \alpha \sin \beta}{\cos \alpha \cos \beta - \sin \alpha \sin \beta} \]

に\(\alpha = \frac{\pi}{4}\)を代入する。また、あとで代入する\(\beta = \frac{\varphi}{2}\)と、\(-\frac{\pi}{4} \lt \frac{\varphi}{2} \lt \frac{\pi}{4}\)から\(\cos \beta + \sin \beta \neq 0\)。

\begin{align} \tan \left(\frac{\pi}{4} + \beta \right) &= \frac{\cos \beta + \sin \beta}{\cos \beta - \sin \beta} \\ &= \frac{(\cos \beta + \sin \beta)(\cos \beta + \sin \beta)}{(\cos \beta - \sin \beta)(\cos \beta + \sin \beta)} \\ &= \frac{1 + 2 \sin \beta \cos \beta}{\cos^2 \beta - \sin^2 \beta} \\ &= \frac{1 + \sin 2 \beta}{\cos 2 \beta} \end{align}

\(\beta = \frac{\varphi}{2}\)も代入すると、

\begin{align} \tan \left(\frac{\pi}{4} + \frac{\varphi}{2} \right) &= \frac{1 + \sin \varphi}{\cos \varphi} \\ &= \sqrt{\frac{(1 + \sin \varphi)^2}{\cos^2 \varphi}} \\ &= \sqrt{\frac{(1 + \sin \varphi)^2}{1 - \sin^2 \varphi}} \\ &= \sqrt{\frac{(1 + \sin \varphi)^2}{(1 + \sin \varphi)(1 - \sin \varphi)}} \\ &= \sqrt{\frac{1 + \sin \varphi}{1 - \sin \varphi}} \end{align}

なので、

\[y = \frac{128}{\pi} \left(\pi - \frac{1}{2} \ln \frac{1 + \sin \varphi}{1 - \sin \varphi} \right)\]

となる。

上で導いた投影式によって変換される座標系を、Google Mapsではworld coordinateと呼んでいる。 特定のズームレベルにおいて、地図画像上のピクセルを指定する座標系がpixel coordinateで、投影式は次のようになる。

\begin{align} x' &= 2^{zoom\ level} x \\ y' &= 2^{zoom\ level} y \end{align}

2017年12月26日火曜日

2点間の方位角を求める公式を導出する

地球を球体とみなす場合、球面三角法の公式を変形するだけで導くことができる。球面三角法の導出は地図における数学で解説されている。

2点を\(A\), \(B\)、極を\(C\)、地球の中心を\(O\)とし、頂角も同じ文字で表す。また、\(\angle AOB\), \(\angle BOC\), \(\angle COA\)をそれぞれ\(c\), \(a\), \(b\)とする。なお、地球の半径を\(1\)とすれば、大円の弧\(AB\), \(BC\), \(CA\)はそれぞれ\(c\), \(a\), \(b\)と等しい。

2点の緯度経度から以下の値がわかっていて、求めたいのは\(A\)の値。

\begin{align} a &= (90 - B_{latitude}) \frac{\pi}{180} \\ b &= (90 - A_{latitude}) \frac{\pi}{180} \\ C &= (B_{longitude} - A_{longitude}) \frac{\pi}{180} \end{align}

球面三角法の余弦定理を次のように変形する。

\begin{align} \cos a &= \cos b \cos c + \sin b \sin c \cos A \\ \sin b \sin c \cos A &= \cos a - \cos b \cos c \\ &= \cos a - \cos b (\cos a \cos b + \sin a \sin b \cos C) \\ &= \cos a - \cos a \cos^2 b - \sin a \sin b \cos b \cos C \\ &= \cos a (1 - \cos^2 b) - \sin a \sin b \cos b \cos C \\ &= \cos a \sin^2 b - \sin a \sin b \cos b \cos C \\ \sin c \cos A &= \cos a \sin b - \sin a \cos b \cos C \end{align}

正弦定理も変形する。

\begin{align} \frac{\sin C}{\sin c} &= \frac{\sin A}{\sin a} \\ \sin c \sin A &= \sin a \sin C \end{align}

上式から\(\tan A\)が求まり、\(A\)の値も求まる。

\begin{align} A &= \arctan(\tan A) \\ &= \arctan \left(\frac{\sin c \sin A}{\sin c \cos A} \right) \\ &= \arctan \left(\frac{\sin a \sin C}{\cos a \sin b - \sin a \cos b \cos C} \right) \end{align}

実際の計算では、先に\(\tan A\)を求めると第一象限と第三象限、第二象限と第四象限の区別がつかない。プログラムを書くときはatan2を使う。なお、\(A\)と\(B\)の経度が等しい場合でも正しい値が得られる。

また、

\begin{align} a' &= B_{latitude} \frac{\pi}{180} \\ b' &= A_{latitude} \frac{\pi}{180} \end{align}

と置くと、

\begin{align} A &= \arctan \left(\frac{\cos a' \sin C}{\sin a' \cos b' - \cos a' \sin b' \cos C} \right) \\ &= \arctan \left(\frac{\sin C}{\tan a' \cos b' - \sin b' \cos C} \right) \end{align}

とも書くことができる。

2017年3月29日水曜日

TBS News iの連続動画再生で新しいニュースのみ選択するブックマークレットとユーザースクリプト

News iの連続動画再生ビューアに何ヶ月も前のニュースが混じっていたりして、毎回新しい動画を選択するのが辛いので、ユーザースクリプトとついでにブックマークレットを書いた。テストはFirefox 52 ESRとGreasemonkeyで行ってある。

(2017/12/11 追記)スクリプトが複数回実行された時に正しく選択されないバグを修正。

ブックマークレット

以下のリンクをブックマークに登録して、連続動画再生ビューアのページで普通のブックマークと同じようにクリックすれば、12時間以内のニュースが選択されるはず。

(2017/3/29 追記)属性値のダブルクォートをエスケープし忘れていたのを修正した。

もっと新しいニュースのみ選択したい場合、例えば6時間以内なら、ブックマークのアドレス部分の末尾を(12)から(6)に変更すればいい。

ユーザースクリプト

// ==UserScript==
// @name        news i selector
// @include     http://news.tbs.co.jp/3snewsi/*
// @grant       none
// ==/UserScript==

var HOURS = 12;

setTimeout(function() {
    var expiredAt = Date.now() / 1000 - HOURS * 60 * 60;
    Array.prototype.forEach.call(document.querySelectorAll(".ls-checkItem"), function(item) {
        var checked = item.querySelector("input[type='checkbox']").checked;
        var time = parseInt(item.querySelector(".gs-metaData>span").getAttribute("data-time"), 10);
        var expired = (time < expiredAt);
        if (checked) {
            item.click();
        }
        if (!expired) {
            item.click();
        }
    });
    // setTimeout(function() {
    //     document.querySelector(".md-3snewsiPlayButton .ls-playSelected").click();
    // }, 0);
}, 3000);

(2017/9/20 追記)最近のアップデートで選択のタイミングが合わなくなってしまったので、setTimeoutを使って遅らせるように修正した。

(2017/12/14 追記)GreasemonkeyがWebExtensionになったことで動かなくなったのを修正。。

最後のコメントアウトをはずすと再生も自動でされる。しかし、連続動画再生ビューアはどういうわけか勝手にリロードすることがあり、その場合、また最初の動画から再生されてしまうので使い勝手は良くない。