Riproiezione polare su griglia cartesiana

Question

Ho una griglia polare (r, theta) (che significa che ogni cella è una sezione di anulus) contenente valori di una certa quantità fisica (ad es. Temperatura), e vorrei ricollocare (o riprogettare o ricampionare questi valori su una griglia cartesiana. Ci sono pacchetti Python che possono farlo?

Non sono interessato a convertire le coordinate dei centri delle celle da polari a cartesiani - questo è molto semplice. Invece, sto cercando un pacchetto che possa effettivamente ricollocare correttamente i dati.

Grazie per eventuali suggerimenti!

Answer 1

Grazie per le vostre risposte - dopo averci pensato un po 'di più su questo mi si avvicinò con il seguente codice:

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

che non è strettamente ri-gridding, ma funziona bene per quello che mi serve. Inserendo semplicemente il codice nel caso in cui sia utile a chiunque altro. Sentiti libero di suggerire miglioramenti!

Answer 2

È possibile eseguire questa operazione in modo più compatto con scipy.ndimage.geometric_transform. Ecco alcuni esempi di codice:

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

È possibile modificare order=0 lo desideri per una migliore interpolazione. L'array di output temperature_cartesian è 2r di 2r qui, ma è possibile specificare qualsiasi dimensione e origine che si desidera.

Answer 3

Sono arrivato a questo post qualche tempo fa quando cercavo di fare qualcosa di simile, questo è, riproiettando i dati polari in una griglia cartesiana e viceversa. La soluzione qui proposta funziona bene. Tuttavia, ci vuole del tempo per eseguire la trasformazione delle coordinate. Volevo solo condividere un altro approccio che può ridurre i tempi di elaborazione fino a 50 volte o più.

L'algoritmo utilizza la funzione scipy.ndimage.interpolation.map_coordinates.

Vediamo un piccolo esempio:

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) 

Spero che questo aiuta qualcuno nella mia stessa situazione.

Answer 4

Are there any Python packages that can do this?

Sì! C'è ora - almeno un - pacchetto Python che ha una funzione per ri-mappare una matrice da coordinate cartesiane a coordinate polari: abel.tools.polar.reproject_image_into_polar(), che fa parte dello PyAbel package.

(Iñigo Hernáez Corres è corretta, scipy.ndimage.interpolation.map_coordinates è il modo più veloce che abbiamo trovato finora riproiettare da cartesiane a coordinate polari.)

PyAbel può essere installato da PyPi inserendo la seguente nella riga di comando:

pip install pyabel 

Poi, in python, è possibile utilizzare il seguente codice di ri-progetto di un'immagine in coordinate polari:

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

[A seconda dell'applicazione, si potrebbe considerare di passare l'argomento jacobian=True, che ridimensiona le intensità della matrice per prendere in considerazione l'allungamento della griglia (cambiando "dimensioni bin") che si verifica quando si trasforma da cartesiano ai coodati polari.]

Ecco un esempio 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: c'è un'altra buona discussione (su immagini a colori ri-mappatura in coordinate polari) su SO: image information along a polar coordinate system