Reprojecting polar to cartesian grid

Question

I have a polar (r,theta) grid (which means that each cell is an annulus section) containing values of some physical quantity (e.g. temperature), and I would like to re-grid (or re-project, or resample) these values onto a cartesian grid. Are there any Python packages that can do this?

I am not interested in converting the coordinates of the centers of the cells from polar to cartesian - this is very easy. Instead, I'm looking for a package that can actually re-grid the data properly.

Thanks for any suggestions!

Recidivism

That's not an easy problem, and it would be both interesting and a huge bear to write. I think it would take me 2-3 days to come up with something horribly inefficient. – Avila

Answer 1

Thanks for your answers - after thinking a bit more about this I came up with the following code:

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

Which is not strictly re-gridding, but works fine for what I need. Just posting the code in case it is useful to anyone else. Feel free to suggest improvements!

Recidivism

Just a tiny correction. I guess, it should be arctan2(Y, X) in your code. – Miso

Answer 2

I came to this post some time ago when trying to do something similar, this is, reprojecting polar data into a cartesian grid and vice-versa. The solution proposed here works fine. However, it takes some time to perform the coordinate transform. I just wanted to share another approach which can reduce the processing time up to 50 times or more.

The algorithm uses the scipy.ndimage.interpolation.map_coordinates function.

Let's see a little example:

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)

I hope this helps someone in the same situation as me.

Chaim

great! ...minor correction is coords should be constructed from Tc and Rc, not Ac and Sc – Expletive

Answer 3

OpenCV 3.4 can do this now pretty easily with warpPolar()

Very simple to call:

import numpy as np
import cv2
from matplotlib import pyplot as plt

# Read in our image from disk
image = cv2.imread('washington_quarter.png',0)
plt.imshow(image),plt.show()

margin = 0.9 # Cut off the outer 10% of the image
# Do the polar rotation along 1024 angular steps with a radius of 256 pixels.
polar_img = cv2.warpPolar(image, (256, 1024), (image.shape[0]/2,image.shape[1]/2), image.shape[1]*margin*0.5, cv2.WARP_POLAR_LINEAR)
# Rotate it sideways to be more visually pleasing
polar_img = cv2.rotate(polar_img, cv2.ROTATE_90_COUNTERCLOCKWISE)
plt.imshow(polar_img),plt.show()

Chaucerian

Answer 4

You can do this more compactly with scipy.ndimage.geometric_transform. Here is some sample 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])})

You can change order=0 as desired for better interpolation. The output array temperature_cartesian is 2r by 2r here, but you can specify any size and origin you like.

Bosket

Answer 5

Are there any Python packages that can do this?

Yes! There is now – at least one – Python package that has a function to re-map a matrix from cartesian to polar coordinates: abel.tools.polar.reproject_image_into_polar(), which is part of the PyAbel package.

(Iñigo Hernáez Corres is correct, scipy.ndimage.interpolation.map_coordinates is the fastest way that we have found so far to reproject from cartesian to polar coordinates.)

PyAbel can be installed from PyPi by entering the following on the command line:

pip install pyabel

Then, in python, you can use the following code to re-project an image into polar coordinates:

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

[Depending on the application, you might consider passing the jacobian=True argument, which re-scales the intensities of the matrix to take into the account the stretching of the grid (changing "bin sizes") that takes place when you transform from Cartesian to polar coodinates.]

Here is a complete example:

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: there is another good discussion (about re-mapping color images to polar coordinates) on SO: image information along a polar coordinate system

Theadora

This is a nice example. However it gives me TypeError: 'numpy.float64' object cannot be interpreted as an integer on python3.4. If you're the maintainer of the code you should check that out. – Timmy

It seems your code transforms images on an index-grid, where the origin is always on a pixel. Sometimes, the origin will need to be in between adjacent pixels. The circle in your example becomes a non-straight line after transformation, partially because the circle on the left isn't centered in the image, partially b/c the origin is shifted half a pixel. You may experiment with the same data which puts the circle in the center abel.tools.analytical.SampleImage(501).image[200:-200, 200:-200] and then interpolated it to an 100x100 grid (was 101x101) and then do your transformation. – Unravel

Answer 6

Here is my take which, to me, is the most straightforward:

from scipy.interpolate import LinearNDInterpolator
import numpy as np
import plotly.graph_objects as go

# List of thetas and rs for which your data is defined
t = np.linspace(0, 2*np.pi, 100)
r = np.linspace(0, 10, 50)

# Create a grid from those values
T, R = np.meshgrid(t, r)

# Convert to cartesian
xs = (R*np.cos(T)).flatten()
ys = (R*np.sin(T)).flatten()

# Define some arbitrary data at the (converted) x, y values 
zs = np.abs(xs + ys)

# Define a new cartesian grid 
X = np.linspace(-10, 10, 100)
Y = np.linspace(-10, 10, 100)
X, Y = np.meshgrid(X, Y) 

# Create the interpolating function
interp = LinearNDInterpolator(list(zip(xs, ys)), zs)
Z = interp(X, Y)
# Get rid of nasty values
Z[np.isnan(Z)] = 0


# Plot the data
fig = go.Figure([
    go.Scatter(x=xs, y=ys, mode='markers', marker=dict(size=zs)),
    go.Scatter(x=X.flatten(), y=Y.flatten(), mode='markers', marker=dict(size=Z.flatten()))
    
    ]
    )
fig.update_layout(showlegend=True, height=1000, width=1000)
fig.show()

The size of the markers reflects the quantity of interest

Mainly