I’ve prepared a code to make it work in Python. But it’s taking too long to finish. I was using similar version of this in MATLAB, but it wasn’t that slow.

import pandas as pd
import numpy as np

# Load the CSV file
csv_file = 'C:/Users/nad/Downloads/Yard/070/interpolated_points.csv'
df = pd.read_csv(csv_file)

# Extract the coordinates (assumed to be in meters)
X_coords = df['X'].values
Z_coords = df['Z'].values
Y_coords = df['Y'].values

# Define the cell size in meters
cell_size = 0.0048  # Cell size in meters

# Define the cube dimensions in terms of number of cells
depth_cells = 104
width_cells = depth_cells * 2
length_cells = depth_cells * 2

# Convert cube dimensions to real-world units (meters)
depth = depth_cells * cell_size
width = width_cells * cell_size
length = length_cells * cell_size

# Determine the bounds of the entire grid
min_x, max_x = X_coords.min() - width / 2, X_coords.max() + width / 2
min_z, max_z = Z_coords.min(), Z_coords.max() + depth
min_y, max_y = Y_coords.min() - length / 2, Y_coords.max() + length / 2

# Convert to integers
grid_size_x = int((max_x - min_x) / cell_size)
grid_size_z = int((max_z - min_z) / cell_size)
grid_size_y = int((max_y - min_y) / cell_size)

# Create the 3D grid with real-world coordinates for each cell
x_coords_grid = np.linspace(min_x, max_x, grid_size_x)
z_coords_grid = np.linspace(min_z, max_z, grid_size_z)
y_coords_grid = np.linspace(min_y, max_y, grid_size_y)

# Create 3D meshgrids for the x, y, and z coordinates of the grid
X_grid, Z_grid, Y_grid = np.meshgrid(x_coords_grid, z_coords_grid, y_coords_grid, indexing='ij')

# Load your GPR data matrix (assumed to have the same number of columns as there are points in the CSV)
gpr_data = np.loadtxt('C:/Users/nad/Downloads/Yard/070/gpr_data.csv', delimiter=',')

# Define the distance interval for GPR data, using the updated cell size as resolution
gpr_data_resolution = cell_size  # Set to cell size

# Initialize an empty 3D array to store the summed GPR data for each grid cell
gpr_grid_data = np.zeros((grid_size_x, grid_size_z, grid_size_y))

# Vectorized distance calculation for each point
for point_index, (x_point, z_point, y_point) in enumerate(zip(X_coords, Z_coords, Y_coords)):
    # Calculate the Euclidean distance from all grid cells to the current point
    distance_grid = 2 * np.sqrt((X_grid - x_point) ** 2 +
                            (Z_grid - z_point) ** 2 +
                            (Y_grid - y_point) ** 2)
    
    # Determine the row in the GPR data that corresponds to this distance
    gpr_data_rows = (distance_grid / gpr_data_resolution).astype(int)

    # Mask to ensure valid row indices
    valid_mask = (gpr_data_rows >= 0) & (gpr_data_rows < gpr_data.shape[0])
    
    # Pull the GPR data for the valid indices using gpr_data_rows, and set invalid to 0
    gpr_data_values = np.zeros_like(distance_grid)
    
    # Extract the GPR data for valid indices
    gpr_data_values[valid_mask] = gpr_data[gpr_data_rows[valid_mask], point_index]

    # Accumulate the GPR data for this point into the final 3D grid
    gpr_grid_data += gpr_data_values

# Now, gpr_grid_data contains the GPR data for each cell after summing all points' contributions

I tried vectorizing all and avoiding using the for loop, then I got memory issue. What can I do for it to make it much much faster? Should I not use Python?