Skip to content

Instantly share code, notes, and snippets.

@AlexandraKapp
Created December 1, 2021 19:29
Show Gist options
  • Save AlexandraKapp/c8fac1fefff933c8dc1452c245606b2f to your computer and use it in GitHub Desktop.
Save AlexandraKapp/c8fac1fefff933c8dc1452c245606b2f to your computer and use it in GitHub Desktop.
from haversine import haversine, Unit
# construct the cost matrix
def get_cost_matrix(gdf1, gdf2):
gdf1_centroids = gdf1.to_crs(3395).centroid.to_crs(4326)
gdf2_centroids = gdf2.to_crs(3395).centroid.to_crs(4326)
coords_sig1 = list(zip(gdf1_centroids.y, gdf1_centroids.x))
coords_sig2 = list(zip(gdf2_centroids.y, gdf2_centroids.x))
#get all potential combinations between all points from sig1 and sig2
grid = np.meshgrid(range(0, len(coords_sig1)),
range(0, len(coords_sig2)))
tile_combinations = np.array([grid[0].flatten(), grid[1].flatten()])
# create an empty cost matrix with the length of all possible combinations
cost_matrix = np.empty(tile_combinations.shape[1], dtype = np.float32) # float32 needed as input for cv2.emd!
# compute haversine distance for all possible combinations
for column in range(0, tile_combinations.shape[1]):
tile_1 = tile_combinations[0, column]
tile_2 = tile_combinations[1, column]
cost_matrix[column] = haversine(coords_sig1[tile_1], coords_sig2[tile_2], unit = Unit.METERS)
# reshape array to matrix
return np.reshape(cost_matrix, (len(coords_sig1),len(coords_sig2)))
# as the coordinates are the same in both scenarios
# you could use the same cost matrix for both scenarios
cost_matrix1 = get_cost_matrix(example, scenario1)
cost_matrix2 = get_cost_matrix(example, scenario2)
# The signature for the custom cost matrix does not need coordinates. It can only consist of an array of values
sig_original = np.array(example.example_data, dtype = np.float32)
sig_scen1 = np.array(scenario1.example_data, dtype = np.float32)
# Compute the EMD
emd_scen1, _ , _ = cv2.EMD(sig_original, sig_scen1, distType = cv2.DIST_USER, cost = cost_matrix1)
emd_scen2, _ , _ = cv2.EMD(sig_original, sig_scen2, distType = cv2.DIST_USER, cost = cost_matrix2)
print("Earth movers distance scenario 1: " + str(round(emd_scen1)) + " meters")
print("Earth movers distance scenario 2: " + str(round(emd_scen2)) + " meters")
sig_scen2 = np.array(scenario2.example_data, dtype = np.float32)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment