python - 用椭圆显示矢量随时间的旋转

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. 将所有椭圆绘制成图形。获得椭圆参数后，可以使用 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()