Question: So I have found similar problems in threads on here, but I haven't been able to find a solution that works for me. I am building a C extension module for Python in Visual Studio 2022 with Python 3.9. The module takes numpy arrays as inputs and returns numpy arrays. Right now, I just have it read the shape of the input array and use that to decide how big of an output array to create. If the input is shape (row, col), then the output will be shape (2, row, col). I dynamically allocate a C array with 2*row*col elements. The following is all within the function static PyObject* transform(PyObject* self, PyObject* args) {} import_array(); PyObject* data_obj; PyArrayObject* data_array; if (!PyArg_ParseTuple(args, "O", &data_obj)) { return NULL; } data_array = (PyArrayObject*)PyArray_FROM_OTF(data_obj, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY); if (data_array == NULL) { PyErr_SetString(PyExc_TypeError, "The data input must be a NumPy Array (2D or 3D)"); return NULL; } int ndim = PyArray_NDIM(data_array); npy_intp* data_shape = PyArray_SHAPE(data_array); size_t rows, columns; rows = data_shape[0]; columns = data_shape[1]; double* r_arr = (double*)malloc((rows * columns * 2) * sizeof(double)); if (r_arr == NULL) { PyErr_SetString(PyExc_ValueError, "Failed to allocate memory to arrays."); Py_DECREF(data_array); }  I then loop over rows and columns and make several calculations to generate two numbers at each row and column and store them as r_arr[rr * rows + cc] and r_arr[rr * rows + cc + rows * columns] for (size_t cc = 0; cc < columns; ++cc) { // some calculation for (size_t rr = 0; rr < rows; ++rr) { // some calculation r_arr[rr * rows + cc] = result1; r_arr[rr * rows + cc + rows * columns] = result2; PySys_WriteStdout("%f, %f\n", r_arr[rr * rows + cc], r_arr[rr * rows + cc + rows * columns]); } }  Finally, I want to turn r_arr into a numpy array to return to Python. I've tried many things here, but this is the current state it's in based on reading other threads here. npy_intp dims[3] = { 2, rows, columns }; PyObject* r_obj = (PyArrayObject *)PyArray_SimpleNewFromData(3, dims, NPY_DOUBLE, r_arr); if (r_obj == NULL) { PyErr_Print("Failed to create internal arrays. Likely due to data input being incorrect shape."); return NULL; } Py_DECREF(data_array); free(x); free(y); //free(r_arr); PyArray_ENABLEFLAGS(r_obj, NPY_ARRAY_OWNDATA); PySys_WriteStdout("Returning result\n"); return Py_BuildValue("O", r_obj);  In Python, when I create an array with shape (5, 10), and use it as an input, I get the following result. The nested loop prints out the proper values, but the array returned seems to be the kind of result I would get if I used np.empty((2, 5, 10)). I've also tried using PyArray_SimpleNewFromData() to create a 1D numpy array with 2*rows*columns elements, and that doesn't make a difference. output: 0.023581, 4.986665 0.021305, 3.986672 0.018796, 2.986680 0.016282, 1.986689 0.014207, 0.986699 -0.986851, 4.986664 -0.986798, 3.986672 -0.986746, 2.986680 -0.986700, 1.986689 -0.986667, 0.986699 -1.986746, 4.986662 -1.986718, 3.986670 -1.986691, 2.986679 -1.986666, 1.986688 -1.986648, 0.986699 -2.986705, 4.986660 -2.986685, 3.986668 -2.986665, 2.986677 -2.986646, 1.986687 -2.986632, 0.986698 -3.986679, 4.986656 -3.986662, 3.986665 -3.986645, 2.986675 -3.986629, 1.986685 -3.986616, 0.986697 -4.986658, 4.986651 -4.986642, 3.986661 -4.986626, 2.986672 -4.986611, 1.986683 -4.986599, 0.986696 -5.986637, 4.986645 -5.986622, 3.986657 -5.986607, 2.986669 -5.986592, 1.986681 -5.986579, 0.986695 -6.986616, 4.986638 -6.986601, 3.986651 -6.986586, 2.986664 -6.986571, 1.986678 -6.986558, 0.986694 -7.986593, 4.986630 -7.986578, 3.986645 -7.986563, 2.986660 -7.986547, 1.986675 -7.986533, 0.986692 -8.986567, 4.986621 -8.986552, 3.986637 -8.986536, 2.986654 -8.986520, 1.986672 -8.986505, 0.986690 Returning result array([[[ 1.14460020e-311, 1.14458244e-311, 1.14458542e-311, 1.14478895e-311, 1.14478881e-311, 1.14478881e-311, 1.14478893e-311, 1.14478894e-311, 1.14458542e-311, 1.14478894e-311], [-8.98656729e+000, 1.14478893e-311, 1.14458542e-311, 1.14478894e-311, 1.14478894e-311, -8.98655189e+000, 1.14458542e-311, 1.14478895e-311, 1.14478895e-311, 1.14478894e-311], [-8.98653594e+000, 1.14478884e-311, 1.14478894e-311, 1.14478894e-311, 1.14478894e-311, -8.98652003e+000, 1.14478884e-311, 1.14478894e-311, 1.14478894e-311, 1.14478894e-311], [-8.98650487e+000, 1.14478893e-311, 1.14478895e-311, 1.14478894e-311, 1.14478894e-311, 1.14478894e-311, 1.14478884e-311, 1.14478894e-311, 1.14478894e-311, 1.14478894e-311], [ 1.14478894e-311, 1.14478893e-311, 1.14478894e-311, 1.14478894e-311, 1.14478894e-311, 1.14478894e-311, 1.14478894e-311, 1.14478895e-311, 1.14478895e-311, 1.14478895e-311]], [[ 1.14478894e-311, 1.14458542e-311, 1.14478895e-311, 1.14478895e-311, 1.14478893e-311, 1.14478895e-311, 1.14478893e-311, 1.14458542e-311, 1.14478884e-311, 1.14478893e-311], [ 4.98662116e+000, 1.14478895e-311, 1.14478893e-311, 1.14478894e-311, 1.14478895e-311, 3.98663750e+000, 1.14478884e-311, 1.14458542e-311, 1.14478893e-311, 1.14458542e-311], [ 2.98665407e+000, 1.14478895e-311, 1.14478895e-311, 1.14478895e-311, 1.14478895e-311, 1.98667152e+000, 1.14458542e-311, 1.14458542e-311, 1.14478895e-311, 1.14478895e-311], [ 9.86690488e-001, 1.14478895e-311, 1.14478895e-311, 1.14458542e-311, 1.14478884e-311, 1.14458542e-311, 1.14478894e-311, 1.14478893e-311, 1.14458542e-311, 1.14478895e-311], [ 1.14478895e-311, 1.14478884e-311, 1.14478893e-311, 1.14478893e-311, 1.14458542e-311, 1.14478895e-311, 9.38662745e-097, 2.07712408e-308, 1.14478658e-311, 1.14458243e-311]]]) 

