Introduction
In this lab, we will learn how to create a custom projection using Matplotlib. We will be showcasing Hammer projection by alleviating many features of Matplotlib. We will be using Python as our programming language.
VM Tips
After the VM startup is done, click the top left corner to switch to the Notebook tab to access Jupyter Notebook for practice.
Sometimes, you may need to wait a few seconds for Jupyter Notebook to finish loading. The validation of operations cannot be automated because of limitations in Jupyter Notebook.
If you face issues during learning, feel free to ask Labby. Provide feedback after the session, and we will promptly resolve the problem for you.
Import Libraries
First, we will import the necessary libraries to create a custom projection.
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
Create GeoAxes Class
We will create an abstract base class for geographic projections called GeoAxes.
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 = 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)
## Do not register xaxis or yaxis with spines -- as done in
## Axes._init_axis() -- until GeoAxes.xaxis.clear() works.
## self.spines['geo'].register_axis(self.yaxis)
def clear(self):
## docstring inherited
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)
## 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)
Create HammerAxes Class
We will create a custom class for the Aitoff-Hammer projection, an equal-area map projection called HammerAxes.
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(projection='custom_hammer')``.
name = 'custom_hammer'
class HammerTransform(Transform):
"""The base Hammer transform."""
input_dims = output_dims = 2
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, latitude = ll.T
## Pre-compute some values
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)
Register Projection
Now we will register the projection with Matplotlib so the user can select it.
register_projection(HammerAxes)
Create Example
Finally, we will create an example using the custom projection.
if __name__ == '__main__':
import matplotlib.pyplot as plt
## Now make a simple example using the custom projection.
fig, ax = plt.subplots(subplot_kw={'projection': 'custom_hammer'})
ax.plot([-1, 1, 1], [-1, -1, 1], "o-")
ax.grid()
plt.show()
Summary
In this lab, we learned how to create a custom projection using Matplotlib. We created a custom projection called Hammer projection by alleviating many features of Matplotlib. We learned how to create GeoAxes class, HammerAxes class, register the projection, and create an example using the custom projection.