Skip to content

Commit de9ae92

Browse files
committed
util to add scale bar to your plot. code from stackoverflow.
1 parent 7482f94 commit de9ae92

File tree

1 file changed

+184
-0
lines changed

1 file changed

+184
-0
lines changed

climada/util/scalebar_plot.py

Lines changed: 184 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,184 @@
1+
"""
2+
Function to plot scale bar. Source code by mephistolotl from:
3+
https://stackoverflow.com/questions/32333870/how-can-i-show-a-km-ruler-on-a-cartopy-matplotlib-plot/50674451#50674451
4+
"""
5+
6+
import numpy as np
7+
import cartopy.crs as ccrs
8+
import cartopy.geodesic as cgeo
9+
10+
def _axes_to_lonlat(ax, coords):
11+
"""(lon, lat) from axes coordinates."""
12+
display = ax.transAxes.transform(coords)
13+
data = ax.transData.inverted().transform(display)
14+
lonlat = ccrs.PlateCarree().transform_point(*data, ax.projection)
15+
16+
return lonlat
17+
18+
def _upper_bound(start, direction, distance, dist_func):
19+
"""A point farther than distance from start, in the given direction.
20+
21+
It doesn't matter which coordinate system start is given in, as long
22+
as dist_func takes points in that coordinate system.
23+
24+
Args:
25+
start: Starting point for the line.
26+
direction Nonzero (2, 1)-shaped array, a direction vector.
27+
distance: Positive distance to go past.
28+
dist_func: A two-argument function which returns distance.
29+
30+
Returns:
31+
Coordinates of a point (a (2, 1)-shaped NumPy array).
32+
"""
33+
if distance <= 0:
34+
raise ValueError(f"Minimum distance is not positive: {distance}")
35+
36+
if np.linalg.norm(direction) == 0:
37+
raise ValueError("Direction vector must not be zero.")
38+
39+
# Exponential search until the distance between start and end is
40+
# greater than the given limit.
41+
length = 0.1
42+
end = start + length * direction
43+
44+
while dist_func(start, end) < distance:
45+
length *= 2
46+
end = start + length * direction
47+
48+
return end
49+
50+
51+
def _distance_along_line(start, end, distance, dist_func, tol):
52+
"""Point at a distance from start on the segment from start to end.
53+
54+
It doesn't matter which coordinate system start is given in, as long
55+
as dist_func takes points in that coordinate system.
56+
57+
Args:
58+
start: Starting point for the line.
59+
end: Outer bound on point's location.
60+
distance: Positive distance to travel.
61+
dist_func: Two-argument function which returns distance.
62+
tol: Relative error in distance to allow.
63+
64+
Returns:
65+
Coordinates of a point (a (2, 1)-shaped NumPy array).
66+
"""
67+
initial_distance = dist_func(start, end)
68+
if initial_distance < distance:
69+
raise ValueError(f"End is closer to start ({initial_distance}) than "
70+
f"given distance ({distance}).")
71+
72+
if tol <= 0:
73+
raise ValueError(f"Tolerance is not positive: {tol}")
74+
75+
# Binary search for a point at the given distance.
76+
left = start
77+
right = end
78+
79+
while not np.isclose(dist_func(start, right), distance, rtol=tol):
80+
midpoint = (left + right) / 2
81+
82+
# If midpoint is too close, search in second half.
83+
if dist_func(start, midpoint) < distance:
84+
left = midpoint
85+
# Otherwise the midpoint is too far, so search in first half.
86+
else:
87+
right = midpoint
88+
89+
return right
90+
91+
92+
def _point_along_line(ax, start, distance, angle=0, tol=0.01):
93+
"""Point at a given distance from start at a given angle.
94+
95+
Args:
96+
ax: CartoPy axes.
97+
start: Starting point for the line in axes coordinates.
98+
distance: Positive physical distance to travel.
99+
angle: Anti-clockwise angle for the bar, in radians. Default: 0
100+
tol: Relative error in distance to allow. Default: 0.01
101+
102+
Returns:
103+
Coordinates of a point (a (2, 1)-shaped NumPy array).
104+
"""
105+
# Direction vector of the line in axes coordinates.
106+
direction = np.array([np.cos(angle), np.sin(angle)])
107+
108+
geodesic = cgeo.Geodesic()
109+
110+
# Physical distance between points.
111+
def dist_func(a_axes, b_axes):
112+
a_phys = _axes_to_lonlat(ax, a_axes)
113+
b_phys = _axes_to_lonlat(ax, b_axes)
114+
115+
# Geodesic().inverse returns a NumPy MemoryView like [[distance,
116+
# start azimuth, end azimuth]].
117+
return geodesic.inverse(a_phys, b_phys).base[0, 0]
118+
119+
end = _upper_bound(start, direction, distance, dist_func)
120+
121+
return _distance_along_line(start, end, distance, dist_func, tol)
122+
123+
124+
def scale_bar(ax, location, length, metres_per_unit=1000, unit_name='km',
125+
tol=0.01, angle=0, color='black', linewidth=3, text_offset=0.005,
126+
ha='center', va='bottom', plot_kwargs=None, text_kwargs=None,
127+
**kwargs):
128+
"""Add a scale bar to CartoPy axes.
129+
130+
For angles between 0 and 90 the text and line may be plotted at
131+
slightly different angles for unknown reasons. To work around this,
132+
override the 'rotation' keyword argument with text_kwargs.
133+
134+
Args:
135+
ax: CartoPy axes.
136+
location: Position of left-side of bar in axes coordinates.
137+
length: Geodesic length of the scale bar.
138+
metres_per_unit: Number of metres in the given unit. Default: 1000
139+
unit_name: Name of the given unit. Default: 'km'
140+
tol: Allowed relative error in length of bar. Default: 0.01
141+
angle: Anti-clockwise rotation of the bar.
142+
color: Color of the bar and text. Default: 'black'
143+
linewidth: Same argument as for plot.
144+
text_offset: Perpendicular offset for text in axes coordinates.
145+
Default: 0.005
146+
ha: Horizontal alignment. Default: 'center'
147+
va: Vertical alignment. Default: 'bottom'
148+
**plot_kwargs: Keyword arguments for plot, overridden by **kwargs.
149+
**text_kwargs: Keyword arguments for text, overridden by **kwargs.
150+
**kwargs: Keyword arguments for both plot and text.
151+
"""
152+
# Setup kwargs, update plot_kwargs and text_kwargs.
153+
if plot_kwargs is None:
154+
plot_kwargs = {}
155+
if text_kwargs is None:
156+
text_kwargs = {}
157+
158+
plot_kwargs = {'linewidth': linewidth, 'color': color, **plot_kwargs,
159+
**kwargs}
160+
text_kwargs = {'ha': ha, 'va': va, 'rotation': angle, 'color': color,
161+
**text_kwargs, **kwargs}
162+
163+
# Convert all units and types.
164+
location = np.asarray(location) # For vector addition.
165+
length_metres = length * metres_per_unit
166+
angle_rad = angle * np.pi / 180
167+
168+
# End-point of bar.
169+
end = _point_along_line(ax, location, length_metres, angle=angle_rad,
170+
tol=tol)
171+
172+
# Coordinates are currently in axes coordinates, so use transAxes to
173+
# put into data coordinates. *zip(a, b) produces a list of x-coords,
174+
# then a list of y-coords.
175+
ax.plot(*zip(location, end), transform=ax.transAxes, **plot_kwargs)
176+
177+
# Push text away from bar in the perpendicular direction.
178+
midpoint = (location + end) / 2
179+
offset = text_offset * np.array([-np.sin(angle_rad), np.cos(angle_rad)])
180+
text_location = midpoint + offset
181+
182+
# 'rotation' keyword argument is in text_kwargs.
183+
ax.text(*text_location, f"{length} {unit_name}", rotation_mode='anchor',
184+
transform=ax.transAxes, **text_kwargs)

0 commit comments

Comments
 (0)