Answer: Here's my solution with a few other fixes. The above code in the OP does not get the data from the array properly (the calculations did not depend on the actual values of the elements, and only the shape). It also does not calculate the index properly. However the crux of the solution is that I am instantiating the array as a PyArrayObject with PyArray_ZEROS() and then getting the data from it to manipulate. PyArrayObject* data_array_obj; if (!PyArg_ParseTuple(args, "O", &data_array_obj)) { return NULL; } if (data_array_obj == NULL) { PyErr_SetString(PyExc_TypeError, "The data input must be a NumPy Array"); return NULL; } // Determine dimensionality and shape int64_t ndim = PyArray_NDIM(data_array_obj); size_t* data_shape = PyArray_SHAPE(data_array_obj); if (ndim != 2) { PyErr_SetString(PyExc_TypeError, "The input NumPy Array must be 2D"); return NULL; } size_t rows, columns; rows = data_shape[0]; columns = data_shape[1]; // Get data from input array double* data = PyArray_DATA(data_array_obj); // Instantiate output array npy_intp dims[3] = { 2, rows, columns }; PyArrayObject* r_obj_array = PyArray_ZEROS(3, dims, NPY_DOUBLE, 0); double* r_data = PyArray_DATA(r_obj_array); for (size_t cc = 0; cc < columns; ++cc) { // some calculation for (size_t rr = 0; rr < rows; ++rr) { // some calculation r_data[rr * columns + cc] = result1; r_data[rr * columns + cc + rows * columns] = result2; } } Py_DECREF(data_array); return r_obj_array;