`pole_density_function()` plots do not appear correct in the center

[yencal]

The center of plots produced by pole_density_function() do not appear to be correct. For example https://orix.readthedocs.io/en/stable/examples/plotting/subplots.html#sphx-glr-examples-plotting-subplots-py. The center of the plot appear distorted and I have similar issues with my pole density plots also.

[hakonanes]

Hi @yencal,

I believe they are correct, as @harripj who added that functionality intended. He recently suggested that we add an option to plot the densities as contours (https://github.com/pyxem/orix/discussions/440#discussioncomment-5470704) instead of the current histogram bins. I'm not sure when that will be added, though...

[yencal]

Hi @hakonanes, that would be great as most people are used to seeing the plots as contours. Thanks

[harripj]

@yencal thanks for opening an issue about this. To properly measure the pole distribution on the unit sphere we use equal area sampling on S2 rather than equal angle sampling. This is to ensure that the probability that a given pole intersects a bin on the S2 sphere is constant. As you point out, this means that the bins near the poles have a larger latitudinal dimension when using equal area sampling, which appears to distort pole region of the distribution.

I've written a short example highlighting the two sampling methods, and the resulting distribution for a set of random vectors (which should then give a uniform distribution):

from orix.vector import Vector3d
from orix import plot
from orix.sampling.S2_sampling import _sample_S2_uv_mesh_coordinates, _sample_S2_equal_area_coordinates
from matplotlib import pyplot as plt
import numpy as np
from orix.projections import StereographicProjection

v = Vector3d(np.random.randn(1_000_000, 3))
azimuth, polar, _ = v.to_polar()

resolution = 2 # deg.
sp = StereographicProjection()

# equal angle grid
ac, pc = _sample_S2_uv_mesh_coordinates(resolution, hemisphere='upper', azimuth_endpoint=True)
hist_uv, *_ = np.histogram2d(
    azimuth,
    polar,
    bins=(ac, pc),
    density=False,
)
uv_grid = Vector3d.from_polar(*np.meshgrid(ac, pc, indexing='ij'))
x_uv, y_uv = sp.vector2xy(uv_grid)

# equal area grid
ac, pc = _sample_S2_equal_area_coordinates(resolution, hemisphere='upper', azimuth_endpoint=True)
hist_ea, *_ = np.histogram2d(
    azimuth,
    polar,
    bins=(ac, pc),
    density=False,
)
ea_grid = Vector3d.from_polar(*np.meshgrid(ac, pc, indexing='ij'))
x_ea, y_ea = sp.vector2xy(ea_grid)

fig, ax = plt.subplots(nrows=2, ncols=2, figsize=(10, 10), subplot_kw=dict(projection='stereographic'))

ax[0, 0].scatter(uv_grid, s=0.5)
ax[0, 1].scatter(ea_grid, s=0.5)
ax[0, 0].set_title('Equal angle grid points')
ax[0, 1].set_title('Equal area grid points')

pc_uv = ax[1, 0].pcolorfast(x_uv.reshape(uv_grid.shape), y_uv.reshape(uv_grid.shape), hist_uv / hist_uv.mean())
pc_ea = ax[1, 1].pcolorfast(x_ea.reshape(ea_grid.shape), y_ea.reshape(ea_grid.shape), hist_ea / hist_ea.mean())
plt.colorbar(pc_uv, ax=ax[1, 0], label='MRD')
plt.colorbar(pc_ea, ax=ax[1, 1], label='MRD')
ax[1, 0].set_title('Distribution on equal angle grid')
ax[1, 1].set_title('Distribution on equal area grid')

From the two lower plots it's clear that the distribution is uniform as expected when using equal area sampling, but is underrepresented at the poles when using equal angle sampling for the reason discussed above. This is discussed in further detail in the paper by Rohrer (RSD+04 in the orix bibliography).

On the topic of contour plots this would definitely be a nice addition to orix, and please feel free to contribute if you already have an implementation written!

[yencal]

Hi @harripj, thanks for the explanation and example code.