

.. _sphx_glr_gallery_api_custom_projection_example.py:


=================
Custom projection
=================

This example showcases the Hammer projection by alleviating many features of
matplotlib.




.. image:: /gallery/api/images/sphx_glr_custom_projection_example_001.png
    :align: center





.. code-block:: python



    from __future__ import unicode_literals

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

    rcParams = matplotlib.rcParams

    # This example projection class is rather long, but it is designed to
    # illustrate many features, not all of which will be used every time.
    # It is also common to factor out a lot of these methods into common
    # code used by a number of projections with similar characteristics
    # (see geo.py).


    class GeoAxes(Axes):
        """
        An abstract base class for geographic projections
        """
        class ThetaFormatter(Formatter):
            """
            Used to format the theta tick labels.  Converts the native
            unit of radians into degrees and adds a degree symbol.
            """
            def __init__(self, round_to=1.0):
                self._round_to = round_to

            def __call__(self, x, pos=None):
                degrees = np.round(np.rad2deg(x) / self._round_to) * self._round_to
                if rcParams['text.usetex'] and not rcParams['text.latex.unicode']:
                    return r"$%0.0f^\circ$" % degrees
                else:
                    return "%0.0f\N{DEGREE SIGN}" % degrees

        RESOLUTION = 75

        def _init_axis(self):
            self.xaxis = maxis.XAxis(self)
            self.yaxis = maxis.YAxis(self)
            # Do not register xaxis or yaxis with spines -- as done in
            # Axes._init_axis() -- until GeoAxes.xaxis.cla() works.
            # self.spines['geo'].register_axis(self.yaxis)
            self._update_transScale()

        def cla(self):
            Axes.cla(self)

            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)
            # Why do we need to turn on yaxis tick labels, but
            # xaxis tick labels are already on?

            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)

        def _set_lim_and_transforms(self):
            # A (possibly non-linear) projection on the (already scaled) data

            # There are three important coordinate spaces going on here:
            #
            #    1. Data space: The space of the data itself
            #
            #    2. Axes space: The unit rectangle (0, 0) to (1, 1)
            #       covering the entire plot area.
            #
            #    3. Display space: The coordinates of the resulting image,
            #       often in pixels or dpi/inch.

            # This function makes heavy use of the Transform classes in
            # ``lib/matplotlib/transforms.py.`` For more information, see
            # the inline documentation there.

            # The goal of the first two transformations is to get from the
            # data space (in this case longitude and latitude) to axes
            # space.  It is separated into a non-affine and affine part so
            # that the non-affine part does not have to be recomputed when
            # a simple affine change to the figure has been made (such as
            # resizing the window or changing the dpi).

            # 1) The core transformation from data space into
            # rectilinear space defined in the HammerTransform class.
            self.transProjection = self._get_core_transform(self.RESOLUTION)

            # 2) The above has an output range that is not in the unit
            # rectangle, so scale and translate it so it fits correctly
            # within the axes.  The peculiar calculations of xscale and
            # yscale are specific to a Aitoff-Hammer projection, so don't
            # worry about them too much.
            self.transAffine = self._get_affine_transform()

            # 3) This is the transformation from axes space to display
            # space.
            self.transAxes = BboxTransformTo(self.bbox)

            # Now put these 3 transforms together -- from data all the way
            # to display coordinates.  Using the '+' operator, these
            # transforms will be applied "in order".  The transforms are
            # automatically simplified, if possible, by the underlying
            # transformation framework.
            self.transData = \
                self.transProjection + \
                self.transAffine + \
                self.transAxes

            # The main data transformation is set up.  Now deal with
            # gridlines and tick labels.

            # Longitude gridlines and ticklabels.  The input to these
            # transforms are in display space in x and axes space in y.
            # Therefore, the input values will be in range (-xmin, 0),
            # (xmax, 1).  The goal of these transforms is to go from that
            # space to display space.  The tick labels will be offset 4
            # pixels from the equator.
            self._xaxis_pretransform = \
                Affine2D() \
                .scale(1.0, self._longitude_cap * 2.0) \
                .translate(0.0, -self._longitude_cap)
            self._xaxis_transform = \
                self._xaxis_pretransform + \
                self.transData
            self._xaxis_text1_transform = \
                Affine2D().scale(1.0, 0.0) + \
                self.transData + \
                Affine2D().translate(0.0, 4.0)
            self._xaxis_text2_transform = \
                Affine2D().scale(1.0, 0.0) + \
                self.transData + \
                Affine2D().translate(0.0, -4.0)

            # Now set up the transforms for the latitude ticks.  The input to
            # these transforms are in axes space in x and display space in
            # y.  Therefore, the input values will be in range (0, -ymin),
            # (1, ymax).  The goal of these transforms is to go from that
            # space to display space.  The tick labels will be offset 4
            # pixels from the edge of the axes ellipse.
            yaxis_stretch = Affine2D().scale(np.pi*2, 1).translate(-np.pi, 0)
            yaxis_space = Affine2D().scale(1.0, 1.1)
            self._yaxis_transform = \
                yaxis_stretch + \
                self.transData
            yaxis_text_base = \
                yaxis_stretch + \
                self.transProjection + \
                (yaxis_space +
                 self.transAffine +
                 self.transAxes)
            self._yaxis_text1_transform = \
                yaxis_text_base + \
                Affine2D().translate(-8.0, 0.0)
            self._yaxis_text2_transform = \
                yaxis_text_base + \
                Affine2D().translate(8.0, 0.0)

        def _get_affine_transform(self):
            transform = self._get_core_transform(1)
            xscale, _ = transform.transform_point((np.pi, 0))
            _, yscale = transform.transform_point((0, np.pi / 2.0))
            return Affine2D() \
                .scale(0.5 / xscale, 0.5 / yscale) \
                .translate(0.5, 0.5)

        def get_xaxis_transform(self, which='grid'):
            """
            Override this method to provide a transformation for the
            x-axis tick labels.

            Returns a tuple of the form (transform, valign, halign)
            """
            if which not in ['tick1', 'tick2', 'grid']:
                msg = "'which' must be on of [ 'tick1' | 'tick2' | 'grid' ]"
                raise ValueError(msg)
            return self._xaxis_transform

        def get_xaxis_text1_transform(self, pad):
            return self._xaxis_text1_transform, 'bottom', 'center'

        def get_xaxis_text2_transform(self, pad):
            """
            Override this method to provide a transformation for the
            secondary x-axis tick labels.

            Returns a tuple of the form (transform, valign, halign)
            """
            return self._xaxis_text2_transform, 'top', 'center'

        def get_yaxis_transform(self, which='grid'):
            """
            Override this method to provide a transformation for the
            y-axis grid and ticks.
            """
            if which not in ['tick1', 'tick2', 'grid']:
                msg = "'which' must be one of [ 'tick1' | 'tick2' | 'grid' ]"
                raise ValueError(msg)
            return self._yaxis_transform

        def get_yaxis_text1_transform(self, pad):
            """
            Override this method to provide a transformation for the
            y-axis tick labels.

            Returns a tuple of the form (transform, valign, halign)
            """
            return self._yaxis_text1_transform, 'center', 'right'

        def get_yaxis_text2_transform(self, pad):
            """
            Override this method to provide a transformation for the
            secondary y-axis tick labels.

            Returns a tuple of the form (transform, valign, halign)
            """
            return self._yaxis_text2_transform, 'center', 'left'

        def _gen_axes_patch(self):
            """
            Override this method to define the shape that is used for the
            background of the plot.  It should be a subclass of Patch.

            In this case, it is a Circle (that may be warped by the axes
            transform into an ellipse).  Any data and gridlines will be
            clipped to this shape.
            """
            return Circle((0.5, 0.5), 0.5)

        def _gen_axes_spines(self):
            return {'geo': mspines.Spine.circular_spine(self, (0.5, 0.5), 0.5)}

        def set_yscale(self, *args, **kwargs):
            if args[0] != 'linear':
                raise NotImplementedError

        # Prevent the user from applying scales to one or both of the
        # axes.  In this particular case, scaling the axes wouldn't make
        # sense, so we don't allow it.
        set_xscale = set_yscale

        # Prevent the user from changing the axes limits.  In our case, we
        # want to display the whole sphere all the time, so we override
        # set_xlim and set_ylim to ignore any input.  This also applies to
        # interactive panning and zooming in the GUI interfaces.
        def set_xlim(self, *args, **kwargs):
            raise TypeError("It is not possible to change axes limits "
                            "for geographic projections. Please consider "
                            "using Basemap or Cartopy.")

        set_ylim = set_xlim

        def format_coord(self, lon, lat):
            """
            Override this method to change how the values are displayed in
            the status bar.

            In this case, we want them to be displayed in degrees N/S/E/W.
            """
            lon, lat = np.rad2deg([lon, lat])
            if lat >= 0.0:
                ns = 'N'
            else:
                ns = 'S'
            if lon >= 0.0:
                ew = 'E'
            else:
                ew = 'W'
            return ('%f\N{DEGREE SIGN}%s, %f\N{DEGREE SIGN}%s'
                    % (abs(lat), ns, abs(lon), ew))

        def set_longitude_grid(self, degrees):
            """
            Set the number of degrees between each longitude grid.

            This is an example method that is specific to this projection
            class -- it provides a more convenient interface to set the
            ticking than set_xticks would.
            """
            # Skip -180 and 180, which are the fixed limits.
            grid = np.arange(-180 + degrees, 180, degrees)
            self.xaxis.set_major_locator(FixedLocator(np.deg2rad(grid)))
            self.xaxis.set_major_formatter(self.ThetaFormatter(degrees))

        def set_latitude_grid(self, degrees):
            """
            Set the number of degrees between each longitude grid.

            This is an example method that is specific to this projection
            class -- it provides a more convenient interface than
            set_yticks would.
            """
            # Skip -90 and 90, which are the fixed limits.
            grid = np.arange(-90 + degrees, 90, degrees)
            self.yaxis.set_major_locator(FixedLocator(np.deg2rad(grid)))
            self.yaxis.set_major_formatter(self.ThetaFormatter(degrees))

        def set_longitude_grid_ends(self, degrees):
            """
            Set the latitude(s) at which to stop drawing the longitude grids.

            Often, in geographic projections, you wouldn't want to draw
            longitude gridlines near the poles.  This allows the user to
            specify the degree at which to stop drawing longitude grids.

            This is an example method that is specific to this projection
            class -- it provides an interface to something that has no
            analogy in the base Axes class.
            """
            self._longitude_cap = np.deg2rad(degrees)
            self._xaxis_pretransform \
                .clear() \
                .scale(1.0, self._longitude_cap * 2.0) \
                .translate(0.0, -self._longitude_cap)

        def get_data_ratio(self):
            """
            Return the aspect ratio of the data itself.

            This method should be overridden by any Axes that have a
            fixed data ratio.
            """
            return 1.0

        # Interactive panning and zooming is not supported with this projection,
        # so we override all of the following methods to disable it.
        def can_zoom(self):
            """
            Return *True* if this axes supports the zoom box button functionality.
            This axes object does not support interactive zoom box.
            """
            return False

        def can_pan(self):
            """
            Return *True* if this axes supports the pan/zoom button functionality.
            This axes object does not support interactive pan/zoom.
            """
            return False

        def start_pan(self, x, y, button):
            pass

        def end_pan(self):
            pass

        def drag_pan(self, button, key, x, y):
            pass


    class HammerAxes(GeoAxes):
        """
        A custom class for the Aitoff-Hammer projection, an equal-area map
        projection.

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

        # The projection must specify a name. This will be used by the
        # user to select the projection,
        # i.e. ``subplot(111, projection='custom_hammer')``.
        name = 'custom_hammer'

        class HammerTransform(Transform):
            """
            The base Hammer transform.
            """
            input_dims = 2
            output_dims = 2
            is_separable = False

            def __init__(self, resolution):
                """
                Create a new Hammer transform.  Resolution is the number of steps
                to interpolate between each input line segment to approximate its
                path in curved Hammer space.
                """
                Transform.__init__(self)
                self._resolution = resolution

            def transform_non_affine(self, ll):
                longitude = ll[:, 0:1]
                latitude = ll[:, 1:2]

                # Pre-compute some values
                half_long = longitude / 2.0
                cos_latitude = np.cos(latitude)
                sqrt2 = np.sqrt(2.0)

                alpha = np.sqrt(1.0 + cos_latitude * np.cos(half_long))
                x = (2.0 * sqrt2) * (cos_latitude * np.sin(half_long)) / alpha
                y = (sqrt2 * np.sin(latitude)) / alpha
                return np.concatenate((x, y), 1)
            transform_non_affine.__doc__ = Transform.transform_non_affine.__doc__

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

            def inverted(self):
                return HammerAxes.InvertedHammerTransform(self._resolution)
            inverted.__doc__ = Transform.inverted.__doc__

        class InvertedHammerTransform(Transform):
            input_dims = 2
            output_dims = 2
            is_separable = False

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

            def transform_non_affine(self, xy):
                x = xy[:, 0:1]
                y = xy[:, 1:2]

                quarter_x = 0.25 * x
                half_y = 0.5 * y
                z = np.sqrt(1.0 - quarter_x*quarter_x - half_y*half_y)
                longitude = 2 * np.arctan((z*x) / (2.0 * (2.0*z*z - 1.0)))
                latitude = np.arcsin(y*z)
                return np.concatenate((longitude, latitude), 1)
            transform_non_affine.__doc__ = Transform.transform_non_affine.__doc__

            def inverted(self):
                return HammerAxes.HammerTransform(self._resolution)
            inverted.__doc__ = Transform.inverted.__doc__

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

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


    # Now register the projection with matplotlib so the user can select
    # it.
    register_projection(HammerAxes)


    if __name__ == '__main__':
        import matplotlib.pyplot as plt
        # Now make a simple example using the custom projection.
        plt.subplot(111, projection="custom_hammer")
        p = plt.plot([-1, 1, 1], [-1, -1, 1], "o-")
        plt.grid(True)

        plt.show()

**Total running time of the script:** ( 0 minutes  0.017 seconds)



.. container:: sphx-glr-footer


  .. container:: sphx-glr-download

     :download:`Download Python source code: custom_projection_example.py <custom_projection_example.py>`



  .. container:: sphx-glr-download

     :download:`Download Jupyter notebook: custom_projection_example.ipynb <custom_projection_example.ipynb>`

.. rst-class:: sphx-glr-signature

    `Generated by Sphinx-Gallery <https://sphinx-gallery.readthedocs.io>`_
