2010-01-29 15 views
16

我有一個包含一些物理量(如溫度)的極性(r,theta)網格(意思是每個單元是一個環形部分),並且我想重新網格化(或重新計劃或重新採樣)這些值到笛卡爾網格上。有沒有可以做到這一點的Python包?極地重新投影到笛卡爾網格

我對將細胞中心的座標從極座標轉換爲笛卡爾不感興趣 - 這很容易。相反,我正在尋找一個能夠正確重新格式化數據的軟件包。

感謝您的任何建議!

+0

這不是一個簡單的問題,這將是有趣和寫作的巨大熊。我想這需要2-3天的時間才能提出可怕的低效率問題。 – Omnifarious 2010-01-29 19:53:47

回答

7

謝謝您的回答 - 思維更加一下這個後,我想出了下面的代碼:

import numpy as np 

import matplotlib 
matplotlib.use('Agg') 
import matplotlib.pyplot as mpl 

from scipy.interpolate import interp1d 
from scipy.ndimage import map_coordinates 


def polar2cartesian(r, t, grid, x, y, order=3): 

    X, Y = np.meshgrid(x, y) 

    new_r = np.sqrt(X*X+Y*Y) 
    new_t = np.arctan2(X, Y) 

    ir = interp1d(r, np.arange(len(r)), bounds_error=False) 
    it = interp1d(t, np.arange(len(t))) 

    new_ir = ir(new_r.ravel()) 
    new_it = it(new_t.ravel()) 

    new_ir[new_r.ravel() > r.max()] = len(r)-1 
    new_ir[new_r.ravel() < r.min()] = 0 

    return map_coordinates(grid, np.array([new_ir, new_it]), 
          order=order).reshape(new_r.shape) 

# Define original polar grid 

nr = 10 
nt = 10 

r = np.linspace(1, 100, nr) 
t = np.linspace(0., np.pi, nt) 
z = np.random.random((nr, nt)) 

# Define new cartesian grid 

nx = 100 
ny = 200 

x = np.linspace(0., 100., nx) 
y = np.linspace(-100., 100., ny) 

# Interpolate polar grid to cartesian grid (nearest neighbor) 

fig = mpl.figure() 
ax = fig.add_subplot(111) 
ax.imshow(polar2cartesian(r, t, z, x, y, order=0), interpolation='nearest') 
fig.savefig('test1.png') 

# Interpolate polar grid to cartesian grid (cubic spline) 

fig = mpl.figure() 
ax = fig.add_subplot(111) 
ax.imshow(polar2cartesian(r, t, z, x, y, order=3), interpolation='nearest') 
fig.savefig('test2.png') 

這是沒有嚴格的重新網格化,但能正常工作我需要什麼。只要發佈代碼,以防其他人有用。隨意提出改進建議!

+0

只是一個小小的更正。我猜,它應該是arctan2(Y,X)在你的代碼中。 – 2014-09-03 09:54:38

3

你可以用scipy.ndimage.geometric_transform更緊湊地做到這一點。下面是一些示例代碼:

import numpy as N 
import scipy as S 
import scipy.ndimage 

temperature = <whatever> 
# This is the data in your polar grid. 
# The 0th and 1st axes correspond to r and θ, respectively. 
# For the sake of simplicity, θ goes from 0 to 2π, 
# and r's units are just its indices. 

def polar2cartesian(outcoords, inputshape, origin): 
    """Coordinate transform for converting a polar array to Cartesian coordinates. 
    inputshape is a tuple containing the shape of the polar array. origin is a 
    tuple containing the x and y indices of where the origin should be in the 
    output array.""" 

    xindex, yindex = outcoords 
    x0, y0 = origin 
    x = xindex - x0 
    y = yindex - y0 

    r = N.sqrt(x**2 + y**2) 
    theta = N.arctan2(y, x) 
    theta_index = N.round((theta + N.pi) * inputshape[1]/(2 * N.pi)) 

    return (r,theta_index) 

temperature_cartesian = S.ndimage.geometric_transform(temperature, polar2cartesian, 
    order=0, 
    output_shape = (temperature.shape[0] * 2, temperature.shape[0] * 2), 
    extra_keywords = {'inputshape':temperature.shape, 
     'center':(temperature.shape[0], temperature.shape[0])}) 

根據需要更好的插值可以更改order=0。輸出數組temperature_cartesian在這裏是2r乘2r,但是你可以指定你喜歡的任何大小和起點。

2

