Scan a spatio-temporal cube and run a cross-correlation on it I have a nested for-loop that looks like the following. The for-loop is used to scan a spatio-temporal (3D) cube and run a specific kind of cross-correlation used in seismology.
# some initialization
i,j,k = np.where(Pg1t &gt; 0.1)
qp = np.sum(Pg1t[:-1,20:40,20:40],axis = 1)/20

corr_arr_x = np.zeros((1750,384,384))
corr_arr_y = np.zeros((1750,384,384))

# for-loop starts
for ii in range (0,1): # this range will actually go (0,1750)
    for jj in range (31,351):
        for kk in range (31,351):
            pulse_2d = np.zeros((39,63,63))
            pulse_2d[i,j,k] = Pressure[ii + i, jj - 31 + j, kk -31 + k]
            corr_arr_x[ii,jj,kk] = np.max(signal.correlate(np.sum(pulse_2d[:,20:40,20:40],axis = 2),qp))
            corr_arr_y[ii,jj,kk] = np.max(signal.correlate(np.sum(pulse_2d[:,20:40,20:40],axis = 1),qp))

Ofcourse this is not very smart way to do it and it takes 1 minute and 36 seconds to be executed. This means, if I select the first loop to run for (0 to 1750), then it will take around 40 hours or more. So I tried the itertools.
for x in itertools.product(np.linspace(0,1,1,dtype = 'int'), np.linspace(31,351-1,351-31,dtype = 'int'), np.linspace(31,351-1,351-31,dtype = 'int')):
    pulse_2d = np.zeros((39,63,63))
    pulse_2d[i,j,k] = Pressure[x[0] + i, x[1] - 31 + j, x[2] -31 + k]
    corr_arr_x[x[0],x[1],x[2]] = np.max(signal.correlate(np.sum(pulse_2d[:,20:40,20:40],axis = 2),qp))
    corr_arr_y[x[0],x[1],x[2]] = np.max(signal.correlate(np.sum(pulse_2d[:,20:40,20:40],axis = 1),qp))

However, this takes about 1 minute 27 seconds to be executed. Not a great improvement, may save an hour at best!
Is there any way to parallelize this code or optimize in any other way so that the execution time can be significantly decreased?
Scale
Since you have
for ii in range (0,1750):
    for jj in range (31,351):
        for kk in range (31,351):
            pulse_2d = np.zeros((39,63,63))

your inner loop will be on the order of
$$
1750 \cdot (351-31)^2 \cdot 39 \cdot 63^2 = 27,738,547,200,000
$$
That is well (well) beyond the point that I would give up Python altogether and drop down to C. Even if you were to vectorize this further with Numpy as much as possible (and I'm not clear on whether that's possible), I doubt that you could get this down to a sane execution time.
Strongly consider GSL and/or BLAS.