python - 使用 Matplotlib 底图在距指定点的大圆距离内绘制阴影区域
问题描述
我想使用 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
我在哪里评论了我想做的事情,但不知道该怎么做。
我尝试了几种方法来做到这一点,但让我感到困惑的是,地图的所有输入都是以度为单位测量的坐标,而我想指定长度点,并将其转换为纬度/经度点进行绘图。我认为这与以度数为单位的纬度/经度数据与地图投影坐标有关。
任何朝着正确方向的轻推将不胜感激谢谢
解决方案
最后我不得不手动实现这个。
简而言之,我使用此处给出的方程来计算给定初始起点的坐标和径向来计算 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 海里的圆圈):
推荐阅读
- python - Inline_Formset 用于特定查询集而不是所有模型?
- reactjs - inferring arguments types based on others
- apache-storm - Apache Storm 2:无法从种子主机中找到领导灵气
- php - 为什么 Bootstrap Modal 在 Laravel 的 PHP 刀片文件中使用方括号不起作用?
- r - 通过 kable 函数将标题添加到行标签
- azure - 写入接收器的 Azure 数据工厂数据流还会创建一个空的 blob 文件
- uwp - UWP 可以在 Xbox 上向 CEC 发送命令吗?
- excel - VBA - 从文件对象获取文件自定义属性
- class - Jmeter groovy类没有卸载
- sql - 存储过程重构