用橢圓顯示隨時間變化的矢量旋轉

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()