Нарисуйте заштрихованную область на расстоянии большого круга от указанной точки с помощью Matplotlib Basemap

Question

Я хочу использовать Python / Matplotlib / Basemap, чтобы нарисовать карту и заштриховать окружность, которая находится в пределах заданного расстояния от указанной точки, подобно этой (Карта, сгенерированная Картографом Великого круга - авторское право © 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

В конце концов мне пришлось реализовать это вручную.

Короче говоря, я использовал уравнение, заданное здесь , чтобы вычислить координаты, заданные начальной начальной точкой, и радиал, чтобы вычислить точки вокруг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 морских миль вокруг Сиэтла):