使用 Matplotlib 底图在距指定点的大圆距离内绘制阴影区域

Question

我想使用 Python/Matplotlib/Basemap 绘制地图并对位于指定点的给定距离内的圆进行着色，类似于此（由 Great Circle Mapper 生成的地图 - 版权所有 © Karl L. Swartz。 ):

我可以让地图生成如下：

from mpl_toolkits.basemap import Basemap
import numpy as np
import matplotlib.pyplot as plt

# create new figure, axes instances.
fig,ax = plt.subplots()

# setup Mercator map projection.
m = Basemap(
            llcrnrlat=47.0,
            llcrnrlon=-126.62,
            urcrnrlat=50.60,
            urcrnrlon=-119.78,
            rsphere=(6378137.00,6356752.3142),
            resolution='i',
            projection='merc',
            lat_0=49.290,
            lon_0=-123.117,
            )

# Latitudes and longitudes of locations of interest
coords = dict()
coords['SEA'] = [47.450, -122.309]

# Plot markers and labels on map
for key in coords:
    lon, lat = coords[key]
    x,y = m(lat, lon)
    m.plot(x, y, 'bo', markersize=5)
    plt.text(x+10000, y+5000, key, color='k')

# Draw in coastlines
m.drawcoastlines()
m.fillcontinents()

m.fillcontinents(color='grey',lake_color='aqua')
m.drawmapboundary(fill_color='aqua')

plt.show()

生成地图：

现在我想在指定点周围创建一个大圆圈，例如顶部地图。

我的尝试是一个函数，它接受地图对象、中心坐标对和距离，并创建两条曲线，然后在它们之间着色，类似于：

def shaded_great_circle(map_, lat_0, lon_0, dist=100, alpha=0.2):  # dist specified in nautical miles
    dist = dist * 1852  # Convert distance to nautical miles
    lat = np.linspace(lat_0-dist/2, lat_0+dist/2,50)
    lon = # Somehow find these points
    # Create curve for longitudes above lon_0
    # Create curve for longitudes below lon_0
    # Shade region between above two curves

我已经评论了我想做的事情，但不知道该怎么做。

我尝试了几种方法来做到这一点，但让我感到困惑的是，地图的所有输入都是以度为单位的坐标，而我想指定长度点，并将其转换为纬度/经度点阴谋。我认为这与以度数为单位的纬度/经度数据与地图投影坐标有关。

任何朝着正确方向的轻推将不胜感激 谢谢

Answer 1

最后我不得不手动实现。

简而言之，我使用给定here 的方程来计算给定的坐标 初始起点和径向计算 360 度左右的点，然后通过这些点绘制一条线。我真的不需要阴影部分，所以我还没有实现。

我认为这是一个有用的功能，所以我是这样实现的：

from mpl_toolkits.basemap import Basemap
import numpy as np
import matplotlib.pyplot as plt


def calc_new_coord(lat1, lon1, rad, dist):
    """
    Calculate coordinate pair given starting point, radial and distance
    Method from: http://www.geomidpoint.com/destination/calculation.html
    """

    flat = 298.257223563
    a = 2 * 6378137.00
    b = 2 * 6356752.3142

    # Calculate the destination point using Vincenty's formula
    f = 1 / flat
    sb = np.sin(rad)
    cb = np.cos(rad)
    tu1 = (1 - f) * np.tan(lat1)
    cu1 = 1 / np.sqrt((1 + tu1*tu1))
    su1 = tu1 * cu1
    s2 = np.arctan2(tu1, cb)
    sa = cu1 * sb
    csa = 1 - sa * sa
    us = csa * (a * a - b * b) / (b * b)
    A = 1 + us / 16384 * (4096 + us * (-768 + us * (320 - 175 * us)))
    B = us / 1024 * (256 + us * (-128 + us * (74 - 47 * us)))
    s1 = dist / (b * A)
    s1p = 2 * np.pi

    while (abs(s1 - s1p) > 1e-12):
        cs1m = np.cos(2 * s2 + s1)
        ss1 = np.sin(s1)
        cs1 = np.cos(s1)
        ds1 = B * ss1 * (cs1m + B / 4 * (cs1 * (- 1 + 2 * cs1m * cs1m) - B / 6 * \
            cs1m * (- 3 + 4 * ss1 * ss1) * (-3 + 4 * cs1m * cs1m)))
        s1p = s1
        s1 = dist / (b * A) + ds1

    t = su1 * ss1 - cu1 * cs1 * cb
    lat2 = np.arctan2(su1 * cs1 + cu1 * ss1 * cb, (1 - f) * np.sqrt(sa * sa + t * t))
    l2 = np.arctan2(ss1 * sb, cu1 * cs1 - su1 * ss1 * cb)
    c = f / 16 * csa * (4 + f * (4 - 3 * csa))
    l = l2 - (1 - c) * f * sa * (s1 + c * ss1 * (cs1m + c * cs1 * (-1 + 2 * cs1m * cs1m)))
    d = np.arctan2(sa, -t)
    finaltc = d + 2 * np.pi
    backtc = d + np.pi
    lon2 = lon1 + l

    return (np.rad2deg(lat2), np.rad2deg(lon2))


def shaded_great_circle(m, lat_0, lon_0, dist=100, alpha=0.2, col='k'):  # dist specified in nautical miles
    dist = dist * 1852  # Convert distance to nautical miles
    theta_arr = np.linspace(0, np.deg2rad(360), 100)
    lat_0 = np.deg2rad(lat_0)
    lon_0 = np.deg2rad(lon_0)

    coords_new = []

    for theta in theta_arr:
        coords_new.append(calc_new_coord(lat_0, lon_0, theta, dist))

    lat = [item[0] for item in coords_new]
    lon = [item[1] for item in coords_new]

    x, y = m(lon, lat)
    m.plot(x, y, col)

# setup Mercator map projection.
m = Basemap(
            llcrnrlat=45.0,
            llcrnrlon=-126.62,
            urcrnrlat=50.60,
            urcrnrlon=-119.78,
            rsphere=(6378137.00,6356752.3142),
            resolution='i',
            projection='merc',
            lat_0=49.290,
            lon_0=-123.117,
            )

# Latitudes and longitudes of locations of interest
coords = dict()
coords['SEA'] = [47.450, -122.309]

# Plot markers and labels on map
for key in coords:
    lon, lat = coords[key]
    x,y = m(lat, lon)
    m.plot(x, y, 'bo', markersize=5)
    plt.text(x+10000, y+5000, key, color='k')

# Draw in coastlines
m.drawcoastlines()
m.fillcontinents()

m.fillcontinents(color='grey',lake_color='aqua')
m.drawmapboundary(fill_color='aqua')

# Draw great circle
shaded_great_circle(m, 47.450, -122.309, 100, col='k')  # Distance specified in nautical miles, i.e. 100 nmi in this case

plt.show()

运行这个应该会给你（西雅图周围 100 海里的圆圈）：