Rearrange data in two-dimensional array according to transformation from polar to Cartesian coordinates

Question

I have a two-dimensional array that represents function values at positions in a polar coordinate system. For example:

import numpy as np

radius = np.linspace(0, 1, 50)
angle = np.linspace(0, 2*np.pi, radius.size)
r_grid, a_grid = np.meshgrid(radius, angle)
data = np.sqrt((r_grid/radius.max())**2
               + (a_grid/angle.max())**2)

Here the data is arranged in a rectangular grid corresponding to the polar coordinates. I want to rearrange the data in the array such that the axes represent the corresponding Cartesian coordinate system. The old versus new layout can be visualized as follows:

import matplotlib.pyplot as plt

fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=plt.figaspect(0.5))
ax1.set(title='Polar coordinates', xlabel='Radius', ylabel='Angle')
ax1.pcolormesh(r_grid, a_grid, data)
ax2.set(title='Cartesian coordinates', xlabel='X', ylabel='Y')
x_grid = r_grid * np.cos(a_grid)
y_grid = r_grid * np.sin(a_grid)
ax2.pcolormesh(x_grid, y_grid, data)

Here the coordinates are explicitly given and the plot is adjusted accordingly. I want the data to be rearranged in the data array itself instead. It should contain all values, optionally filling with zeros to fit the shape (similar to scipy.ndimage.rotate(..., reshape=True)).

If I manually loop over the polar arrays to compute the Cartesian coordinates, the result contains empty regions which ideally should be filled as well:

new = np.zeros_like(data)
visits = np.zeros_like(new)
for r, a, d in np.nditer((r_grid, a_grid, data)):
    i = 0.5 * (1 + r * np.sin(a)) * new.shape[0]
    j = 0.5 * (1 + r * np.cos(a)) * new.shape[1]
    i = min(int(i), new.shape[0] - 1)
    j = min(int(j), new.shape[1] - 1)
    new[i, j] += d
    visits[i, j] += 1
new /= np.maximum(visits, 1)
ax2.imshow(new, origin='lower')

Is there a way to achieve the transformation while avoiding empty regions in the resulting data array?

Answer 1

tl;dr: No, not without changing some conditions of your problem.

The artefact you are seeing is a property of the transformation. It is not due to the fixed resolution in angle for all radii. Hence, it is not due to a wrong or bad implementation of the transformation. The Cartesian grid simply implies a higher special resolution at these areas as there are resolved points from the polar map.

• The only "clean" way (that I can think of right now) to handle this is to have an adjustable resolution in the polar coordinates to account for the 1/r scaling. (If you input data allows it)

• A somewhat cheating way of visualizing this without the gaps would to randomly distribute them over the gaps. The argument here is, that you do not have the resolution to decide in which bin they were to begin with. So you could just randomly throw them in one which might have been a possible origin and not throw them all in same one (as you are doing right now). However, I would like discourage this stronlgy. It just gives you a prettier plot. Note, that this is somewhat equivalent to the behaviour of the upper right plot in your question.

Answer 2

This doesn't really give the expected result, but maybe will help you in achieving a solution after some needed correction...

import numpy as np

radius = np.linspace(0, 1, 50)
angle = np.linspace(0, 2*np.pi, radius.size)
r_grid, a_grid = np.meshgrid(radius, angle)
data = np.sqrt((r_grid/radius.max())**2
               + (a_grid/angle.max())**2)


def polar_to_cartesian(data):
    new = np.zeros_like(data) * np.nan
    x = np.linspace(-1, 1, new.shape[1])
    y = np.linspace(-1, 1, new.shape[0])
    for i in range(new.shape[0]):
        for j in range(new.shape[1]):
            x0, y0 = x[j], y[i]
            r, a = np.sqrt(x0**2 + y0**2), np.arctan2(y0, x0)
            data_i = np.argmin(np.abs(a_grid[:, 0] - a))
            data_j = np.argmin(np.abs(r_grid[0, :] - r))
            val = data[data_i, data_j]

            if r <= 1:
                new[i, j] = val

    return new

new = polar_to_cartesian(data)
fig, ax = plt.subplots()
ax.imshow(new, origin='lower')

EDIT: Modified using np.arctan2 according to the suggestions of OP.