Repository: geofileops/geofileops
Pull request #634 (2025-02-17 15:02:38)

BUG: When calculating join-nearest, the distance returned is not correct for polygons
When calculating join-nearest and with spatialite verion >= 5.1 the distance calculated before this PR was not correct for polygons. It was the distance between one polygon and the centroid of the other.

codecov[bot]: ## [Codecov](https://app.codecov.io/gh/geofileops/geofileops/pull/634?dropdown=coverage&src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=geofileops) Report
All modified and coverable lines are covered by tests :white_check_mark:
> Project coverage is 86.27%. Comparing base [(`5f08842`)](https://app.codecov.io/gh/geofileops/geofileops/commit/5f08842e475a5aaff2f4ace75e7bc77c9c2fe023?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=geofileops) to head [(`5205d2a`)](https://app.codecov.io/gh/geofileops/geofileops/commit/5205d2a2cf0836470270c875e71f08bb7a23c6a7?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=geofileops).
> Report is 109 commits behind head on main.

<details><summary>Additional details and impacted files</summary>


```diff
@@            Coverage Diff             @@
##             main     #634      +/-   ##
==========================================
- Coverage   93.11%   86.27%   -6.85%     
==========================================
  Files          21       22       +1     
  Lines        4113     4489     +376     
==========================================
+ Hits         3830     3873      +43     
- Misses        283      616     +333     
```

</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/geofileops/geofileops/pull/634?dropdown=coverage&src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=geofileops).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=geofileops).


Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index e866fbb..24a5aaf 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -29,6 +29,9 @@
 - Don't convert to multi-part geometries by default in `copy_layer`,... (#570)
 - Add configuration option to only warn on dissolve errors (#561)
 - Add some pre-flight checks when geofileops is imported (#573, #627)
+- When using join_nearest with spatialite version >= 5.1,
+  show ST_distance between the two geometries instead of 
+  the distance between the centroid of the two geometries (#634)
 
 ### Bugs fixed
 
diff --git a/geofileops/util/_geoops_sql.py b/geofileops/util/_geoops_sql.py
index 9171590..7cf8953 100644
--- a/geofileops/util/_geoops_sql.py
+++ b/geofileops/util/_geoops_sql.py
@@ -2351,7 +2351,11 @@ def join_nearest(
             SELECT layer1.{{input1_geometrycolumn}} as geom
                   {{layer1_columns_prefix_alias_str}}
                   {{layer2_columns_prefix_alias_str}}
-                  ,k.pos, k.distance_m AS distance, k.distance_crs
+                  ,k.pos
+                  ,ST_Distance(
+                    layer1.{{input1_geometrycolumn}}, layer2.{{input2_geometrycolumn}}
+                  ) AS distance
+                  ,k.distance_crs
               FROM "{{input1_layer}}" layer1
               JOIN knn2 k
               JOIN "{{input2_layer}}" layer2 ON layer2.rowid = k.fid


Test patch:
diff --git a/tests/test_geofileops_twolayers.py b/tests/test_geofileops_twolayers.py
index a1c3791..d9cc6b4 100644
--- a/tests/test_geofileops_twolayers.py
+++ b/tests/test_geofileops_twolayers.py
@@ -1546,6 +1546,39 @@ def test_join_nearest_invalid_params(
         )
 
 
+def test_join_nearest_distance(tmp_path):
+    geoms = [
+        "POLYGON ((0 0, 3 0, 3 3, 0 3, 0 0))",
+        "POLYGON ((10 1, 13 1, 13 4, 10 4, 10 1))",
+    ]
+    for index, geom in enumerate(geoms):
+        box_geom = shapely.from_wkt(geom)
+        box_geom = shapely.segmentize(box_geom, 1)
+        gdf_geom = gpd.GeoDataFrame(geometry=[box_geom], crs="EPSG:31370")
+        geom_path = tmp_path / f"geom{index+1}.gpkg"
+        gfo.to_file(gdf=gdf_geom, path=geom_path)
+
+    # Test
+    output_path = tmp_path / "geom_join_nearest.gpkg"
+    gfo.join_nearest(
+        input1_path=tmp_path / "geom1.gpkg",
+        input2_path=tmp_path / "geom2.gpkg",
+        output_path=output_path,
+        nb_nearest=1,
+        distance=50,
+        expand=True,
+        force=True,
+    )
+
+    # Check if the output file is correctly created
+    assert output_path.exists()
+    output_layerinfo = gfo.get_layerinfo(output_path)
+    assert output_layerinfo.featurecount == 1
+    # Check the contents of the result file
+    output_gdf = gfo.read_file(output_path)
+    assert output_gdf["distance"][0] == 7
+
+
 @pytest.mark.parametrize(
     "suffix, epsg, gridsize",
     [(".gpkg", 31370, 0.01), (".gpkg", 4326, 0.0), (".shp", 31370, 0.01)],