地図投影法を開発してみたシリーズ4

Equal-area ovalシリーズ3部作その2。

https://research.tableau.com/paper/bevy-area-preserving-transforms-map-projection-designersから

Directional path offset、前回のように擬円筒図法を詰めてつなぎ合わせるだけじゃなく、経緯度関係なく詰めたりしてもいいんだなー思いつかなかった。例えばこんな計算ということか。

f:id:dummy_index:20200831233553p:plain — Offsetted azimuthal equal-area double hemisphere

エイトフ変換をせずに席詰め（以降directional path offsetのことをこう呼ぶことにします）をすると緯線が90度を超えて反ってしまうな（水平から90度を超えるかどうかより地図の端で経線と交差する角度のほうが力学的に気になる）エイトフ変換をしても極付近の緯線はほぼ垂直になるから他になにか挟んだだけでまた90度超えしそうだし、もっとエイトフ変換について考えないと…

http://www.cartography.oregonstate.edu/pdf/2014_Savric_Jenny_ANewPseudocylindricalEqual-areaProjection.pdfの3ページ

$1 \lt B \lt 2$ は緯線が両方の意味で反りが強すぎるのでそのままでは採用できません。 $B=4$ (Eckert-Greifendorff projection) は緯線を水平に近づけるのには成功しているけどそれが同時にHammerから高緯度の歪で知られるSinusoidalに近づけるような真似になってしまいます。 $B \lt 2$ にするのではない方法で逆方向に変化させるには…楕円からひし形に近づけるのの逆をするには…

こうじゃ！

f:id:dummy_index:20200901233522p:plain — Eckert-Greifendorff projection uses -pi/4 ~ +pi/4

f:id:dummy_index:20200901233643p:plain — Use the portion out of Eckert-Greifendorff projection, then do the directional path offset

f:id:dummy_index:20200901234243p:plain — Complementary Eckert-Greifendorff projection

Y座標からoffset量を求めるのに最初はニュートン法を使っていたのですが、Wolfram|Alphaを頼って閉形式の式を得ました。ヒントになったのがhttps://research.tableau.com/paper/bevy-area-preserving-transforms-map-projection-designersのHamusoidalの $y_c$ と $\phi_{CM}$ の関係式。これなら $\lambda=0$ 以外でも明示的に解けそう…

solve y = (2/(1+cos(phi)*c))^0.5 * sin(phi) for phi - Wolfram|Alpha

そうしたら前提の $\lambda$ と得られた $\phi$ から改めてforward formulaで $x$ を得ます*1。式を整理して

$x = \frac{\sqrt{(\cos^2(\lambda)-2)y^2 - \cos(\lambda)\cdot \sqrt{\cos^2(\lambda)\cdot y^4-8y^2+16}+4}}{\sqrt 2} \cdot \mbox{sgn}(\lambda) \cdot \mbox{sgn}(\cos^2(\lambda)\cdot y^4-8y^2+16)$