如何在cartopy/matplotlib图上显示km标尺？
gau*_*teh 14 python matplotlib cartopy
如何为地图的放大部分显示公里标尺,无论是插入图像还是作为绘图侧面的标尺？
例如,侧面(左侧)的50公里栏或mi(右侧)的插图:
(问题:cartopy#490)
这是我为自己编写的Cartopy比例尺函数,它使用了更简单的pp-mo的答案:编辑:修改代码创建一个新的投影居中,使得比例尺平行于许多坐标系的轴,包括一些正交和更大的地图,并删除指定utm系统的需要.如果未指定比例尺长度,还会添加代码来计算比例尺长度.
import cartopy.crs as ccrs
import numpy as np
def scale_bar(ax, length=None, location=(0.5, 0.05), linewidth=3):
"""
ax is the axes to draw the scalebar on.
length is the length of the scalebar in km.
location is center of the scalebar in axis coordinates.
(ie. 0.5 is the middle of the plot)
linewidth is the thickness of the scalebar.
"""
#Get the limits of the axis in lat long
llx0, llx1, lly0, lly1 = ax.get_extent(ccrs.PlateCarree())
#Make tmc horizontally centred on the middle of the map,
#vertically at scale bar location
sbllx = (llx1 + llx0) / 2
sblly = lly0 + (lly1 - lly0) * location[1]
tmc = ccrs.TransverseMercator(sbllx, sblly)
#Get the extent of the plotted area in coordinates in metres
x0, x1, y0, y1 = ax.get_extent(tmc)
#Turn the specified scalebar location into coordinates in metres
sbx = x0 + (x1 - x0) * location[0]
sby = y0 + (y1 - y0) * location[1]
#Calculate a scale bar length if none has been given
#(Theres probably a more pythonic way of rounding the number but this works)
if not length:
length = (x1 - x0) / 5000 #in km
ndim = int(np.floor(np.log10(length))) #number of digits in number
length = round(length, -ndim) #round to 1sf
#Returns numbers starting with the list
def scale_number(x):
if str(x)[0] in ['1', '2', '5']: return int(x)
else: return scale_number(x - 10 ** ndim)
length = scale_number(length)
#Generate the x coordinate for the ends of the scalebar
bar_xs = [sbx - length * 500, sbx + length * 500]
#Plot the scalebar
ax.plot(bar_xs, [sby, sby], transform=tmc, color='k', linewidth=linewidth)
#Plot the scalebar label
ax.text(sbx, sby, str(length) + ' km', transform=tmc,
horizontalalignment='center', verticalalignment='bottom')
它有一些限制,但相对简单,所以我希望你能看到如果你想要不同的东西如何改变它.
用法示例:
import matplotlib.pyplot as plt
ax = plt.axes(projection=ccrs.Mercator())
plt.title('Cyprus')
ax.set_extent([31, 35.5, 34, 36], ccrs.Geodetic())
ax.coastlines(resolution='10m')
scale_bar(ax, 100)
plt.show()
小智 9
通过在CartoPy 0.15中添加测地模块,我们现在可以非常轻松地在地图上计算精确的长度.想弄清楚如何在直线上的地图线上找到两个点是一个有点棘手的问题.一旦指定了地图上的方向,我就会执行指数搜索以找到足够远的点.然后我执行二分搜索以找到足够接近所需距离的点.
该scale_bar功能是很简单的,但它有很多的选择.基本签名是scale_bar(ax, location, length).ax是任何CartoPy轴,location是轴坐标中条形左侧的位置(因此每个坐标从0到1),并且length是以公里为单位的条形长度.支持其他长度,如metres_per_unit和unit_name关键字参数一样.
额外的关键字参数(如color)简单地传递给text和plot.但是,特定于text或plot(例如family或path_effects)的关键字参数必须通过text_kwargs和作为字典传递plot_kwargs.
我已经包含了我认为常见用例的示例.
请分享任何问题,评论或批评.
scalebar.py
import numpy as np
import cartopy.crs as ccrs
import cartopy.geodesic as cgeo
def _axes_to_lonlat(ax, coords):
"""(lon, lat) from axes coordinates."""
display = ax.transAxes.transform(coords)
data = ax.transData.inverted().transform(display)
lonlat = ccrs.PlateCarree().transform_point(*data, ax.projection)
return lonlat
def _upper_bound(start, direction, distance, dist_func):
"""A point farther than distance from start, in the given direction.
It doesn't matter which coordinate system start is given in, as long
as dist_func takes points in that coordinate system.
Args:
start: Starting point for the line.
direction Nonzero (2, 1)-shaped array, a direction vector.
distance: Positive distance to go past.
dist_func: A two-argument function which returns distance.
Returns:
Coordinates of a point (a (2, 1)-shaped NumPy array).
"""
if distance <= 0:
raise ValueError(f"Minimum distance is not positive: {distance}")
if np.linalg.norm(direction) == 0:
raise ValueError("Direction vector must not be zero.")
# Exponential search until the distance between start and end is
# greater than the given limit.
length = 0.1
end = start + length * direction
while dist_func(start, end) < distance:
length *= 2
end = start + length * direction
return end
def _distance_along_line(start, end, distance, dist_func, tol):
"""Point at a distance from start on the segment from start to end.
It doesn't matter which coordinate system start is given in, as long
as dist_func takes points in that coordinate system.
Args:
start: Starting point for the line.
end: Outer bound on point's location.
distance: Positive distance to travel.
dist_func: Two-argument function which returns distance.
tol: Relative error in distance to allow.
Returns:
Coordinates of a point (a (2, 1)-shaped NumPy array).
"""
initial_distance = dist_func(start, end)
if initial_distance < distance:
raise ValueError(f"End is closer to start ({initial_distance}) than "
f"given distance ({distance}).")
if tol <= 0:
raise ValueError(f"Tolerance is not positive: {tol}")
# Binary search for a point at the given distance.
left = start
right = end
while not np.isclose(dist_func(start, right), distance, rtol=tol):
midpoint = (left + right) / 2
# If midpoint is too close, search in second half.
if dist_func(start, midpoint) < distance:
left = midpoint
# Otherwise the midpoint is too far, so search in first half.
else:
right = midpoint
return right
def _point_along_line(ax, start, distance, angle=0, tol=0.01):
"""Point at a given distance from start at a given angle.
Args:
ax: CartoPy axes.
start: Starting point for the line in axes coordinates.
distance: Positive physical distance to travel.
angle: Anti-clockwise angle for the bar, in radians. Default: 0
tol: Relative error in distance to allow. Default: 0.01
Returns:
Coordinates of a point (a (2, 1)-shaped NumPy array).
"""
# Direction vector of the line in axes coordinates.
direction = np.array([np.cos(angle), np.sin(angle)])
geodesic = cgeo.Geodesic()
# Physical distance between points.
def dist_func(a_axes, b_axes):
a_phys = _axes_to_lonlat(ax, a_axes)
b_phys = _axes_to_lonlat(ax, b_axes)
# Geodesic().inverse returns a NumPy MemoryView like [[distance,
# start azimuth, end azimuth]].
return geodesic.inverse(a_phys, b_phys).base[0, 0]
end = _upper_bound(start, direction, distance, dist_func)
return _distance_along_line(start, end, distance, dist_func, tol)
def scale_bar(ax, location, length, metres_per_unit=1000, unit_name='km',
tol=0.01, angle=0, color='black', linewidth=3, text_offset=0.005,
ha='center', va='bottom', plot_kwargs=None, text_kwargs=None,
**kwargs):
"""Add a scale bar to CartoPy axes.
For angles between 0 and 90 the text and line may be plotted at
slightly different angles for unknown reasons. To work around this,
override the 'rotation' keyword argument with text_kwargs.
Args:
ax: CartoPy axes.
location: Position of left-side of bar in axes coordinates.
length: Geodesic length of the scale bar.
metres_per_unit: Number of metres in the given unit. Default: 1000
unit_name: Name of the given unit. Default: 'km'
tol: Allowed relative error in length of bar. Default: 0.01
angle: Anti-clockwise rotation of the bar.
color: Color of the bar and text. Default: 'black'
linewidth: Same argument as for plot.
text_offset: Perpendicular offset for text in axes coordinates.
Default: 0.005
ha: Horizontal alignment. Default: 'center'
va: Vertical alignment. Default: 'bottom'
**plot_kwargs: Keyword arguments for plot, overridden by **kwargs.
**text_kwargs: Keyword arguments for text, overridden by **kwargs.
**kwargs: Keyword arguments for both plot and text.
"""
# Setup kwargs, update plot_kwargs and text_kwargs.
if plot_kwargs is None:
plot_kwargs = {}
if text_kwargs is None:
text_kwargs = {}
plot_kwargs = {'linewidth': linewidth, 'color': color, **plot_kwargs,
**kwargs}
text_kwargs = {'ha': ha, 'va': va, 'rotation': angle, 'color': color,
**text_kwargs, **kwargs}
# Convert all units and types.
location = np.asarray(location) # For vector addition.
length_metres = length * metres_per_unit
angle_rad = angle * np.pi / 180
# End-point of bar.
end = _point_along_line(ax, location, length_metres, angle=angle_rad,
tol=tol)
# Coordinates are currently in axes coordinates, so use transAxes to
# put into data coordinates. *zip(a, b) produces a list of x-coords,
# then a list of y-coords.
ax.plot(*zip(location, end), transform=ax.transAxes, **plot_kwargs)
# Push text away from bar in the perpendicular direction.
midpoint = (location + end) / 2
offset = text_offset * np.array([-np.sin(angle_rad), np.cos(angle_rad)])
text_location = midpoint + offset
# 'rotation' keyword argument is in text_kwargs.
ax.text(*text_location, f"{length} {unit_name}", rotation_mode='anchor',
transform=ax.transAxes, **text_kwargs)
demo.py
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
from scalebar import scale_bar
fig = plt.figure(1, figsize=(10, 10))
ax = fig.add_subplot(111, projection=ccrs.Mercator())
ax.set_extent([-180, 180, -85, 85])
ax.coastlines(facecolor='black')
ax.add_feature(cfeature.LAND)
# Standard 6,000 km scale bar.
scale_bar(ax, (0.65, 0.4), 6_000)
# Length of the bar reflects its position on the map.
scale_bar(ax, (0.55, 0.7), 6_000, color='green')
# Bar can be placed at any angle. Any units can be used.
scale_bar(ax, (0.4, 0.4), 3_000, metres_per_unit=1609, angle=-90,
unit_name='mi', color='red')
# Text and line can be styled separately. Keywords are simply passed to
# text or plot.
text_kwargs = dict(family='serif', size='xx-large', color='red')
plot_kwargs = dict(linestyle='dashed', color='blue')
scale_bar(ax, (0.05, 0.3), 6_000, text_kwargs=text_kwargs,
plot_kwargs=plot_kwargs)
# Angles between 0 and 90 may result in the text and line plotted at
# slightly different angles for an unknown reason.
scale_bar(ax, (0.45, 0.15), 5_000, color='purple', angle=45, text_offset=0)
# To get around this override the text's angle and fiddle manually.
scale_bar(ax, (0.55, 0.15), 5_000, color='orange', angle=45, text_offset=0,
text_kwargs={'rotation': 41})
plt.show()
这是@Siyh答案的精炼版本，其中添加了：
- 自动UTM区域选择
- 文本/栏后面的缓冲区，使其在背景下显示
- 北箭
笔记：
- 如果您的轴不使用UTM，则条形图会弯曲
- 北箭头假定北朝上
码：
import os
import cartopy.crs as ccrs
from math import floor
import matplotlib.pyplot as plt
from matplotlib import patheffects
import matplotlib
if os.name == 'nt':
matplotlib.rc('font', family='Arial')
else: # might need tweaking, must support black triangle for N arrow
matplotlib.rc('font', family='DejaVu Sans')
def utm_from_lon(lon):
"""
utm_from_lon - UTM zone for a longitude
Not right for some polar regions (Norway, Svalbard, Antartica)
:param float lon: longitude
:return: UTM zone number
:rtype: int
"""
return floor( ( lon + 180 ) / 6) + 1
def scale_bar(ax, proj, length, location=(0.5, 0.05), linewidth=3,
units='km', m_per_unit=1000):
"""
http://stackoverflow.com/a/35705477/1072212
ax is the axes to draw the scalebar on.
proj is the projection the axes are in
location is center of the scalebar in axis coordinates ie. 0.5 is the middle of the plot
length is the length of the scalebar in km.
linewidth is the thickness of the scalebar.
units is the name of the unit
m_per_unit is the number of meters in a unit
"""
# find lat/lon center to find best UTM zone
x0, x1, y0, y1 = ax.get_extent(proj.as_geodetic())
# Projection in metres
utm = ccrs.UTM(utm_from_lon((x0+x1)/2))
# Get the extent of the plotted area in coordinates in metres
x0, x1, y0, y1 = ax.get_extent(utm)
# Turn the specified scalebar location into coordinates in metres
sbcx, sbcy = x0 + (x1 - x0) * location[0], y0 + (y1 - y0) * location[1]
# Generate the x coordinate for the ends of the scalebar
bar_xs = [sbcx - length * m_per_unit/2, sbcx + length * m_per_unit/2]
# buffer for scalebar
buffer = [patheffects.withStroke(linewidth=5, foreground="w")]
# Plot the scalebar with buffer
ax.plot(bar_xs, [sbcy, sbcy], transform=utm, color='k',
linewidth=linewidth, path_effects=buffer)
# buffer for text
buffer = [patheffects.withStroke(linewidth=3, foreground="w")]
# Plot the scalebar label
t0 = ax.text(sbcx, sbcy, str(length) + ' ' + units, transform=utm,
horizontalalignment='center', verticalalignment='bottom',
path_effects=buffer, zorder=2)
left = x0+(x1-x0)*0.05
# Plot the N arrow
t1 = ax.text(left, sbcy, u'\u25B2\nN', transform=utm,
horizontalalignment='center', verticalalignment='bottom',
path_effects=buffer, zorder=2)
# Plot the scalebar without buffer, in case covered by text buffer
ax.plot(bar_xs, [sbcy, sbcy], transform=utm, color='k',
linewidth=linewidth, zorder=3)
if __name__ == '__main__':
ax = plt.axes(projection=ccrs.Mercator())
plt.title('Cyprus')
ax.set_extent([31, 35.5, 34, 36], ccrs.Geodetic())
ax.stock_img()
ax.coastlines(resolution='10m')
scale_bar(ax, ccrs.Mercator(), 100) # 100 km scale bar
# or to use m instead of km
# scale_bar(ax, ccrs.Mercator(), 100000, m_per_unit=1, units='m')
# or to use miles instead of km
# scale_bar(ax, ccrs.Mercator(), 60, m_per_unit=1609.34, units='miles')
plt.show()
小智 2
我认为对此没有简单的盆栽解决方案：您必须使用图形元素自己将其绘制出来。
几年前，我编写了一些自适应代码，将比例尺添加到任意比例的操作系统网格地图中。
我认为这并不是您想要的，但它显示了必要的技术：
def add_osgb_scalebar(ax, at_x=(0.1, 0.4), at_y=(0.05, 0.075), max_stripes=5):
"""
Add a scalebar to a GeoAxes of type cartopy.crs.OSGB (only).
Args:
* at_x : (float, float)
target axes X coordinates (0..1) of box (= left, right)
* at_y : (float, float)
axes Y coordinates (0..1) of box (= lower, upper)
* max_stripes
typical/maximum number of black+white regions
"""
# ensure axis is an OSGB map (meaning coords are just metres)
assert isinstance(ax.projection, ccrs.OSGB)
# fetch axes coordinate mins+maxes
x0, x1 = ax.get_xlim()
y0, y1 = ax.get_ylim()
# set target rectangle in-visible-area (aka 'Axes') coordinates
ax0, ax1 = at_x
ay0, ay1 = at_y
# choose exact X points as sensible grid ticks with Axis 'ticker' helper
x_targets = [x0 + ax * (x1 - x0) for ax in (ax0, ax1)]
ll = mpl.ticker.MaxNLocator(nbins=max_stripes, steps=[1,2,4,5,10])
x_vals = ll.tick_values(*x_targets)
# grab min+max for limits
xl0, xl1 = x_vals[0], x_vals[-1]
# calculate Axes Y coordinates of box top+bottom
yl0, yl1 = [y0 + ay * (y1 - y0) for ay in [ay0, ay1]]
# calculate Axes Y distance of ticks + label margins
y_margin = (yl1-yl0)*0.25
# fill black/white 'stripes' and draw their boundaries
fill_colors = ['black', 'white']
i_color = 0
for xi0, xi1 in zip(x_vals[:-1],x_vals[1:]):
# fill region
plt.fill((xi0, xi1, xi1, xi0, xi0), (yl0, yl0, yl1, yl1, yl0),
fill_colors[i_color])
# draw boundary
plt.plot((xi0, xi1, xi1, xi0, xi0), (yl0, yl0, yl1, yl1, yl0),
'black')
i_color = 1 - i_color
# add short tick lines
for x in x_vals:
plt.plot((x, x), (yl0, yl0-y_margin), 'black')
# add a scale legend 'Km'
font_props = mfonts.FontProperties(size='medium', weight='bold')
plt.text(
0.5 * (xl0 + xl1),
yl1 + y_margin,
'Km',
verticalalignment='bottom',
horizontalalignment='center',
fontproperties=font_props)
# add numeric labels
for x in x_vals:
plt.text(x,
yl0 - 2 * y_margin,
'{:g}'.format((x - xl0) * 0.001),
verticalalignment='top',
horizontalalignment='center',
fontproperties=font_props)
虽然很乱，不是吗？
您可能认为可以为此添加某种“浮动轴对象”，以提供自动自动缩放图形，但我无法找到一种方法来做到这一点（我想我仍然无法做到这一点） t）。
华泰
| 归档时间： |
|
| 查看次数： |
3949 次 |
| 最近记录： |