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

2024-05-19 • 问答

我可以按如下方式生成地图：

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

我已在其中评论了我想做什么，但不确定如何做。

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

任何朝着正确方向前进的人都会受到赞赏谢谢

最后，我必须手动执行此操作。

简而言之，我使用给定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,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,'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海里的圈）：

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

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

大家都在问