grid_helper_curvelinear.py 12 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336
  1. """
  2. An experimental support for curvilinear grid.
  3. """
  4. import functools
  5. from itertools import chain
  6. import numpy as np
  7. import matplotlib as mpl
  8. from matplotlib.path import Path
  9. from matplotlib.transforms import Affine2D, IdentityTransform
  10. from .axislines import (
  11. _FixedAxisArtistHelperBase, _FloatingAxisArtistHelperBase, GridHelperBase)
  12. from .axis_artist import AxisArtist
  13. from .grid_finder import GridFinder
  14. def _value_and_jacobian(func, xs, ys, xlims, ylims):
  15. """
  16. Compute *func* and its derivatives along x and y at positions *xs*, *ys*,
  17. while ensuring that finite difference calculations don't try to evaluate
  18. values outside of *xlims*, *ylims*.
  19. """
  20. eps = np.finfo(float).eps ** (1/2) # see e.g. scipy.optimize.approx_fprime
  21. val = func(xs, ys)
  22. # Take the finite difference step in the direction where the bound is the
  23. # furthest; the step size is min of epsilon and distance to that bound.
  24. xlo, xhi = sorted(xlims)
  25. dxlo = xs - xlo
  26. dxhi = xhi - xs
  27. xeps = (np.take([-1, 1], dxhi >= dxlo)
  28. * np.minimum(eps, np.maximum(dxlo, dxhi)))
  29. val_dx = func(xs + xeps, ys)
  30. ylo, yhi = sorted(ylims)
  31. dylo = ys - ylo
  32. dyhi = yhi - ys
  33. yeps = (np.take([-1, 1], dyhi >= dylo)
  34. * np.minimum(eps, np.maximum(dylo, dyhi)))
  35. val_dy = func(xs, ys + yeps)
  36. return (val, (val_dx - val) / xeps, (val_dy - val) / yeps)
  37. class FixedAxisArtistHelper(_FixedAxisArtistHelperBase):
  38. """
  39. Helper class for a fixed axis.
  40. """
  41. def __init__(self, grid_helper, side, nth_coord_ticks=None):
  42. """
  43. nth_coord = along which coordinate value varies.
  44. nth_coord = 0 -> x axis, nth_coord = 1 -> y axis
  45. """
  46. super().__init__(loc=side)
  47. self.grid_helper = grid_helper
  48. if nth_coord_ticks is None:
  49. nth_coord_ticks = self.nth_coord
  50. self.nth_coord_ticks = nth_coord_ticks
  51. self.side = side
  52. def update_lim(self, axes):
  53. self.grid_helper.update_lim(axes)
  54. def get_tick_transform(self, axes):
  55. return axes.transData
  56. def get_tick_iterators(self, axes):
  57. """tick_loc, tick_angle, tick_label"""
  58. v1, v2 = axes.get_ylim() if self.nth_coord == 0 else axes.get_xlim()
  59. if v1 > v2: # Inverted limits.
  60. side = {"left": "right", "right": "left",
  61. "top": "bottom", "bottom": "top"}[self.side]
  62. else:
  63. side = self.side
  64. g = self.grid_helper
  65. ti1 = g.get_tick_iterator(self.nth_coord_ticks, side)
  66. ti2 = g.get_tick_iterator(1-self.nth_coord_ticks, side, minor=True)
  67. return chain(ti1, ti2), iter([])
  68. class FloatingAxisArtistHelper(_FloatingAxisArtistHelperBase):
  69. def __init__(self, grid_helper, nth_coord, value, axis_direction=None):
  70. """
  71. nth_coord = along which coordinate value varies.
  72. nth_coord = 0 -> x axis, nth_coord = 1 -> y axis
  73. """
  74. super().__init__(nth_coord, value)
  75. self.value = value
  76. self.grid_helper = grid_helper
  77. self._extremes = -np.inf, np.inf
  78. self._line_num_points = 100 # number of points to create a line
  79. def set_extremes(self, e1, e2):
  80. if e1 is None:
  81. e1 = -np.inf
  82. if e2 is None:
  83. e2 = np.inf
  84. self._extremes = e1, e2
  85. def update_lim(self, axes):
  86. self.grid_helper.update_lim(axes)
  87. x1, x2 = axes.get_xlim()
  88. y1, y2 = axes.get_ylim()
  89. grid_finder = self.grid_helper.grid_finder
  90. extremes = grid_finder.extreme_finder(grid_finder.inv_transform_xy,
  91. x1, y1, x2, y2)
  92. lon_min, lon_max, lat_min, lat_max = extremes
  93. e_min, e_max = self._extremes # ranges of other coordinates
  94. if self.nth_coord == 0:
  95. lat_min = max(e_min, lat_min)
  96. lat_max = min(e_max, lat_max)
  97. elif self.nth_coord == 1:
  98. lon_min = max(e_min, lon_min)
  99. lon_max = min(e_max, lon_max)
  100. lon_levs, lon_n, lon_factor = \
  101. grid_finder.grid_locator1(lon_min, lon_max)
  102. lat_levs, lat_n, lat_factor = \
  103. grid_finder.grid_locator2(lat_min, lat_max)
  104. if self.nth_coord == 0:
  105. xx0 = np.full(self._line_num_points, self.value)
  106. yy0 = np.linspace(lat_min, lat_max, self._line_num_points)
  107. xx, yy = grid_finder.transform_xy(xx0, yy0)
  108. elif self.nth_coord == 1:
  109. xx0 = np.linspace(lon_min, lon_max, self._line_num_points)
  110. yy0 = np.full(self._line_num_points, self.value)
  111. xx, yy = grid_finder.transform_xy(xx0, yy0)
  112. self._grid_info = {
  113. "extremes": (lon_min, lon_max, lat_min, lat_max),
  114. "lon_info": (lon_levs, lon_n, np.asarray(lon_factor)),
  115. "lat_info": (lat_levs, lat_n, np.asarray(lat_factor)),
  116. "lon_labels": grid_finder.tick_formatter1(
  117. "bottom", lon_factor, lon_levs),
  118. "lat_labels": grid_finder.tick_formatter2(
  119. "bottom", lat_factor, lat_levs),
  120. "line_xy": (xx, yy),
  121. }
  122. def get_axislabel_transform(self, axes):
  123. return Affine2D() # axes.transData
  124. def get_axislabel_pos_angle(self, axes):
  125. def trf_xy(x, y):
  126. trf = self.grid_helper.grid_finder.get_transform() + axes.transData
  127. return trf.transform([x, y]).T
  128. xmin, xmax, ymin, ymax = self._grid_info["extremes"]
  129. if self.nth_coord == 0:
  130. xx0 = self.value
  131. yy0 = (ymin + ymax) / 2
  132. elif self.nth_coord == 1:
  133. xx0 = (xmin + xmax) / 2
  134. yy0 = self.value
  135. xy1, dxy1_dx, dxy1_dy = _value_and_jacobian(
  136. trf_xy, xx0, yy0, (xmin, xmax), (ymin, ymax))
  137. p = axes.transAxes.inverted().transform(xy1)
  138. if 0 <= p[0] <= 1 and 0 <= p[1] <= 1:
  139. d = [dxy1_dy, dxy1_dx][self.nth_coord]
  140. return xy1, np.rad2deg(np.arctan2(*d[::-1]))
  141. else:
  142. return None, None
  143. def get_tick_transform(self, axes):
  144. return IdentityTransform() # axes.transData
  145. def get_tick_iterators(self, axes):
  146. """tick_loc, tick_angle, tick_label, (optionally) tick_label"""
  147. lat_levs, lat_n, lat_factor = self._grid_info["lat_info"]
  148. yy0 = lat_levs / lat_factor
  149. lon_levs, lon_n, lon_factor = self._grid_info["lon_info"]
  150. xx0 = lon_levs / lon_factor
  151. e0, e1 = self._extremes
  152. def trf_xy(x, y):
  153. trf = self.grid_helper.grid_finder.get_transform() + axes.transData
  154. return trf.transform(np.column_stack(np.broadcast_arrays(x, y))).T
  155. # find angles
  156. if self.nth_coord == 0:
  157. mask = (e0 <= yy0) & (yy0 <= e1)
  158. (xx1, yy1), (dxx1, dyy1), (dxx2, dyy2) = _value_and_jacobian(
  159. trf_xy, self.value, yy0[mask], (-np.inf, np.inf), (e0, e1))
  160. labels = self._grid_info["lat_labels"]
  161. elif self.nth_coord == 1:
  162. mask = (e0 <= xx0) & (xx0 <= e1)
  163. (xx1, yy1), (dxx2, dyy2), (dxx1, dyy1) = _value_and_jacobian(
  164. trf_xy, xx0[mask], self.value, (-np.inf, np.inf), (e0, e1))
  165. labels = self._grid_info["lon_labels"]
  166. labels = [l for l, m in zip(labels, mask) if m]
  167. angle_normal = np.arctan2(dyy1, dxx1)
  168. angle_tangent = np.arctan2(dyy2, dxx2)
  169. mm = (dyy1 == 0) & (dxx1 == 0) # points with degenerate normal
  170. angle_normal[mm] = angle_tangent[mm] + np.pi / 2
  171. tick_to_axes = self.get_tick_transform(axes) - axes.transAxes
  172. in_01 = functools.partial(
  173. mpl.transforms._interval_contains_close, (0, 1))
  174. def f1():
  175. for x, y, normal, tangent, lab \
  176. in zip(xx1, yy1, angle_normal, angle_tangent, labels):
  177. c2 = tick_to_axes.transform((x, y))
  178. if in_01(c2[0]) and in_01(c2[1]):
  179. yield [x, y], *np.rad2deg([normal, tangent]), lab
  180. return f1(), iter([])
  181. def get_line_transform(self, axes):
  182. return axes.transData
  183. def get_line(self, axes):
  184. self.update_lim(axes)
  185. x, y = self._grid_info["line_xy"]
  186. return Path(np.column_stack([x, y]))
  187. class GridHelperCurveLinear(GridHelperBase):
  188. def __init__(self, aux_trans,
  189. extreme_finder=None,
  190. grid_locator1=None,
  191. grid_locator2=None,
  192. tick_formatter1=None,
  193. tick_formatter2=None):
  194. """
  195. Parameters
  196. ----------
  197. aux_trans : `.Transform` or tuple[Callable, Callable]
  198. The transform from curved coordinates to rectilinear coordinate:
  199. either a `.Transform` instance (which provides also its inverse),
  200. or a pair of callables ``(trans, inv_trans)`` that define the
  201. transform and its inverse. The callables should have signature::
  202. x_rect, y_rect = trans(x_curved, y_curved)
  203. x_curved, y_curved = inv_trans(x_rect, y_rect)
  204. extreme_finder
  205. grid_locator1, grid_locator2
  206. Grid locators for each axis.
  207. tick_formatter1, tick_formatter2
  208. Tick formatters for each axis.
  209. """
  210. super().__init__()
  211. self._grid_info = None
  212. self.grid_finder = GridFinder(aux_trans,
  213. extreme_finder,
  214. grid_locator1,
  215. grid_locator2,
  216. tick_formatter1,
  217. tick_formatter2)
  218. def update_grid_finder(self, aux_trans=None, **kwargs):
  219. if aux_trans is not None:
  220. self.grid_finder.update_transform(aux_trans)
  221. self.grid_finder.update(**kwargs)
  222. self._old_limits = None # Force revalidation.
  223. def new_fixed_axis(self, loc,
  224. nth_coord=None,
  225. axis_direction=None,
  226. offset=None,
  227. axes=None):
  228. if axes is None:
  229. axes = self.axes
  230. if axis_direction is None:
  231. axis_direction = loc
  232. helper = FixedAxisArtistHelper(self, loc, nth_coord_ticks=nth_coord)
  233. axisline = AxisArtist(axes, helper, axis_direction=axis_direction)
  234. # Why is clip not set on axisline, unlike in new_floating_axis or in
  235. # the floating_axig.GridHelperCurveLinear subclass?
  236. return axisline
  237. def new_floating_axis(self, nth_coord,
  238. value,
  239. axes=None,
  240. axis_direction="bottom"
  241. ):
  242. if axes is None:
  243. axes = self.axes
  244. helper = FloatingAxisArtistHelper(
  245. self, nth_coord, value, axis_direction)
  246. axisline = AxisArtist(axes, helper)
  247. axisline.line.set_clip_on(True)
  248. axisline.line.set_clip_box(axisline.axes.bbox)
  249. # axisline.major_ticklabels.set_visible(True)
  250. # axisline.minor_ticklabels.set_visible(False)
  251. return axisline
  252. def _update_grid(self, x1, y1, x2, y2):
  253. self._grid_info = self.grid_finder.get_grid_info(x1, y1, x2, y2)
  254. def get_gridlines(self, which="major", axis="both"):
  255. grid_lines = []
  256. if axis in ["both", "x"]:
  257. for gl in self._grid_info["lon"]["lines"]:
  258. grid_lines.extend(gl)
  259. if axis in ["both", "y"]:
  260. for gl in self._grid_info["lat"]["lines"]:
  261. grid_lines.extend(gl)
  262. return grid_lines
  263. def get_tick_iterator(self, nth_coord, axis_side, minor=False):
  264. # axisnr = dict(left=0, bottom=1, right=2, top=3)[axis_side]
  265. angle_tangent = dict(left=90, right=90, bottom=0, top=0)[axis_side]
  266. # angle = [0, 90, 180, 270][axisnr]
  267. lon_or_lat = ["lon", "lat"][nth_coord]
  268. if not minor: # major ticks
  269. for (xy, a), l in zip(
  270. self._grid_info[lon_or_lat]["tick_locs"][axis_side],
  271. self._grid_info[lon_or_lat]["tick_labels"][axis_side]):
  272. angle_normal = a
  273. yield xy, angle_normal, angle_tangent, l
  274. else:
  275. for (xy, a), l in zip(
  276. self._grid_info[lon_or_lat]["tick_locs"][axis_side],
  277. self._grid_info[lon_or_lat]["tick_labels"][axis_side]):
  278. angle_normal = a
  279. yield xy, angle_normal, angle_tangent, ""
  280. # for xy, a, l in self._grid_info[lon_or_lat]["ticks"][axis_side]:
  281. # yield xy, a, ""