Intercept Coincident points

libpysal/graph/_triangulation.py

@@ -83,6 +83,28 @@ def _delaunay(coordinates, ids=None, bandwidth=numpy.inf, kernel="boxcar"):
         coordinates, ids=ids, valid_geometry_types=_VALID_GEOMETRY_TYPES
     )
 
+    n_coincident, coincident_lut = _validate_coincident(geoms)
+    if n_coincident > 0:
+        if coincident == "raise":
+            raise ValueError(
+                f"There are {(coincident_lut[geoms.index.name].str.len()>1).sum()}"
+                f"unique locations in the dataset, but {len(geoms)} observations."
+                "This means there are multiple points in the same location, which"
+                " is undefined for this graph type. To address this issue, consider setting "
+                " `coincident='clique' or consult the documentation about coincident points."
+            )
+        elif coincident == "jitter":
+            coordinates, geoms = _jitter_geoms(coordinates, geoms)
+        elif coincident == "clique":
+            input_coordinates, input_ids, input_geoms = coordinates, ids, geoms
+            coordinates, ids, geoms = _validate_geometry_input(
+                coincident_lut.geometry, ids=coincident_lut.index, valid_geometry_types=_VALID_GEOMETRY_TYPES
+            )
+        else:
+            raise ValueError(
+                f"Recieved option `coincident='{coincident}', but only options 'raise','clique','jitter' are suppported."
+            )
+
     edges, _ = _voronoi_edges(coordinates)
 
     ids = numpy.asarray(ids)
@@ -100,6 +122,9 @@ def _delaunay(coordinates, ids=None, bandwidth=numpy.inf, kernel="boxcar"):
     # dropped points from the triangulation and the
     # misalignment of the weights and the attribute array
 
+    if (n_coincident > 0) & (coincident == "clique"):
+       graph = _induce_cliques(graph, coincident_lut, fill_value=1)
+
     return head, tail, weights
 
 
@@ -151,6 +176,27 @@ def _gabriel(coordinates, ids=None, bandwidth=numpy.inf, kernel="boxcar"):
     coordinates, ids, _ = _validate_geometry_input(
         coordinates, ids=ids, valid_geometry_types=_VALID_GEOMETRY_TYPES
     )
