Cartopy 投影上的叠加线，线上有 n 个点答案

【问题标题】：Overlay line on Cartopy projection with n dots on lineCartopy 投影上的叠加线，线上有 n 个点
【发布时间】：2020-10-15 01:39:06
【问题描述】：

我试图在从指定点 A 到指定点 B 的 cartopy 投影上覆盖一条线，然后让这条线以设定的间隔沿路径具有 n=10 个点。我目前没有这些点的确切位置，这就是为什么我希望它们只是在设定的间隔长度上。我最接近的是通过将 x1 和 y1 设置为 nplinspace(start lat, endlat, npoints) 并使用 matplotlib 覆盖它。但是，这画了一条直线，我希望它是弯曲的（使用 transform=ccrs.Geodetic()）。如果我不使用 np.linspace，我会在我想要的线上得到曲线，但线上只有两个点而不是 10 个。有没有办法指定这种类型的线？

这是我目前的代码（仅显示两点）：

plt.figure()
ax = plt.axes（投影=ccrs.PlateCarree（））
ax.set_extent([-125,-60,15,65], ccrs.PlateCarree())
ax.add_feature(cfeature.LAND, color='lightgrey')

plt.plot([-120, -64], [20, 60],'o-', color='blue', transform=ccrs.Geodetic())

【问题讨论】：

标签： python matplotlib cartopy

【解决方案1】：

您正在沿着测地线（或大致大圆弧）处理点。 geographiclib 包对于沿任何测地线确定点非常方便。这是一个可运行的代码（及其输出），您可以尝试根据需要进行调整。

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from geographiclib.geodesic import Geodesic
import numpy as np

plt.figure()
proj = ccrs.PlateCarree()
proj._threshold /= 20.  #allows fine grain plot

fig = plt.figure(figsize=(8,6))
ax = fig.add_subplot(1, 1, 1, projection=proj)

ax.set_extent([-125,-60,15,65], ccrs.PlateCarree())
ax.add_feature(cfeature.LAND, color='lightgrey')

plt.plot([-120, -64], [20, 60], 'o-', color='blue',  transform=ccrs.Geodetic())

# start location
lon_fr = -120
lat_fr = 20
# end location
lon_to = -64
lat_to = 60

# This gets geodesic between the two points
# WGS84 ellipsoid is used
gl = Geodesic.WGS84.InverseLine(lat_fr, lon_fr, lat_to, lon_to)

num_points = 10  #for points on geodesic
print("distance latitude longitude azimuth")

# Compute points on the geodesic, and plot them as red dots
# gl.s13 is the total length of the geodesic
# the points are equally spaced by true distances, but not on the map
# due to the projection distortion
for ea in np.linspace(0, gl.s13, num_points):
    g = gl.Position(ea, Geodesic.STANDARD | Geodesic.LONG_UNROLL)
    print("{:.0f} {:.5f} {:.5f} {:.5f}".format(g['s12'], g['lat2'], g['lon2'], g['azi2']))
    lon2 = g['lon2']
    lat2 = g['lat2']
    ax.plot(lon2, lat2, "ro", transform=ccrs.PlateCarree())

ax.gridlines(draw_labels=True)
plt.show()

打印输出：

distance latitude longitude azimuth
0 20.00000 -120.00000 30.08493
692435 25.37542 -116.55578 31.41521
1384869 30.65898 -112.79470 33.18430
2077304 35.81710 -108.60549 35.48354
2769738 40.80499 -103.84610 38.43788
3462173 45.56121 -98.33485 42.21422
4154607 50.00000 -91.84447 47.02679
4847042 54.00165 -84.10986 53.12905
5539476 57.40348 -74.87293 60.76851
6231911 60.00000 -64.00000 70.06907

剧情：

【讨论】：