Skip to content

Instantly share code, notes, and snippets.

@AlexandraKapp
Created December 1, 2021 19:25
Show Gist options
  • Save AlexandraKapp/20c51dd788985fd8c059b507cf23875e to your computer and use it in GitHub Desktop.
Save AlexandraKapp/20c51dd788985fd8c059b507cf23875e to your computer and use it in GitHub Desktop.
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")
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment