Des polaires à la grille cartésienne

Merci pour vos réponses - après avoir réfléchi un peu plus à propos de ce que je suis venu avec le code suivant:

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

Qui n'est pas strictement re-maillage, mais fonctionne très bien pour ce dont j'ai besoin. Simplement poster le code au cas où il est utile à quelqu'un d'autre. N'hésitez pas à suggérer des améliorations!

vous pouvez le faire de façon plus compacte avec scipy.ndimage.geometric_transform. Voici un exemple de code:

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])})

vous pouvez changer order=0 comme souhaité pour une meilleure interpolation. Le tableau de sortie temperature_cartesian est-2r par 2r ici, mais vous pouvez spécifier la taille et l'origine vous le souhaitez.

je suis venu à ce post il y a quelque temps en essayant de faire quelque chose de similaire, c'est-à-dire, en réprimandant les données polaires dans une grille cartésienne et vice-versa. La solution proposée ici fonctionne très bien. Cependant, il faut un certain temps pour effectuer la transformation des coordonnées. Je voulais juste partager une autre approche qui permet de réduire le temps de traitement jusqu'à 50 fois ou plus.

L'algorithme utilise les scipy.ndimage.interpolation.map_coordinates fonction.

voyons un petit exemple:

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)

j'espère que ce aide à quelqu'un dans la même situation que moi.

y a-t-il des paquets Python qui peuvent faire cela?

Oui! Il y a maintenant – au moins un paquet Python qui a une fonction pour re-mapper une matrice de coordonnées cartésiennes aux coordonnées polaires: abel.tools.polar.reproject_image_into_polar(), qui fait partie de la paquet PyAbel.

(Iñigo Hernáez Corres est correct,scipy.ndimage.interpolation.map_coordinates est le moyen le plus rapide que nous avons trouvé jusqu'à présent pour reprojeter de cartésiennes en coordonnées polaires.)

PyAbel peut être installé à partir de PyPi en entrant ce qui suit sur la ligne de commande:

pip install pyabel

Puis, en python, vous pouvez utiliser le code suivant pour re-projeter une image en coordonnées polaires:

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

[Selon l'application, vous pouvez envisager de passer le jacobian=True argument, qui Re-Balance les intensités de la matrice pour prendre en compte l'étirement de la grille (changer "bin sizes") qui se produit lorsque vous transformez de cartésien à polaire coordonnées.]

Voici un exemple complet:

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

Note: il y a un bon sujet de discussion (à propos de re-mappage des images en couleur en coordonnées polaires): information d'image le long d'un système de coordonnées polaires