Matplotlib を使ったカスタム投影 | Python プログラミング

はじめに

この実験では、Matplotlibを使ってカスタム投影を作成する方法を学びます。Matplotlibの多くの機能を活かして、ハンマー投影を紹介します。プログラミング言語としてPythonを使用します。

VMのヒント

VMの起動が完了したら、左上隅をクリックしてノートブックタブに切り替え、Jupyter Notebookを使って練習しましょう。

Jupyter Notebookが読み込み終わるまで数秒待つことがあります。Jupyter Notebookの制限により、操作の検証を自動化することはできません。

学習中に問題がある場合は、Labbyにお問い合わせください。セッション後にフィードバックを提供してください。すぐに問題を解決いたします。

ライブラリのインポート

まず、カスタム投影を作成するために必要なライブラリをインポートします。

import numpy as np
import matplotlib
from matplotlib.axes import Axes
import matplotlib.axis as maxis
from matplotlib.patches import Circle
from matplotlib.path import Path
from matplotlib.projections import register_projection
import matplotlib.spines as mspines
from matplotlib.ticker import FixedLocator, Formatter, NullLocator
from matplotlib.transforms import Affine2D, BboxTransformTo, Transform

GeoAxes クラスの作成

地理投影用の抽象基底クラスである GeoAxes を作成します。

class GeoAxes(Axes):
    """
    地理投影用の抽象基底クラス
    """

    class ThetaFormatter(Formatter):
        """
        ゼータ目盛のラベルをフォーマットするために使用されます。
        ラジアンの元の単位を度に変換し、度の記号を追加します。
        """
        def __init__(self, round_to=1.0):
            self._round_to = round_to

        def __call__(self, x, pos=None):
            degrees = round(np.rad2deg(x) / self._round_to) * self._round_to
            return f"{degrees:0.0f}\N{DEGREE SIGN}"

    RESOLUTION = 75

    def _init_axis(self):
        self.xaxis = maxis.XAxis(self)
        self.yaxis = maxis.YAxis(self)
        ## GeoAxes.xaxis.clear() が機能するまで、Axes._init_axis() で行われているように、
        ## xaxis または yaxis をスパインに登録しないでください。
        ## self.spines['geo'].register_axis(self.yaxis)

    def clear(self):
        ## ドキュメント文字列は継承されます
        super().clear()

        self.set_longitude_grid(30)
        self.set_latitude_grid(15)
        self.set_longitude_grid_ends(75)
        self.xaxis.set_minor_locator(NullLocator())
        self.yaxis.set_minor_locator(NullLocator())
        self.xaxis.set_ticks_position('none')
        self.yaxis.set_ticks_position('none')
        self.yaxis.set_tick_params(label1On=True)
        ## なぜ yaxis の目盛ラベルをオンにする必要があるのでしょうか。
        ## 一方、xaxis の目盛ラベルは既にオンになっています。

        self.grid(rcParams['axes.grid'])

        Axes.set_xlim(self, -np.pi, np.pi)
        Axes.set_ylim(self, -np.pi / 2.0, np.pi / 2.0)

注：「ゼータ」は原文の「theta」の訳語候補の一つですが、この文脈では「ゼータ」が正しい訳語かどうかは疑問です。もし「theta」が本来の角度を表す用語であることが確認できれば、「ゼータ」の代わりに「ゼータ（theta）」のように表記するか、正しい角度名に置き換える必要があります。また、「rcParams」はそのままのまま残しておくことが一般的ですが、もしこれが何か特定の設定を表すものであることが分かっている場合、それに合わせた適切な日本語訳を考える必要があります。ここではそのまま残しています。

ハンマー座標軸（HammerAxes）クラスの作成

等面積地図投影であるアイトフ・ハンマー投影用のカスタムクラスである「HammerAxes」を作成します。

class HammerAxes(GeoAxes):
    """
    等面積地図投影であるアイトフ・ハンマー投影用のカスタムクラス。

    https://en.wikipedia.org/wiki/Hammer_projection
    """

    ## 投影には名前を指定する必要があります。これはユーザーが投影を選択する際に使用されます。
    ## すなわち、``subplot(projection='custom_hammer')`` のように。
    name = 'custom_hammer'

    class HammerTransform(Transform):
        """基本的なハンマー変換。"""
        input_dims = output_dims = 2

        def __init__(self, resolution):
            """
            新しいハンマー変換を作成します。解像度は、曲線状のハンマー空間における各入力線分の間を補間するためのステップ数です。
            """
            Transform.__init__(self)
            self._resolution = resolution

        def transform_non_affine(self, ll):
            longitude, latitude = ll.T

            ## いくつかの値を事前計算します
            half_long = longitude / 2
            cos_latitude = np.cos(latitude)
            sqrt2 = np.sqrt(2)

            alpha = np.sqrt(1 + cos_latitude * np.cos(half_long))
            x = (2 * sqrt2) * (cos_latitude * np.sin(half_long)) / alpha
            y = (sqrt2 * np.sin(latitude)) / alpha
            return np.column_stack([x, y])

        def transform_path_non_affine(self, path):
            ## vertices = path.vertices
            ipath = path.interpolated(self._resolution)
            return Path(self.transform(ipath.vertices), ipath.codes)

        def inverted(self):
            return HammerAxes.InvertedHammerTransform(self._resolution)

    class InvertedHammerTransform(Transform):
        input_dims = output_dims = 2

        def __init__(self, resolution):
            Transform.__init__(self)
            self._resolution = resolution

        def transform_non_affine(self, xy):
            x, y = xy.T
            z = np.sqrt(1 - (x / 4) ** 2 - (y / 2) ** 2)
            longitude = 2 * np.arctan((z * x) / (2 * (2 * z ** 2 - 1)))
            latitude = np.arcsin(y*z)
            return np.column_stack([longitude, latitude])

        def inverted(self):
            return HammerAxes.HammerTransform(self._resolution)

    def __init__(self, *args, **kwargs):
        self._longitude_cap = np.pi / 2.0
        super().__init__(*args, **kwargs)
        self.set_aspect(0.5, adjustable='box', anchor='C')
        self.clear()

    def _get_core_transform(self, resolution):
        return self.HammerTransform(resolution)

投影の登録

これで、Matplotlibに投影を登録して、ユーザーが選択できるようにします。

register_projection(HammerAxes)

サンプルの作成

最後に、カスタム投影を使用したサンプルを作成します。

if __name__ == '__main__':
    import matplotlib.pyplot as plt

    ## ここでは、カスタム投影を使用した簡単なサンプルを作成します。
    fig, ax = plt.subplots(subplot_kw={'projection': 'custom_hammer'})
    ax.plot([-1, 1, 1], [-1, -1, 1], "o-")
    ax.grid()

    plt.show()

まとめ

この実験では、Matplotlibを使用してカスタム投影を作成する方法を学びました。Matplotlibの多くの機能を利用して、ハンマー投影と呼ばれるカスタム投影を作成しました。また、GeoAxesクラス、HammerAxesクラスを作成する方法、投影を登録する方法、およびカスタム投影を使用したサンプルを作成する方法を学びました。

カスタム Matplotlib 投影の紹介