Cartopy 投影上的叠加线,线上有 n 个点
Overlay line on Cartopy projection with n dots on line
我正在尝试在从指定点 A 到指定点 B 的 cartopy 投影上叠加一条线,然后让该线以设定的间隔沿路径具有 n=10 个点。我目前没有点所在的确切位置,这就是为什么我希望它们只是在设定的间隔长度上。我最接近的是通过将 x1 和 y1 设置为 nplinspace(start lat, endlat, npoints) 并使用 matplotlib 来覆盖它。然而,这画了一条直线,我希望它是弯曲的(使用 transform=ccrs.Geodetic())。如果我不使用 np.linspace,我得到了我想要的线中的曲线,但线上只有两个点而不是 10 个点。有没有办法指定这种类型的线?
这是我目前的代码(只显示了两点):
plt.figure()
ax = plt.axes(projection=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())
您正在处理沿测地线(或大致为大圆弧)的点。包 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
剧情:
我正在尝试在从指定点 A 到指定点 B 的 cartopy 投影上叠加一条线,然后让该线以设定的间隔沿路径具有 n=10 个点。我目前没有点所在的确切位置,这就是为什么我希望它们只是在设定的间隔长度上。我最接近的是通过将 x1 和 y1 设置为 nplinspace(start lat, endlat, npoints) 并使用 matplotlib 来覆盖它。然而,这画了一条直线,我希望它是弯曲的(使用 transform=ccrs.Geodetic())。如果我不使用 np.linspace,我得到了我想要的线中的曲线,但线上只有两个点而不是 10 个点。有没有办法指定这种类型的线?
这是我目前的代码(只显示了两点):
plt.figure() ax = plt.axes(projection=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())
您正在处理沿测地线(或大致为大圆弧)的点。包 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
剧情: