affinity.py 8.1 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266
  1. """Affine transforms, both in general and specific, named transforms."""
  2. from math import cos, pi, sin, tan
  3. import numpy as np
  4. import shapely
  5. __all__ = ["affine_transform", "rotate", "scale", "skew", "translate"]
  6. def affine_transform(geom, matrix):
  7. r"""Return a transformed geometry using an affine transformation matrix.
  8. The coefficient matrix is provided as a list or tuple with 6 or 12 items
  9. for 2D or 3D transformations, respectively.
  10. For 2D affine transformations, the 6 parameter matrix is::
  11. [a, b, d, e, xoff, yoff]
  12. which represents the augmented matrix::
  13. [x'] / a b xoff \ [x]
  14. [y'] = | d e yoff | [y]
  15. [1 ] \ 0 0 1 / [1]
  16. or the equations for the transformed coordinates::
  17. x' = a * x + b * y + xoff
  18. y' = d * x + e * y + yoff
  19. For 3D affine transformations, the 12 parameter matrix is::
  20. [a, b, c, d, e, f, g, h, i, xoff, yoff, zoff]
  21. which represents the augmented matrix::
  22. [x'] / a b c xoff \ [x]
  23. [y'] = | d e f yoff | [y]
  24. [z'] | g h i zoff | [z]
  25. [1 ] \ 0 0 0 1 / [1]
  26. or the equations for the transformed coordinates::
  27. x' = a * x + b * y + c * z + xoff
  28. y' = d * x + e * y + f * z + yoff
  29. z' = g * x + h * y + i * z + zoff
  30. """
  31. if len(matrix) == 6:
  32. ndim = 2
  33. a, b, d, e, xoff, yoff = matrix
  34. if geom.has_z:
  35. ndim = 3
  36. i = 1.0
  37. c = f = g = h = zoff = 0.0
  38. elif len(matrix) == 12:
  39. ndim = 3
  40. a, b, c, d, e, f, g, h, i, xoff, yoff, zoff = matrix
  41. if not geom.has_z:
  42. ndim = 2
  43. else:
  44. raise ValueError("'matrix' expects either 6 or 12 coefficients")
  45. # if ndim == 2:
  46. # A = np.array([[a, b], [d, e]], dtype=float)
  47. # off = np.array([xoff, yoff], dtype=float)
  48. # else:
  49. # A = np.array([[a, b, c], [d, e, f], [g, h, i]], dtype=float)
  50. # off = np.array([xoff, yoff, zoff], dtype=float)
  51. def _affine_coords(coords):
  52. # These are equivalent, but unfortunately not robust
  53. # result = np.matmul(coords, A.T) + off
  54. # result = np.matmul(A, coords.T).T + off
  55. # Therefore, manual matrix multiplication is needed
  56. if ndim == 2:
  57. x, y = coords.T
  58. xp = a * x + b * y + xoff
  59. yp = d * x + e * y + yoff
  60. result = np.stack([xp, yp]).T
  61. elif ndim == 3:
  62. x, y, z = coords.T
  63. xp = a * x + b * y + c * z + xoff
  64. yp = d * x + e * y + f * z + yoff
  65. zp = g * x + h * y + i * z + zoff
  66. result = np.stack([xp, yp, zp]).T
  67. return result
  68. return shapely.transform(geom, _affine_coords, include_z=ndim == 3)
  69. def interpret_origin(geom, origin, ndim):
  70. """Return interpreted coordinate tuple for origin parameter.
  71. This is a helper function for other transform functions.
  72. The point of origin can be a keyword 'center' for the 2D bounding box
  73. center, 'centroid' for the geometry's 2D centroid, a Point object or a
  74. coordinate tuple (x0, y0, z0).
  75. """
  76. # get coordinate tuple from 'origin' from keyword or Point type
  77. if origin == "center":
  78. # bounding box center
  79. minx, miny, maxx, maxy = geom.bounds
  80. origin = ((maxx + minx) / 2.0, (maxy + miny) / 2.0)
  81. elif origin == "centroid":
  82. origin = geom.centroid.coords[0]
  83. elif isinstance(origin, str):
  84. raise ValueError(f"'origin' keyword {origin!r} is not recognized")
  85. elif getattr(origin, "geom_type", None) == "Point":
  86. origin = origin.coords[0]
  87. # origin should now be tuple-like
  88. if len(origin) not in (2, 3):
  89. raise ValueError("Expected number of items in 'origin' to be either 2 or 3")
  90. if ndim == 2:
  91. return origin[0:2]
  92. else: # 3D coordinate
  93. if len(origin) == 2:
  94. return origin + (0.0,)
  95. else:
  96. return origin
  97. def rotate(geom, angle, origin="center", use_radians=False):
  98. r"""Return a rotated geometry on a 2D plane.
  99. The angle of rotation can be specified in either degrees (default) or
  100. radians by setting ``use_radians=True``. Positive angles are
  101. counter-clockwise and negative are clockwise rotations.
  102. The point of origin can be a keyword 'center' for the bounding box
  103. center (default), 'centroid' for the geometry's centroid, a Point object
  104. or a coordinate tuple (x0, y0).
  105. The affine transformation matrix for 2D rotation is:
  106. / cos(r) -sin(r) xoff \
  107. | sin(r) cos(r) yoff |
  108. \ 0 0 1 /
  109. where the offsets are calculated from the origin Point(x0, y0):
  110. xoff = x0 - x0 * cos(r) + y0 * sin(r)
  111. yoff = y0 - x0 * sin(r) - y0 * cos(r)
  112. """
  113. if geom.is_empty:
  114. return geom
  115. if not use_radians: # convert from degrees
  116. angle = angle * pi / 180.0
  117. cosp = cos(angle)
  118. sinp = sin(angle)
  119. if abs(cosp) < 2.5e-16:
  120. cosp = 0.0
  121. if abs(sinp) < 2.5e-16:
  122. sinp = 0.0
  123. x0, y0 = interpret_origin(geom, origin, 2)
  124. # fmt: off
  125. matrix = (cosp, -sinp, 0.0,
  126. sinp, cosp, 0.0,
  127. 0.0, 0.0, 1.0,
  128. x0 - x0 * cosp + y0 * sinp, y0 - x0 * sinp - y0 * cosp, 0.0)
  129. # fmt: on
  130. return affine_transform(geom, matrix)
  131. def scale(geom, xfact=1.0, yfact=1.0, zfact=1.0, origin="center"):
  132. r"""Return a scaled geometry, scaled by factors along each dimension.
  133. The point of origin can be a keyword 'center' for the 2D bounding box
  134. center (default), 'centroid' for the geometry's 2D centroid, a Point
  135. object or a coordinate tuple (x0, y0, z0).
  136. Negative scale factors will mirror or reflect coordinates.
  137. The general 3D affine transformation matrix for scaling is:
  138. / xfact 0 0 xoff \
  139. | 0 yfact 0 yoff |
  140. | 0 0 zfact zoff |
  141. \ 0 0 0 1 /
  142. where the offsets are calculated from the origin Point(x0, y0, z0):
  143. xoff = x0 - x0 * xfact
  144. yoff = y0 - y0 * yfact
  145. zoff = z0 - z0 * zfact
  146. """
  147. if geom.is_empty:
  148. return geom
  149. x0, y0, z0 = interpret_origin(geom, origin, 3)
  150. # fmt: off
  151. matrix = (xfact, 0.0, 0.0,
  152. 0.0, yfact, 0.0,
  153. 0.0, 0.0, zfact,
  154. x0 - x0 * xfact, y0 - y0 * yfact, z0 - z0 * zfact)
  155. # fmt: on
  156. return affine_transform(geom, matrix)
  157. def skew(geom, xs=0.0, ys=0.0, origin="center", use_radians=False):
  158. r"""Return a skewed geometry, sheared by angles along x and y dimensions.
  159. The shear angle can be specified in either degrees (default) or radians
  160. by setting ``use_radians=True``.
  161. The point of origin can be a keyword 'center' for the bounding box
  162. center (default), 'centroid' for the geometry's centroid, a Point object
  163. or a coordinate tuple (x0, y0).
  164. The general 2D affine transformation matrix for skewing is:
  165. / 1 tan(xs) xoff \
  166. | tan(ys) 1 yoff |
  167. \ 0 0 1 /
  168. where the offsets are calculated from the origin Point(x0, y0):
  169. xoff = -y0 * tan(xs)
  170. yoff = -x0 * tan(ys)
  171. """
  172. if geom.is_empty:
  173. return geom
  174. if not use_radians: # convert from degrees
  175. xs = xs * pi / 180.0
  176. ys = ys * pi / 180.0
  177. tanx = tan(xs)
  178. tany = tan(ys)
  179. if abs(tanx) < 2.5e-16:
  180. tanx = 0.0
  181. if abs(tany) < 2.5e-16:
  182. tany = 0.0
  183. x0, y0 = interpret_origin(geom, origin, 2)
  184. # fmt: off
  185. matrix = (1.0, tanx, 0.0,
  186. tany, 1.0, 0.0,
  187. 0.0, 0.0, 1.0,
  188. -y0 * tanx, -x0 * tany, 0.0)
  189. # fmt: on
  190. return affine_transform(geom, matrix)
  191. def translate(geom, xoff=0.0, yoff=0.0, zoff=0.0):
  192. r"""Return a translated geometry shifted by offsets along each dimension.
  193. The general 3D affine transformation matrix for translation is:
  194. / 1 0 0 xoff \
  195. | 0 1 0 yoff |
  196. | 0 0 1 zoff |
  197. \ 0 0 0 1 /
  198. """
  199. if geom.is_empty:
  200. return geom
  201. # fmt: off
  202. matrix = (1.0, 0.0, 0.0,
  203. 0.0, 1.0, 0.0,
  204. 0.0, 0.0, 1.0,
  205. xoff, yoff, zoff)
  206. # fmt: on
  207. return affine_transform(geom, matrix)