angle_helper.py 13 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394
  1. import numpy as np
  2. import math
  3. from mpl_toolkits.axisartist.grid_finder import ExtremeFinderSimple
  4. def select_step_degree(dv):
  5. degree_limits_ = [1.5, 3, 7, 13, 20, 40, 70, 120, 270, 520]
  6. degree_steps_ = [1, 2, 5, 10, 15, 30, 45, 90, 180, 360]
  7. degree_factors = [1.] * len(degree_steps_)
  8. minsec_limits_ = [1.5, 2.5, 3.5, 8, 11, 18, 25, 45]
  9. minsec_steps_ = [1, 2, 3, 5, 10, 15, 20, 30]
  10. minute_limits_ = np.array(minsec_limits_) / 60
  11. minute_factors = [60.] * len(minute_limits_)
  12. second_limits_ = np.array(minsec_limits_) / 3600
  13. second_factors = [3600.] * len(second_limits_)
  14. degree_limits = [*second_limits_, *minute_limits_, *degree_limits_]
  15. degree_steps = [*minsec_steps_, *minsec_steps_, *degree_steps_]
  16. degree_factors = [*second_factors, *minute_factors, *degree_factors]
  17. n = np.searchsorted(degree_limits, dv)
  18. step = degree_steps[n]
  19. factor = degree_factors[n]
  20. return step, factor
  21. def select_step_hour(dv):
  22. hour_limits_ = [1.5, 2.5, 3.5, 5, 7, 10, 15, 21, 36]
  23. hour_steps_ = [1, 2, 3, 4, 6, 8, 12, 18, 24]
  24. hour_factors = [1.] * len(hour_steps_)
  25. minsec_limits_ = [1.5, 2.5, 3.5, 4.5, 5.5, 8, 11, 14, 18, 25, 45]
  26. minsec_steps_ = [1, 2, 3, 4, 5, 6, 10, 12, 15, 20, 30]
  27. minute_limits_ = np.array(minsec_limits_) / 60
  28. minute_factors = [60.] * len(minute_limits_)
  29. second_limits_ = np.array(minsec_limits_) / 3600
  30. second_factors = [3600.] * len(second_limits_)
  31. hour_limits = [*second_limits_, *minute_limits_, *hour_limits_]
  32. hour_steps = [*minsec_steps_, *minsec_steps_, *hour_steps_]
  33. hour_factors = [*second_factors, *minute_factors, *hour_factors]
  34. n = np.searchsorted(hour_limits, dv)
  35. step = hour_steps[n]
  36. factor = hour_factors[n]
  37. return step, factor
  38. def select_step_sub(dv):
  39. # subarcsec or degree
  40. tmp = 10.**(int(math.log10(dv))-1.)
  41. factor = 1./tmp
  42. if 1.5*tmp >= dv:
  43. step = 1
  44. elif 3.*tmp >= dv:
  45. step = 2
  46. elif 7.*tmp >= dv:
  47. step = 5
  48. else:
  49. step = 1
  50. factor = 0.1*factor
  51. return step, factor
  52. def select_step(v1, v2, nv, hour=False, include_last=True,
  53. threshold_factor=3600.):
  54. if v1 > v2:
  55. v1, v2 = v2, v1
  56. dv = (v2 - v1) / nv
  57. if hour:
  58. _select_step = select_step_hour
  59. cycle = 24.
  60. else:
  61. _select_step = select_step_degree
  62. cycle = 360.
  63. # for degree
  64. if dv > 1 / threshold_factor:
  65. step, factor = _select_step(dv)
  66. else:
  67. step, factor = select_step_sub(dv*threshold_factor)
  68. factor = factor * threshold_factor
  69. levs = np.arange(np.floor(v1 * factor / step),
  70. np.ceil(v2 * factor / step) + 0.5,
  71. dtype=int) * step
  72. # n : number of valid levels. If there is a cycle, e.g., [0, 90, 180,
  73. # 270, 360], the grid line needs to be extended from 0 to 360, so
  74. # we need to return the whole array. However, the last level (360)
  75. # needs to be ignored often. In this case, so we return n=4.
  76. n = len(levs)
  77. # we need to check the range of values
  78. # for example, -90 to 90, 0 to 360,
  79. if factor == 1. and levs[-1] >= levs[0] + cycle: # check for cycle
  80. nv = int(cycle / step)
  81. if include_last:
  82. levs = levs[0] + np.arange(0, nv+1, 1) * step
  83. else:
  84. levs = levs[0] + np.arange(0, nv, 1) * step
  85. n = len(levs)
  86. return np.array(levs), n, factor
  87. def select_step24(v1, v2, nv, include_last=True, threshold_factor=3600):
  88. v1, v2 = v1 / 15, v2 / 15
  89. levs, n, factor = select_step(v1, v2, nv, hour=True,
  90. include_last=include_last,
  91. threshold_factor=threshold_factor)
  92. return levs * 15, n, factor
  93. def select_step360(v1, v2, nv, include_last=True, threshold_factor=3600):
  94. return select_step(v1, v2, nv, hour=False,
  95. include_last=include_last,
  96. threshold_factor=threshold_factor)
  97. class LocatorBase:
  98. def __init__(self, nbins, include_last=True):
  99. self.nbins = nbins
  100. self._include_last = include_last
  101. def set_params(self, nbins=None):
  102. if nbins is not None:
  103. self.nbins = int(nbins)
  104. class LocatorHMS(LocatorBase):
  105. def __call__(self, v1, v2):
  106. return select_step24(v1, v2, self.nbins, self._include_last)
  107. class LocatorHM(LocatorBase):
  108. def __call__(self, v1, v2):
  109. return select_step24(v1, v2, self.nbins, self._include_last,
  110. threshold_factor=60)
  111. class LocatorH(LocatorBase):
  112. def __call__(self, v1, v2):
  113. return select_step24(v1, v2, self.nbins, self._include_last,
  114. threshold_factor=1)
  115. class LocatorDMS(LocatorBase):
  116. def __call__(self, v1, v2):
  117. return select_step360(v1, v2, self.nbins, self._include_last)
  118. class LocatorDM(LocatorBase):
  119. def __call__(self, v1, v2):
  120. return select_step360(v1, v2, self.nbins, self._include_last,
  121. threshold_factor=60)
  122. class LocatorD(LocatorBase):
  123. def __call__(self, v1, v2):
  124. return select_step360(v1, v2, self.nbins, self._include_last,
  125. threshold_factor=1)
  126. class FormatterDMS:
  127. deg_mark = r"^{\circ}"
  128. min_mark = r"^{\prime}"
  129. sec_mark = r"^{\prime\prime}"
  130. fmt_d = "$%d" + deg_mark + "$"
  131. fmt_ds = r"$%d.%s" + deg_mark + "$"
  132. # %s for sign
  133. fmt_d_m = r"$%s%d" + deg_mark + r"\,%02d" + min_mark + "$"
  134. fmt_d_ms = r"$%s%d" + deg_mark + r"\,%02d.%s" + min_mark + "$"
  135. fmt_d_m_partial = "$%s%d" + deg_mark + r"\,%02d" + min_mark + r"\,"
  136. fmt_s_partial = "%02d" + sec_mark + "$"
  137. fmt_ss_partial = "%02d.%s" + sec_mark + "$"
  138. def _get_number_fraction(self, factor):
  139. ## check for fractional numbers
  140. number_fraction = None
  141. # check for 60
  142. for threshold in [1, 60, 3600]:
  143. if factor <= threshold:
  144. break
  145. d = factor // threshold
  146. int_log_d = int(np.floor(np.log10(d)))
  147. if 10**int_log_d == d and d != 1:
  148. number_fraction = int_log_d
  149. factor = factor // 10**int_log_d
  150. return factor, number_fraction
  151. return factor, number_fraction
  152. def __call__(self, direction, factor, values):
  153. if len(values) == 0:
  154. return []
  155. ss = np.sign(values)
  156. signs = ["-" if v < 0 else "" for v in values]
  157. factor, number_fraction = self._get_number_fraction(factor)
  158. values = np.abs(values)
  159. if number_fraction is not None:
  160. values, frac_part = divmod(values, 10 ** number_fraction)
  161. frac_fmt = "%%0%dd" % (number_fraction,)
  162. frac_str = [frac_fmt % (f1,) for f1 in frac_part]
  163. if factor == 1:
  164. if number_fraction is None:
  165. return [self.fmt_d % (s * int(v),) for s, v in zip(ss, values)]
  166. else:
  167. return [self.fmt_ds % (s * int(v), f1)
  168. for s, v, f1 in zip(ss, values, frac_str)]
  169. elif factor == 60:
  170. deg_part, min_part = divmod(values, 60)
  171. if number_fraction is None:
  172. return [self.fmt_d_m % (s1, d1, m1)
  173. for s1, d1, m1 in zip(signs, deg_part, min_part)]
  174. else:
  175. return [self.fmt_d_ms % (s, d1, m1, f1)
  176. for s, d1, m1, f1
  177. in zip(signs, deg_part, min_part, frac_str)]
  178. elif factor == 3600:
  179. if ss[-1] == -1:
  180. inverse_order = True
  181. values = values[::-1]
  182. signs = signs[::-1]
  183. else:
  184. inverse_order = False
  185. l_hm_old = ""
  186. r = []
  187. deg_part, min_part_ = divmod(values, 3600)
  188. min_part, sec_part = divmod(min_part_, 60)
  189. if number_fraction is None:
  190. sec_str = [self.fmt_s_partial % (s1,) for s1 in sec_part]
  191. else:
  192. sec_str = [self.fmt_ss_partial % (s1, f1)
  193. for s1, f1 in zip(sec_part, frac_str)]
  194. for s, d1, m1, s1 in zip(signs, deg_part, min_part, sec_str):
  195. l_hm = self.fmt_d_m_partial % (s, d1, m1)
  196. if l_hm != l_hm_old:
  197. l_hm_old = l_hm
  198. l = l_hm + s1
  199. else:
  200. l = "$" + s + s1
  201. r.append(l)
  202. if inverse_order:
  203. return r[::-1]
  204. else:
  205. return r
  206. else: # factor > 3600.
  207. return [r"$%s^{\circ}$" % v for v in ss*values]
  208. class FormatterHMS(FormatterDMS):
  209. deg_mark = r"^\mathrm{h}"
  210. min_mark = r"^\mathrm{m}"
  211. sec_mark = r"^\mathrm{s}"
  212. fmt_d = "$%d" + deg_mark + "$"
  213. fmt_ds = r"$%d.%s" + deg_mark + "$"
  214. # %s for sign
  215. fmt_d_m = r"$%s%d" + deg_mark + r"\,%02d" + min_mark+"$"
  216. fmt_d_ms = r"$%s%d" + deg_mark + r"\,%02d.%s" + min_mark+"$"
  217. fmt_d_m_partial = "$%s%d" + deg_mark + r"\,%02d" + min_mark + r"\,"
  218. fmt_s_partial = "%02d" + sec_mark + "$"
  219. fmt_ss_partial = "%02d.%s" + sec_mark + "$"
  220. def __call__(self, direction, factor, values): # hour
  221. return super().__call__(direction, factor, np.asarray(values) / 15)
  222. class ExtremeFinderCycle(ExtremeFinderSimple):
  223. # docstring inherited
  224. def __init__(self, nx, ny,
  225. lon_cycle=360., lat_cycle=None,
  226. lon_minmax=None, lat_minmax=(-90, 90)):
  227. """
  228. This subclass handles the case where one or both coordinates should be
  229. taken modulo 360, or be restricted to not exceed a specific range.
  230. Parameters
  231. ----------
  232. nx, ny : int
  233. The number of samples in each direction.
  234. lon_cycle, lat_cycle : 360 or None
  235. If not None, values in the corresponding direction are taken modulo
  236. *lon_cycle* or *lat_cycle*; in theory this can be any number but
  237. the implementation actually assumes that it is 360 (if not None);
  238. other values give nonsensical results.
  239. This is done by "unwrapping" the transformed grid coordinates so
  240. that jumps are less than a half-cycle; then normalizing the span to
  241. no more than a full cycle.
  242. For example, if values are in the union of the [0, 2] and
  243. [358, 360] intervals (typically, angles measured modulo 360), the
  244. values in the second interval are normalized to [-2, 0] instead so
  245. that the values now cover [-2, 2]. If values are in a range of
  246. [5, 1000], this gets normalized to [5, 365].
  247. lon_minmax, lat_minmax : (float, float) or None
  248. If not None, the computed bounding box is clipped to the given
  249. range in the corresponding direction.
  250. """
  251. self.nx, self.ny = nx, ny
  252. self.lon_cycle, self.lat_cycle = lon_cycle, lat_cycle
  253. self.lon_minmax = lon_minmax
  254. self.lat_minmax = lat_minmax
  255. def __call__(self, transform_xy, x1, y1, x2, y2):
  256. # docstring inherited
  257. x, y = np.meshgrid(
  258. np.linspace(x1, x2, self.nx), np.linspace(y1, y2, self.ny))
  259. lon, lat = transform_xy(np.ravel(x), np.ravel(y))
  260. # iron out jumps, but algorithm should be improved.
  261. # This is just naive way of doing and my fail for some cases.
  262. # Consider replacing this with numpy.unwrap
  263. # We are ignoring invalid warnings. They are triggered when
  264. # comparing arrays with NaNs using > We are already handling
  265. # that correctly using np.nanmin and np.nanmax
  266. with np.errstate(invalid='ignore'):
  267. if self.lon_cycle is not None:
  268. lon0 = np.nanmin(lon)
  269. lon -= 360. * ((lon - lon0) > 180.)
  270. if self.lat_cycle is not None:
  271. lat0 = np.nanmin(lat)
  272. lat -= 360. * ((lat - lat0) > 180.)
  273. lon_min, lon_max = np.nanmin(lon), np.nanmax(lon)
  274. lat_min, lat_max = np.nanmin(lat), np.nanmax(lat)
  275. lon_min, lon_max, lat_min, lat_max = \
  276. self._add_pad(lon_min, lon_max, lat_min, lat_max)
  277. # check cycle
  278. if self.lon_cycle:
  279. lon_max = min(lon_max, lon_min + self.lon_cycle)
  280. if self.lat_cycle:
  281. lat_max = min(lat_max, lat_min + self.lat_cycle)
  282. if self.lon_minmax is not None:
  283. min0 = self.lon_minmax[0]
  284. lon_min = max(min0, lon_min)
  285. max0 = self.lon_minmax[1]
  286. lon_max = min(max0, lon_max)
  287. if self.lat_minmax is not None:
  288. min0 = self.lat_minmax[0]
  289. lat_min = max(min0, lat_min)
  290. max0 = self.lat_minmax[1]
  291. lat_max = min(max0, lat_max)
  292. return lon_min, lon_max, lat_min, lat_max