Reproyección de rejilla polar a cartesiana

Question

Tengo una rejilla polar (r, theta) (lo que significa que cada celda es una sección de anillo) que contiene valores de alguna cantidad física (por ejemplo, temperatura) y me gustaría volver a cuadrícula (o volver a proyectar o volver a muestrear estos valores en una cuadrícula cartesiana. ¿Hay algún paquete de Python que pueda hacer esto?

No estoy interesado en convertir las coordenadas de los centros de las celdas de polar a cartesiano, esto es muy fácil. En cambio, estoy buscando un paquete que realmente pueda volver a cuadrar los datos correctamente.

¡Gracias por cualquier sugerencia!

Answer 1

Gracias por sus respuestas - después de pensar un poco más sobre esto me ocurrió con el siguiente código:

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

Lo cual no es re-grillado estrictamente, pero funciona muy bien para lo que necesito. Simplemente publicando el código en caso de que sea útil para cualquier otra persona. ¡Siéntase libre de sugerir mejoras!

Answer 2

Puede hacer esto de forma más compacta con scipy.ndimage.geometric_transform. Aquí hay un código de ejemplo:

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

Puede cambiar order=0 como se desee para una mejor interpolación. La matriz de salida temperature_cartesian es 2r por 2r aquí, pero puede especificar cualquier tamaño y origen que desee.

Answer 3

Llegué a esta publicación hace algún tiempo cuando intentaba hacer algo similar, esto es, reproyectar datos polares en una grilla cartesiana y viceversa. La solución propuesta aquí funciona bien. Sin embargo, lleva algo de tiempo realizar la transformación de coordenadas. Solo quería compartir otro enfoque que puede reducir el tiempo de procesamiento hasta 50 veces o más.

El algoritmo utiliza la función scipy.ndimage.interpolation.map_coordinates.

Veamos un pequeño ejemplo:

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) 

espero que esto ayude a alguien en la misma situación que yo.

Answer 4

Are there any Python packages that can do this?

Sí! Ahora hay, al menos, un paquete de Python que tiene una función para volver a mapear una matriz desde coordenadas cartesianas a coordenadas polares: abel.tools.polar.reproject_image_into_polar(), que es parte de PyAbel package.

(Iñigo Hernáez Corres es correcta, scipy.ndimage.interpolation.map_coordinates es la manera más rápida que hemos encontrado hasta ahora reproyectar desde cartesianas a coordenadas polares.)

PyAbel se puede instalar desde PyPi introduciendo lo siguiente en la línea de comandos:

pip install pyabel 

Luego, en Python, puede utilizar el código siguiente para volver a proyectar una imagen en coordenadas polares:

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

[Dependiendo de la aplicación, puede considerar pasar el argumento jacobian=True, que vuelve a escalar las intensidades de la matriz para incluir en la cuenta el estiramiento de la cuadrícula (cambiando el "tamaño del contenedor") que se produce cuando se transforma de cartesiano a coodinates polares.]

Aquí es un ejemplo completo:

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

Nota: hay otra buena discusión (sobre imágenes en color re-mapeo de coordenadas polares) en SO: image information along a polar coordinate system