How do I make the calculation for this distance matrix faster?

I am working on a clustering task with geospatial data. I want to compute my own distance matrix that combines both geographical and temporal distance. My data (np.array) contains latitude, longitude, and timestamp. A sample of my DataFrame df (dict to reproduce):

	    latitude	longitude	timestamp
412671	52.506136	6.068709	2017-01-01 00:00:23.518
412672	52.503316	6.071496	2017-01-01 00:01:30.764
412673	52.505122	6.068912	2017-01-01 00:02:30.858
412674	52.501792	6.068605	2017-01-01 00:03:38.194
412675	52.508105	6.075160	2017-01-01 00:06:41.116

I currently use the following code:

np_data = df.to_numpy()

# convert latitudes and longitudes to radians
lat_lon_rad = np.radians(np_data[:,:2].astype(float))

# compute Haversine distance matrix
haversine_matrix = haversine_distances(lat_lon_rad)
haversine_matrix /= np.max(haversine_matrix)

# compute time difference matrix
timestamps = np_data[:,2]
time_matrix = np.abs(np.subtract.outer(timestamps, timestamps)) # This line is SLOW
time_matrix /= np.max(time_matrix)

combined_matrix = 0.5 * haversine_matrix + 0.5 * time_matrix

This produces the desired result. However, when my data set is 1000 rows, this code takes +- 25 seconds to complete, mainly due to the calculation of the time_matrix (the haversine matrix is very fast). The problem is: I have to work with data sets of +- 200-500k rows. Using only the Haversine function is then still fine, but calculating my time_matrix will take way too long.

My question: how do I speed up the calculation of the time_matrix? I cannot find any way to perform the np.subtract.outer(timestamps, timestamps) calculation faster.

How about converting the timestamps to floats and using Numpy's broadcasting feature to compute your time differences? The integer representation of time avoids the costly overhead associated using pandas timestamps. i.e

timestamps_sec = np.array([(ts - pd.Timestamp("1970-01-01")) // pd.Timedelta("1s") for ts in np_data[:, 2]])
timestamps_sec = timestamps_sec[:, np.newaxis]

Now you can compute the time difference matrix using broadcasting

time_matrix = np.abs(timestamps_sec - timestamps_sec.T)
time_matrix = time_matrix.astype(float)  # Convert to float to avoid integer division
time_matrix /= np.max(time_matrix)