2010-01-29 72 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