用椭圆显示随时间变化的矢量旋转

Question

我有一个随时间旋转的风矢量网格。我有这些矢量的u和v组件（东西和南北）。

这里是一个例子在一个网格点覆盖所有时间的矢量。

quiver(zeros((8,1)), zeros((8,1)),u[:,1,1], v[:,1,1]) 

我想总结这些向量的旋转在一个情节，通过在每个网格点，基本上跟踪向量的路径随时间绘制椭圆。

我基本上找什么做的是在这个情节在这里完成： https://mdc.coaps.fsu.edu/scatterometry/meeting/docs/2015/NewProductsAndApplications/gille_ovwst15.pdf

的椭圆稍微淡淡的，但他们在那里。

我猜我应该以某种方式使用matplotlib.patches.ellipse，但我不知道如何让椭圆的角度我的数据。

Answer 1

这个问题有两个主要组成部分。

1. 拟合椭圆。 在this site上发现的Python中，实际上有一个很好的例子来说明如何将椭圆拟合到数据点。所以我们可以使用它从数据中获得椭圆的旋转角度和2维。

2. 将所有椭圆图绘制成图。一旦获得椭圆的参数，椭圆可以使用matplotlib.patches.Ellipse

绘制下面是完整的代码：

import numpy as np 
from numpy.linalg import eig, inv 
import matplotlib.pyplot as plt 
from matplotlib.patches import Ellipse 

###################### 
### Ellipse fitting ## 
###################### 
# taken from 
# http://nicky.vanforeest.com/misc/fitEllipse/fitEllipse.html 

def fitEllipse(x,y): 
    x = x[:,np.newaxis] 
    y = y[:,np.newaxis] 
    D = np.hstack((x*x, x*y, y*y, x, y, np.ones_like(x))) 
    S = np.dot(D.T,D) 
    C = np.zeros([6,6]) 
    C[0,2] = C[2,0] = 2; C[1,1] = -1 
    try: 
     E, V = eig(np.dot(inv(S), C)) 
     n = np.argmax(np.abs(E)) 
     a = V[:,n] 
     return a 
    except: 
     return [np.nan]*5 

def ellipse_center(a): 
    b,c,d,f,g,a = a[1]/2, a[2], a[3]/2, a[4]/2, a[5], a[0] 
    num = b*b-a*c 
    x0=(c*d-b*f)/num 
    y0=(a*f-b*d)/num 
    return np.real(np.array([x0,y0])) 


def ellipse_angle_of_rotation(a): 
    b,c,d,f,g,a = a[1]/2, a[2], a[3]/2, a[4]/2, a[5], a[0] 
    return np.real(0.5*np.arctan(2*b/(a-c))) 


def ellipse_axis_length(a): 
    b,c,d,f,g,a = a[1]/2, a[2], a[3]/2, a[4]/2, a[5], a[0] 
    up = 2*(a*f*f+c*d*d+g*b*b-2*b*d*f-a*c*g) 
    down1=(b*b-a*c)*((c-a)*np.sqrt(1+4*b*b/((a-c)*(a-c)))-(c+a)) 
    down2=(b*b-a*c)*((a-c)*np.sqrt(1+4*b*b/((a-c)*(a-c)))-(c+a)) 
    res1=np.sqrt(up/down1) 
    res2=np.sqrt(up/down2) 
    return np.real(np.array([res1, res2])) 


######################## 
### Data Generation ### 
######################## 

n_el = 8 # number of ellipse points 
# define grid 
x = np.linspace(-7,7, 15) 
y = np.linspace(4,18, 15) 
# data (2 for x,y (west, north), n_el, dimensions of grid in x and y) 
data = np.zeros((2, n_el,len(x), len(y))) 
for i in range(len(y)): 
    for j in range(len(x)): 
     #generate n_el points on an ellipse 
     r = np.linspace(0,2*np.pi, n_el) 
     data[0,:,j,i] = 0.5*(0.9*np.random.random(1)+0.1) * np.cos(r+2*np.random.random(1)*np.pi) 
     data[1,:,j,i] = 0.5*(0.9*np.random.random(1)+0.1) * np.sin(r) 

# Test case: fit an ellipse and print the parameters 
a = fitEllipse(data[0,:,0,0], data[1,:,0,0]) 
ang = ellipse_angle_of_rotation(a) 
l = ellipse_axis_length(a) 
center = ellipse_center(a) 
print "\tangle: {}\n\tlength: {}\n\tcenter: {}".format(ang, l, center) 




###################### 
####### plotting ### 
###################### 
fig, (ax, ax2) = plt.subplots(ncols=2, figsize=(11,5)) 

# First, draw the test case ellipse 
# raw data 
ax.scatter(data[0,:,0,0], data[1,:,0,0], s=30, c="r", zorder=10) 

# Fitted Ellipse 
# matplotlib.patches.Ellipse 
# http://matplotlib.org/api/patches_api.html#matplotlib.patches.Ellipse 
# takes width and height as diameter instead of half width and rotation in degrees 
e = Ellipse(xy=(0,0), width=2*l[0], height=2*l[1], angle=ang*180./np.pi, facecolor="b", alpha=0.2, zorder=0) 
ec = Ellipse(xy=(0,0), width=2*l[0], height=2*l[1], angle=ang*180./np.pi, fill=False, zorder=1) 
ax.add_artist(e) 
ax.add_artist(ec) 

ax.set_aspect("equal") 
ax.set_xlim([-1,1]) 
ax.set_ylim([-1,1]) 

# Fit ellipse for every datapoint on grid and place in figure 
for i in range(len(y)): 
    for j in range(len(x)): 
     a = fitEllipse(data[0,:,j,i], data[1,:,j,i]) 
     ang = ellipse_angle_of_rotation(a) 
     l = ellipse_axis_length(a) 
     e = Ellipse(xy=(x[j],y[i]), width=2*l[0], height=2*l[1], angle=ang*180./np.pi, fill=False, zorder=1) 
     ax2.add_artist(e) 

ax2.set_ylim([y.min()-0.5, y.max()+0.5 ]) 
ax2.set_xlim([x.min()-0.5, x.max()+0.5 ]) 
ax2.set_aspect("equal") 

# load some background image. 
image = "https://upload.wikimedia.org/wikipedia/commons/thumb/c/ca/Singapore-OutlineMap-20050606.png/600px-Singapore-OutlineMap-20050606.png" 
image = np.rot90(plt.imread(image)) 
im = ax2.imshow(image,extent=[x.min()-0.5, x.max()+0.5, y.min()-0.5, y.max()+0.5, ])   

plt.show()