+    n_coincident, coincident_lut = _validate_coincident(geoms)
+    if n_coincident > 0:
+        if coincident == "raise":
+            raise ValueError(
+                f"There are {(coincident_lut[geoms.index.name].str.len()>1).sum()}"
+                f"unique locations in the dataset, but {len(geoms)} observations."
+                "This means there are multiple points in the same location, which"
+                " is undefined for this graph type. To address this issue, consider setting "
+                " `coincident='clique' or consult the documentation about coincident points."
+            )
+        elif coincident == "jitter":
+            coordinates, geoms = _jitter_geoms(coordinates, geoms)
+        elif coincident == "clique":
+            input_coordinates, input_ids, input_geoms = coordinates, ids, geoms
+            coordinates, ids, geoms = _validate_geometry_input(
+                coincident_lut.geometry, ids=coincident_lut.index, valid_geometry_types=_VALID_GEOMETRY_TYPES
+            )
+        else:
+            raise ValueError(
+                f"Recieved option `coincident='{coincident}', but only options 'raise','clique','jitter' are suppported."
+            )
 
     edges, dt = _voronoi_edges(coordinates)
     droplist = _filter_gabriel(
@@ -222,6 +268,27 @@ def _relative_neighborhood(coordinates, ids=None, bandwidth=numpy.inf, kernel="b
     coordinates, ids, _ = _validate_geometry_input(
         coordinates, ids=ids, valid_geometry_types=_VALID_GEOMETRY_TYPES
     )
+    n_coincident, coincident_lut = _validate_coincident(geoms)
+    if n_coincident > 0:
+        if coincident == "raise":
+            raise ValueError(
+                f"There are {(coincident_lut[geoms.index.name].str.len()>1).sum()}"
+                f"unique locations in the dataset, but {len(geoms)} observations."
+                "This means there are multiple points in the same location, which"
+                " is undefined for this graph type. To address this issue, consider setting "
+                " `coincident='clique' or consult the documentation about coincident points."
+            )
+        elif coincident == "jitter":
+            coordinates, geoms = _jitter_geoms(coordinates, geoms)
+        elif coincident == "clique":
+            input_coordinates, input_ids, input_geoms = coordinates, ids, geoms
+            coordinates, ids, geoms = _validate_geometry_input(
+                coincident_lut.geometry, ids=coincident_lut.index, valid_geometry_types=_VALID_GEOMETRY_TYPES
+            )
+        else:
+            raise ValueError(
+                f"Recieved option `coincident='{coincident}', but only options 'raise','clique','jitter' are suppported."
+            )
 
     edges, dt = _voronoi_edges(coordinates)
     output, _ = _filter_relativehood(edges, dt.points, return_dkmax=False)
@@ -293,9 +360,36 @@ def _voronoi(coordinates, ids=None, clip="extent", rook=True):
     coordinates, ids, _ = _validate_geometry_input(
         coordinates, ids=ids, valid_geometry_types=_VALID_GEOMETRY_TYPES
     )
+    n_coincident, coincident_lut = _validate_coincident(geoms)
+    if n_coincident > 0:
+        if coincident == "raise":
+            raise ValueError(
+                f"There are {(coincident_lut[geoms.index.name].str.len()>1).sum()}"
+                f"unique locations in the dataset, but {len(geoms)} observations."
+                "This means there are multiple points in the same location, which"
+                " is undefined for this graph type. To address this issue, consider setting "
+                " `coincident='clique' or consult the documentation about coincident points."
+            )
+        elif coincident == "jitter":
+            coordinates, geoms = _jitter_geoms(coordinates, geoms)
+        elif coincident == "clique":
+            input_coordinates, input_ids, input_geoms = coordinates, ids, geoms
+            coordinates, ids, geoms = _validate_geometry_input(
+                coincident_lut.geometry, ids=coincident_lut.index, valid_geometry_types=_VALID_GEOMETRY_TYPES
+            )
+        else:
+            raise ValueError(
+                f"Recieved option `coincident='{coincident}', but only options 'raise','clique','jitter' are suppported."
+            )
 
     cells, _ = voronoi_frames(coordinates, clip=clip)
-    return _vertex_set_intersection(cells, rook=rook, ids=ids)
+    graph = _vertex_set_intersection(cells, rook=rook, ids=ids)
+
+    if (n_coincident > 0) & (coincident == "clique"):
+        fill_value = _kernel_functions[kernel](0, bandwidth)
+        new_adjtable = _induce_cliques(graph, coincident_lut, fill_value=fill_value)
+        # from here, how to ensure ordering?
+    return graph
 
 
 #### utilities

libpysal/graph/_utils.py

@@ -3,6 +3,47 @@
 import pandas as pd
 import shapely
 
+def _jitter_geoms(*args,**kwargs):
+    raise NotImplementedError()
+
+def _induce_cliques(graph, clique_to_members, fill_value=1):
+    """
+    HERE BE DRAGONS! the "members" array must be in the same order as
+    the intended output graph. 
+    """
+    across_clique = graph.adjacency.merge(
+        clique_to_members['index'], left_index=True, right_index=True
+    ).explode('index').rename(
+        columns=dict(index='subclique_focal')
+    ).merge(
+        clique_to_members['index'], left_on='neighbor', right_index=True
+    ).explode('index').rename(
+        columns=dict(index='subclique_neighbor')
+    ).reset_index().drop(
+        ['index','neighbor'], axis=1
+    ).rename(
+        columns=dict(subclique_focal="focal", subclique_neighbor='neighbor')
+    )
+    is_clique = lut['index'].str.len()>1
+    within_clique = lut[
+        is_clique
+    ]['index'].apply(
+        lambda x: list(combinations(x, 2))
+    ).explode().apply(
+        pandas.Series
+    ).rename(columns={0:"focal", 1:"neighbor"}).assign(weight=fill_value)
+
+    new_adj = pandas.concat((across_clique, within_clique), ignore_index=True, axis=0).reset_index(drop=True)
+
+    return new_adj
+
+
+
+    # try building on top of earlier blockweights code for cliques
+    # then use index tricks iterating over cliques
+
+
+
 
 def _neighbor_dict_to_edges(neighbors, weights=None):
     """
@@ -22,6 +63,24 @@ def _neighbor_dict_to_edges(neighbors, weights=None):
         data_array[heads == tails] = 0
     return heads, tails, data_array
 
+def _validate_coincident(geoms):
+    """
+    Identify coincident points and create a look-up table for the coincident geometries. 
+    """
+    valid_coincident_geom_types = set(("Point",))
+    if not set(geoms.geom_type) <= valid_coincident_geom_types:
+        raise ValueError(
+            f"coindicence checks are only well-defined for geom_types: {valid_coincident_geom_types}"
+            )
+    wkb = geoms.to_frame('geometry').assign(hash = geoms.apply(lambda x: x.wkb)).reset_index(names='index')
+    if len(wkb) > len(wkb.hash.unique()): 
+        grouper = wkb.groupby("hash")
+        ids = grouper['index'].agg(list)
+        max_coincident = ids.str.len().max()
+        unique_locs = grouper['geometry'].agg("first")
+        return max_coincident, geopandas.GeoDataFrame(ids.to_frame("index"), geometry=unique_locs)
+    else:
+        return 0, wkb
 
 def _validate_geometry_input(geoms, ids=None, valid_geometry_types=None):
     """