如何将向量旋转到Basemap地图投影上？

Question

我正在使用Python中的matplotlib库和basemap工具制作风向矢量图。

我有一组风的观测数据，这些数据是在任意的纬度和经度上收集的，也就是说，它们并不是在一个规则的网格上。

在绘图之前，我需要把这些矢量旋转到地图的投影上，否则风向的箭头会指错方向。请问有什么好的方法可以做到这一点呢？

例如：

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

# Define locations of my vectors
lat = numpy.array([50.1,46.2,51.6,52.2,54.4])
lon = numpy.array([-3.3,-1.0,-5.2,-1.2,0.2])

# Define some east-west vectors to illustrate the problem
u = numpy.array([5,5,5,5,5])
v = numpy.array([0,0,0,0,0])

# Set up map projection
m = Basemap(llcrnrlon=-15.,llcrnrlat=46.,urcrnrlon=15.,urcrnrlat=59.,
            projection='lcc',lat_1=40.,lat_2=50.,lon_0=-50.,
            resolution ='l')

# Calculate positions of vectors on map projection 
x,y = m(lon,lat)

# Draw barbs
m.barbs(x,y,u,v, length=7, color='red')

# Draw some grid lines for reference
parallels = numpy.arange(-80.,90,20.)
meridians = numpy.arange(0.,360.,20.)
m.drawparallels(parallels)
m.drawmeridians(meridians)
m.drawcoastlines(linewidth=0.5)

plt.show()

注意，在这个图中，矢量并不是指向东西方向。

我尝试过使用rotate_vector和transform_vector这两个函数，但它们只适用于有规则网格的数据。

有没有什么方法可以把这些任意的纬度、经度和风速、风向的配对数据旋转到地图投影上呢？

任何帮助都将不胜感激！

Answer 1

对于那些有网格数据的人，如果你碰到了这个问题

建议使用内置的函数rotate_vector，你可以在这里找到它：

http://matplotlib.org/basemap/api/basemap_api.html

Answer 2

你的问题在于，你把 u 和 v 用的是经纬度，而 x 和 y 用的是地图坐标。barbs 似乎希望这两者都用地图坐标，而不是混合使用。

最简单的方法就是计算出端点来获取这些分量。（我之前的描述可能不太清楚，所以这里给你一个更清晰的思路：）

x, y = m(lon, lat)
x1, y1 = m(lon+u, lat+v)
u_map, v_map = x1-x, y1-y

你还需要重新调整大小，因为幅度也需要处理。下面是一个完整的例子：

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

# Define locations of my vectors
lat = numpy.array([50.1,46.2,51.6,52.2,54.4])
lon = numpy.array([-3.3,-1.0,-5.2,-1.2,0.2])

# Define some east-west vectors to illustrate the problem
u = numpy.array([5,5,5,5,5])
v = numpy.array([0,0,0,0,0])

# Set up map projection
m = Basemap(llcrnrlon=-15.,llcrnrlat=46.,urcrnrlon=15.,urcrnrlat=59.,
            projection='lcc',lat_1=40.,lat_2=50.,lon_0=-50.,
            resolution ='l')

# Calculate positions of vectors on map projection 
x,y = m(lon,lat)

# Calculate the orientation of the vectors
x1, y1 = m(lon+u, lat+v)
u_map, v_map = x1-x, y1-y

# Rescale the magnitudes of the vectors...
mag_scale = np.hypot(u_map, v_map) / np.hypot(u, v)
u_map /= mag_scale
v_map /= mag_scale

# Draw barbs
m.barbs(x,y,u_map,v_map, length=7, color='red')

# Draw some grid lines for reference
parallels = numpy.arange(-80.,90,20.)
meridians = numpy.arange(0.,360.,20.)
m.drawparallels(parallels)
m.drawmeridians(meridians)
m.drawcoastlines(linewidth=0.5)

plt.show()