Skip to content

Navigation Menu

Sign in
Appearance settings

Search code, repositories, users, issues, pull requests...

Provide feedback

We read every piece of feedback, and take your input very seriously.

Saved searches

Use saved searches to filter your results more quickly
Appearance settings

Commit 5ca9137

Browse filesBrowse files
anntzeranntzer
committed
Improve handling of degenerate jacobians in non-rectilinear grids.
grid_helper_curvelinear and floating_axes have code to specifically handle the case where the transform from rectlinear to non-rectilinear axes has null derivatives in one of the directions, inferring the angle of the jacobian from the derivative in the other direction. (This angle defines the rotation applied to axis labels, ticks, and tick labels.) This approach, however, is insufficient if the derivatives in both directions are zero. A classical example is e.g. the ``exp(-1/x**2)`` transform, for which all derivatives are zero. To handle this case more robustly (and also to better encapsulate the angle calculation, which is currently repeated at a few places), instead, one can increase the step size of the numerical differentiation until the gradient becomes nonzero. This amounts to moving along the corresponding gridline until one actually leaves the position of the tick, and thus is indeed a justifiable approach to compute the tick rotation. Full examples: import numpy as np from matplotlib import pyplot as plt from matplotlib.projections.polar import PolarTransform from matplotlib.transforms import Affine2D from mpl_toolkits.axisartist import ( angle_helper, GridHelperCurveLinear, HostAxes) import mpl_toolkits.axisartist.floating_axes as floating_axes # def tr(x, y): return x - y, x + y # def inv_tr(u, v): return (u + v) / 2, (v - u) / 2 @np.errstate(divide="ignore") # at x=0, exp(-1/x**2)=0; div-by-zero can be ignored. def tr(x, y): return np.exp(-x**-2) - np.exp(-y**-2), np.exp(-x**-2) + np.exp(-y**-2) def inv_tr(u, v): return (-np.log((u+v)/2))**(1/2), (-np.log((v-u)/2))**(1/2) plt.subplot( 121, axes_class=floating_axes.FloatingAxes, grid_helper=floating_axes.GridHelperCurveLinear( (tr, inv_tr), extremes=(0, 10, 0, 10))) ax = plt.subplot( 122, axes_class=HostAxes, grid_helper=GridHelperCurveLinear( Affine2D().scale(np.pi / 180, 1) + PolarTransform() + Affine2D().scale(2, 1), extreme_finder=angle_helper.ExtremeFinderCycle( 20, 20, lon_cycle=360, lat_cycle=None, lon_minmax=None, lat_minmax=(0, np.inf), ), grid_locator1=angle_helper.LocatorDMS(12), tick_formatter1=angle_helper.FormatterDMS(), ), aspect=1, xlim=(-5, 12), ylim=(-5, 10)) ax.axis["lat"] = axis = ax.new_floating_axis(0, 40) axis.label.set_text(r"$\theta = 40^{\circ}$") axis.label.set_visible(True) ax.grid(True) plt.show()
1 parent 07002c2 commit 5ca9137
Copy full SHA for 5ca9137

File tree

Expand file treeCollapse file tree

3 files changed

+102
-41
lines changed
Filter options
Expand file treeCollapse file tree

3 files changed

+102
-41
lines changed

‎lib/mpl_toolkits/axisartist/floating_axes.py

Copy file name to clipboardExpand all lines: lib/mpl_toolkits/axisartist/floating_axes.py
+4-10Lines changed: 4 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -71,25 +71,19 @@ def trf_xy(x, y):
7171