前段時間,當我嘗試做類似的事情時,我來到這篇文章,這是重新將極座標數據重新映射到笛卡爾網格,反之亦然。這裏提出的解決方案工作正常。但是,執行座標變換需要一些時間。我只是想分享另一種方法,可以將處理時間縮短50倍甚至更多。

該算法使用scipy.ndimage.interpolation.map_coordinates函數。

讓我們看一個小例子:

import numpy as np 

# Auxiliary function to map polar data to a cartesian plane 
def polar_to_cart(polar_data, theta_step, range_step, x, y, order=3): 

    from scipy.ndimage.interpolation import map_coordinates as mp 

    # "x" and "y" are numpy arrays with the desired cartesian coordinates 
    # we make a meshgrid with them 
    X, Y = np.meshgrid(x, y) 

    # Now that we have the X and Y coordinates of each point in the output plane 
    # we can calculate their corresponding theta and range 
    Tc = np.degrees(np.arctan2(Y, X)).ravel() 
    Rc = (np.sqrt(X**2 + Y**2)).ravel() 

    # Negative angles are corrected 
    Tc[Tc < 0] = 360 + Tc[Tc < 0] 

    # Using the known theta and range steps, the coordinates are mapped to 
    # those of the data grid 
    Tc = Tc/theta_step 
    Rc = Rc/range_step 

    # An array of polar coordinates is created stacking the previous arrays 
    coords = np.vstack((Ac, Sc)) 

    # To avoid holes in the 360º - 0º boundary, the last column of the data 
    # copied in the begining 
    polar_data = np.vstack((polar_data, polar_data[-1,:])) 

    # The data is mapped to the new coordinates 
    # Values outside range are substituted with nans 
    cart_data = mp(polar_data, coords, order=order, mode='constant', cval=np.nan) 

    # The data is reshaped and returned 
    return(cart_data.reshape(len(y), len(x)).T) 

polar_data = ... # Here a 2D array of data is assumed, with shape thetas x ranges 

# We create the x and y axes of the output cartesian data 
x = y = np.arange(-100000, 100000, 1000) 

# We call the mapping function assuming 1 degree of theta step and 500 meters of 
# range step. The default order of 3 is used. 
cart_data = polar_to_cart(polar_data, 1, 500, x, y) 

我希望這可以幫助別人在同樣的情況我。

0

Are there any Python packages that can do this?

是的!現在至少有一個Python包可以將矩陣從笛卡兒座標重新映射到極座標:abel.tools.polar.reproject_image_into_polar(),它是PyAbel package的一部分。

(伊尼戈Hernáez受文者是正確的,scipy.ndimage.interpolation.map_coordinates是我們迄今發現從笛卡爾重新投影到極座標的最快方式。)

PyAbel可以從PyPi通過輸入命令行下面的安裝:

pip install pyabel 

然後,在Python中,你可以使用下面的代碼重新投影的圖像爲極座標:

import abel 
abel.tools.polar.reproject_image_into_polar(MyImage) 

[根據應用程序的不同,您可能會考慮通過jacobian=True參數,該參數重新縮放矩陣的強度,以考慮從笛卡兒轉換時發生的網格拉伸(更改「bin大小」)極地coodinates。]

下面是一個完整的例子:

import numpy as np 
import matplotlib.pyplot as plt 
import abel 

CartImage = abel.tools.analytical.sample_image(501)[201:-200, 201:-200] 

PolarImage, r_grid, theta_grid = abel.tools.polar.reproject_image_into_polar(CartImage) 

fig, axs = plt.subplots(1,2, figsize=(7,3.5)) 
axs[0].imshow(CartImage , aspect='auto', origin='lower') 
axs[1].imshow(PolarImage, aspect='auto', origin='lower', 
       extent=(np.min(theta_grid), np.max(theta_grid), np.min(r_grid), np.max(r_grid))) 

axs[0].set_title('Cartesian') 
axs[0].set_xlabel('x') 
axs[0].set_ylabel('y') 

axs[1].set_title('Polar') 
axs[1].set_xlabel('Theta') 
axs[1].set_ylabel('r') 

plt.tight_layout() 
plt.show() 

enter image description here

注:還有另外一個很好的討論(有關重新映射彩色圖像極座標)上SO:image information along a polar coordinate system

+0

這是一個很好的例子。然而它給了我'TypeError:'numpy.float64'對象不能被解釋爲一個整數'在python3.4上。如果你是代碼的維護者,你應該檢查一下。 – TomCho 2017-06-05 20:31:46