blog_emd_option1.py

import numpy as np 
import geopandas as gpd 
import cv2 

# read in data
example = gpd.read_file("data/example.gpkg") 
scenario1 = gpd.read_file("data/example_similar.gpkg")
scenario2 = gpd.read_file("data/example_different.gpkg")

def signature_opt1(gdf, crs):
    
    centroids = gdf.geometry.to_crs(crs).centroid
    sig = np.empty((len(gdf), 3), dtype=np.float32) 
    # float32 needed as input for cv2.emd!

    # we need to normalize the data in case the total 
    # n of the two compared distributions are not equal
    sig[:,0] = gdf.example_data /
                    gdf.example_data.sum()
    sig[:,1] = centroids.x
    sig[:,2] = centroids.y
    return sig
    
sig_original = signature_opt1(example, 3035)
sig_scen1 = signature_opt1(scenario1, 3035)
sig_scen2 = signature_opt1(scenario2, 3035)

emd_scen1, _ , _ = cv2.EMD(sig_original, sig_scen1, distType = cv2.DIST_L2)

emd_scen2, _ , _ = cv2.EMD(sig_original, sig_scen2, distType = cv2.DIST_L2)

print("Earth movers distance scenario 1 (CRS: EGSG:3035): " + str(round(emd_scen1)) + " meters")
print("Earth movers distance scenario 2 (CRS: EGSG:3035): " + str(round(emd_scen2)) + " meters")