7272
if self.nth_coord == 0:
7373
mask = (ymin <= yy0) & (yy0 <= ymax)
74-
(xx1, yy1), (dxx1, dyy1), (dxx2, dyy2) = \
75-
grid_helper_curvelinear._value_and_jacobian(
74+
(xx1, yy1), angle_normal, angle_tangent = \
75+
grid_helper_curvelinear._value_and_jac_angle(
7676
trf_xy, self.value, yy0[mask], (xmin, xmax), (ymin, ymax))
7777
labels = self._grid_info["lat_labels"]
7878

7979
elif self.nth_coord == 1:
8080
mask = (xmin <= xx0) & (xx0 <= xmax)
81-
(xx1, yy1), (dxx2, dyy2), (dxx1, dyy1) = \
82-
grid_helper_curvelinear._value_and_jacobian(
81+
(xx1, yy1), angle_tangent, angle_normal = \
82+
grid_helper_curvelinear._value_and_jac_angle(
8383
trf_xy, xx0[mask], self.value, (xmin, xmax), (ymin, ymax))
8484
labels = self._grid_info["lon_labels"]
8585

8686
labels = [l for l, m in zip(labels, mask) if m]
87-
88-
angle_normal = np.arctan2(dyy1, dxx1)
89-
angle_tangent = np.arctan2(dyy2, dxx2)
90-
mm = (dyy1 == 0) & (dxx1 == 0) # points with degenerate normal
91-
angle_normal[mm] = angle_tangent[mm] + np.pi / 2
92-
9387
tick_to_axes = self.get_tick_transform(axes) - axes.transAxes
9488
in_01 = functools.partial(
9589
mpl.transforms._interval_contains_close, (0, 1))

‎lib/mpl_toolkits/axisartist/grid_helper_curvelinear.py

Copy file name to clipboardExpand all lines: lib/mpl_toolkits/axisartist/grid_helper_curvelinear.py
+70-31Lines changed: 70 additions & 31 deletions
Original file line numberDiff line numberDiff line change
@@ -16,29 +16,75 @@
1616
from .grid_finder import GridFinder
1717

1818

19-
def _value_and_jacobian(func, xs, ys, xlims, ylims):
19+
def _value_and_jac_angle(func, xs, ys, xlim, ylim):
2020
"""
21-
Compute *func* and its derivatives along x and y at positions *xs*, *ys*,
22-
while ensuring that finite difference calculations don't try to evaluate
23-
values outside of *xlims*, *ylims*.
21+
Parameters
22+
----------
23+
func : callable
24+
A function that transforms the coordinates of a point (x, y) to a new coordinate
25+
system (u, v), and which can also take x and y as arrays of shape *shape* and
26+
returns (u, v) as a ``(2, shape)`` array.
27+
xs, ys : array-likes
28+
Points where *func* and its derivatives will be evaluated.
29+
xlim, ylim : pairs of floats
30+
(min, max) beyond which *func* should not be evaluated.
31+
32+
Returns
33+
-------
34+
val
35+
Value of *func* at each point of ``(xs, ys)``.
36+
thetas_dx
37+
Angles (in radians) defined by the (u, v) components of the numerically
38+
differentiated df/dx vector, at each point of ``(xs, ys)``. If needed, the
39+
differentiation step size is increased until at least one component of df/dx
40+
is nonzero, under the constraint of not going out of the *xlims*, *ylims*
41+
bounds. If the gridline at a point is actually null (and the angle is thus not
42+
well defined), the derivatives are evaluated after taking a small step along y;
43+
this ensures e.g. that the tick at r=0 on a radial axis of a polar plot is
44+
parallel with the ticks at r!=0.
45+
thetas_dy
46+
Like *thetas_dx*, but for df/dy.
2447
"""
25-
eps = np.finfo(float).eps ** (1/2) # see e.g. scipy.optimize.approx_fprime
48+
49+
shape = np.broadcast_shapes(np.shape(xs), np.shape(ys))
2650
val = func(xs, ys)
27-
# Take the finite difference step in the direction where the bound is the
28-
# furthest; the step size is min of epsilon and distance to that bound.
29-
xlo, xhi = sorted(xlims)
30-
dxlo = xs - xlo
31-
dxhi = xhi - xs
32-
xeps = (np.take([-1, 1], dxhi >= dxlo)
33-
* np.minimum(eps, np.maximum(dxlo, dxhi)))
34-
val_dx = func(xs + xeps, ys)
35-
ylo, yhi = sorted(ylims)
36-
dylo = ys - ylo
37-
dyhi = yhi - ys
38-
yeps = (np.take([-1, 1], dyhi >= dylo)
39-
* np.minimum(eps, np.maximum(dylo, dyhi)))
40-
val_dy = func(xs, ys + yeps)
41-
return (val, (val_dx - val) / xeps, (val_dy - val) / yeps)
51+
52+
# Take finite difference steps towards the furthest bound; the step size will be the
53+
# min of epsilon and the distance to that bound.
54+
eps0 = np.finfo(float).eps ** (1/2) # cf. scipy.optimize.approx_fprime
55+
56+
def calc_eps(vals, lim):
57+
lo, hi = sorted(lim)
58+
dlo = vals - lo
59+
dhi = hi - vals
60+
eps_max = np.maximum(dlo, dhi)
61+
eps = np.where(dhi >= dlo, 1, -1) * np.minimum(eps0, eps_max)
62+
return eps, eps_max
63+
64+
xeps, xeps_max = calc_eps(xs, xlim)
65+
yeps, yeps_max = calc_eps(ys, ylim)
66+
67+
def calc_thetas(dfunc, ps, eps_p0, eps_max, eps_q):
68+
thetas_dp = np.full(shape, np.nan)
69+
missing = np.full(shape, True)
70+
eps_p = eps_p0
71+
for it, eps_q in enumerate([0, eps_q]):
72+
while missing.any() and (abs(eps_p) < eps_max).any():
73+
if it == 0 and (eps_p > 1).any():
74+
break # Degenerate derivative, move a bit along the other coord.
75+
eps_p = np.minimum(eps_p, eps_max)
76+
df_x, df_y = (dfunc(eps_p, eps_q) - dfunc(0, eps_q)) / eps_p
77+
good = missing & ((df_x != 0) | (df_y != 0))
78+
thetas_dp[good] = np.arctan2(df_y, df_x)[good]
79+
missing &= ~good
80+
eps_p *= 2
81+
return thetas_dp
82+
83+
thetas_dx = calc_thetas(lambda eps_p, eps_q: func(xs + eps_p, ys + eps_q),
84+
xs, xeps, xeps_max, yeps)
85+
thetas_dy = calc_thetas(lambda eps_p, eps_q: func(xs + eps_q, ys + eps_p),
86+
ys, yeps, yeps_max, xeps)
87+
return (val, thetas_dx, thetas_dy)
4288

4389

4490
class FixedAxisArtistHelper(_FixedAxisArtistHelperBase):
@@ -167,12 +213,11 @@ def trf_xy(x, y):
167213
elif self.nth_coord == 1:
168214
xx0 = (xmin + xmax) / 2
169215
yy0 = self.value
170-
xy1, dxy1_dx, dxy1_dy = _value_and_jacobian(
216+
xy1, angle_dx, angle_dy = _value_and_jac_angle(
171217
trf_xy, xx0, yy0, (xmin, xmax), (ymin, ymax))
172218
p = axes.transAxes.inverted().transform(xy1)
173219
if 0 <= p[0] <= 1 and 0 <= p[1] <= 1:
174-
d = [dxy1_dy, dxy1_dx][self.nth_coord]
175-
return xy1, np.rad2deg(np.arctan2(*d[::-1]))
220+
return xy1, np.rad2deg([angle_dy, angle_dx][self.nth_coord])
176221
else:
177222
return None, None
178223

@@ -197,23 +242,17 @@ def trf_xy(x, y):
197242
# find angles
198243
if self.nth_coord == 0:
199244
mask = (e0 <= yy0) & (yy0 <= e1)
200-
(xx1, yy1), (dxx1, dyy1), (dxx2, dyy2) = _value_and_jacobian(
245+
(xx1, yy1), angle_normal, angle_tangent = _value_and_jac_angle(
201246
trf_xy, self.value, yy0[mask], (-np.inf, np.inf), (e0, e1))
202247
labels = self._grid_info["lat_labels"]
203248

204249
elif self.nth_coord == 1:
205250
mask = (e0 <= xx0) & (xx0 <= e1)
206-
(xx1, yy1), (dxx2, dyy2), (dxx1, dyy1) = _value_and_jacobian(
251+
(xx1, yy1), angle_tangent, angle_normal = _value_and_jac_angle(
207252
trf_xy, xx0[mask], self.value, (-np.inf, np.inf), (e0, e1))
208253
labels = self._grid_info["lon_labels"]
209254

210255
labels = [l for l, m in zip(labels, mask) if m]
211-
212-
angle_normal = np.arctan2(dyy1, dxx1)
213-
angle_tangent = np.arctan2(dyy2, dxx2)
214-
mm = (dyy1 == 0) & (dxx1 == 0) # points with degenerate normal
215-
angle_normal[mm] = angle_tangent[mm] + np.pi / 2
216-
217256
tick_to_axes = self.get_tick_transform(axes) - axes.transAxes
218257
in_01 = functools.partial(
219258
mpl.transforms._interval_contains_close, (0, 1))

‎lib/mpl_toolkits/axisartist/tests/test_floating_axes.py

Copy file name to clipboardExpand all lines: lib/mpl_toolkits/axisartist/tests/test_floating_axes.py
+28Lines changed: 28 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,7 @@
11
import numpy as np
22

3+
import pytest
4+
35
import matplotlib.pyplot as plt
46
import matplotlib.projections as mprojections
57
import matplotlib.transforms as mtransforms
@@ -113,3 +115,29 @@ def test_axis_direction():
113115
ax.axis['y'] = ax.new_floating_axis(nth_coord=1, value=0,
114116
axis_direction='left')
115117
assert ax.axis['y']._axis_direction == 'left'
118+
119+
120+
def test_transform_with_zero_derivatives():
121+
# The transform is really a 45° rotation
122+
# tr(x, y) = x-y, x+y; inv_tr(u, v) = (u+v)/2, (u-v)/2
123+
# with an additional x->exp(-x**-2) on each coordinate.
124+
# Therefore all ticks should be at +/-45°, even the one at zero where the
125+
# transform derivatives are zero.
126+
127+
# at x=0, exp(-x**-2)=0; div-by-zero can be ignored.
128+
@np.errstate(divide="ignore")
129+
def tr(x, y):
130+
return np.exp(-x**-2) - np.exp(-y**-2), np.exp(-x**-2) + np.exp(-y**-2)
131+
132+
def inv_tr(u, v):
133+
return (-np.log((u+v)/2))**(1/2), (-np.log((v-u)/2))**(1/2)
134+
135+
fig = plt.figure()
136+
ax = fig.add_subplot(
137+
axes_class=FloatingAxes, grid_helper=GridHelperCurveLinear(
138+
(tr, inv_tr), extremes=(0, 10, 0, 10)))
139+
fig.canvas.draw()
140+
141+
for k in ax.axis:
142+
for l, a in ax.axis[k].major_ticks.locs_angles:
143+
assert a % 90 == pytest.approx(45, abs=1e-3)

0 commit comments

Comments
0 (0)
Morty Proxy This is a proxified and sanitized view of the page, visit original site.