Converting coordinate tuple information to numpy arrays

Question

I have results from a finite element program that give various measures of interest (e.g., temperature, density, pressure) at regularly-spaced grid locations in three-dimensional space.

Values are equally spaced along each coordinate, but this spacing may be different for the different coordinates. For example,

x1 = [0, 0.1, 0.2, ..., 1.0]      (a total of NX1 pts) 
x2 = [0, 0.5, 1.0, ..., 20]       (a total of NX2 pts) 
x3 = [0, 0.2, 0.4, ..., 15]       (a total of NX3 pts)

The results output from the software are in the following form:

x1_1, x2_1, x3_1, f_x, g_x, h_x
x1_1, x2_1, x3_2, f_x, g_x, h_x
x1_1, x2_1, x3_3, f_x, g_x, h_x
...
x1_1, x2_2, x3_1, f_x, g_x, h_x
x1_1, x2_2, x3_2, f_x, g_x, h_x
x1_1, x2_2, x3_3, f_x, g_x, h_x
...
x1_2, x2_1, x3_1, f_x, g_x, h_x
x1_2, x2_1, x3_2, f_x, g_x, h_x
x1_2, x2_1, x3_3, f_x, g_x, h_x
...

where the f_x, g_x, h_x are the measures of interest at the particular grid point.

I would like to convert the above data format and obtain (NX1 x NX2 x NX3) numpy arrays for f, g, and h.

Some of the results sets are fairly large (80 x 120 x 100).

Does anyone have any hints for doing this conversion in an efficient manner?

Answer 1

Lets say you read your whole array into memory as an array data of shape (Nx1 * Nx2 * Nx3, 6).

data = np.loadtxt('data.txt', dtype=float, delimiter=',')

If, as your, example suggests, the points are generated in lexicographical order, you only need to grab the columns to f, g and h and reshape them:

f = data[:, 3].reshape(Nx1, Nx2, Nx3)
g = data[:, 4].reshape(Nx1, Nx2, Nx3)
h = data[:, 5].reshape(Nx1, Nx2, Nx3)

If you need to figure out what Nx1, Nx2 and Nx3 are, you can use np.unique:

Nx1 = np.unique(data[:, 0]).shape[0]
Nx2 = np.unique(data[:, 1]).shape[0]
Nx3 = np.unique(data[:, 2]).shape[0]

A more robust solution in case the order of the points is not guaranteed, would be to use np.unique to extract indices to the grid values:

Nx1, idx1 = np.unique(data[:, 0], return_inverse=True)
Nx1 = Nx1.shape[0]
Nx2, idx2 = np.unique(data[:, 1], return_inverse=True)
Nx2 = Nx2.shape[1]
Nx3, idx3 = np.unique(data[:, 2], return_inverse=True)
Nx3 = Nx3.shape[0]

f = np.empty((Nx1, Nx2, Nx3))
f[idx1, idx2, idx3] = data[:, 3]
g = np.empty((Nx1, Nx2, Nx3))
g[idx1, idx2, idx3] = data[:, 4]
h = np.empty((Nx1, Nx2, Nx3))
h[idx1, idx2, idx3] = data[:, 5]

This will create new arrays for f, g and h, not views into the original data array, so it will use more memory.

And of course instead of my ugly code above repeating everything three times, you should use a loop, or a list comprehension!