Repository: geofileops/geofileops
Pull request #570 (2024-10-22 17:57:10)

ENH: don't write multi geometries by default in copy_layer,...
The reason was that "PROMOTE_TO_MULTI" was always applied if `explodecollection` was `False`.

But apparently a side effect of a `geometrytype` being specified lead to arrow not being used in `gdal.VectorTranslate`... hence removing it lead to arrow being used. But, apparently datetime columns weren't always handled correctly for for e.g. Geopackages by arrow... (https://github.com/OSGeo/gdal/pull/11213).

So, for gdal versions < GDAL 3.11, arrow is now disabled using config option `OGR2OGR_USE_ARROW_API NO`.

codecov[bot]: ## [Codecov](https://app.codecov.io/gh/geofileops/geofileops/pull/570?dropdown=coverage&src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=geofileops) Report
Attention: Patch coverage is `42.85714%` with `8 lines` in your changes missing coverage. Please review.
> Project coverage is 90.52%. 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 [(`07023fa`)](https://app.codecov.io/gh/geofileops/geofileops/commit/07023fad0d731d65de2b76b1647cba344ef10609?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=geofileops).
> Report is 62 commits behind head on main.

| [Files with missing lines](https://app.codecov.io/gh/geofileops/geofileops/pull/570?dropdown=coverage&src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=geofileops) | Patch % | Lines |
|---|---|---|
| [geofileops/util/\_ogr\_util.py](https://app.codecov.io/gh/geofileops/geofileops/pull/570?src=pr&el=tree&filepath=geofileops%2Futil%2F_ogr_util.py&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=geofileops#diff-Z2VvZmlsZW9wcy91dGlsL19vZ3JfdXRpbC5weQ==) | 27.27% | [8 Missing :warning: ](https://app.codecov.io/gh/geofileops/geofileops/pull/570?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=geofileops) |

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


```diff
@@            Coverage Diff             @@
##             main     #570      +/-   ##
==========================================
- Coverage   93.11%   90.52%   -2.60%     
==========================================
  Files          21       21              
  Lines        4113     4241     +128     
==========================================
+ Hits         3830     3839       +9     
- Misses        283      402     +119     
```

</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/geofileops/geofileops/pull/570?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 948e40e..7a86d9f 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -9,6 +9,7 @@
   column should be retained (#523)
 - Enable "CURVE" geometrytype files to be processed in the general file and
   layer operations (#558)
+- Don't convert to multi-part geometries by default in `copy_layer`,... (#570)
 - Add configuration option to only warn on dissolve errors (#561)
 
 ## 0.9.1 (2024-07-18)
diff --git a/geofileops/_compat.py b/geofileops/_compat.py
index c4f2726..db631f7 100644
--- a/geofileops/_compat.py
+++ b/geofileops/_compat.py
@@ -18,3 +18,4 @@ spatialite_version = spatialite_version_info["spatialite_version"]
 SPATIALITE_GTE_51 = version.parse(spatialite_version) >= version.parse("5.1")
 
 GDAL_GTE_38 = version.parse(gdal.__version__) >= version.parse("3.8")
+GDAL_STE_310 = version.parse(gdal.__version__) <= version.parse("3.10")
diff --git a/geofileops/fileops.py b/geofileops/fileops.py
index b150a0f..b4525b5 100644
--- a/geofileops/fileops.py
+++ b/geofileops/fileops.py
@@ -2273,7 +2273,7 @@ def _append_to_nolock(
     reproject: bool = False,
     explodecollections: bool = False,
     create_spatial_index: Optional[bool] = None,
-    force_output_geometrytype: Union[GeometryType, str, None] = None,
+    force_output_geometrytype: Union[GeometryType, str, Iterable[str], None] = None,
     transaction_size: int = 50000,
     preserve_fid: Optional[bool] = None,
     dst_dimensions: Optional[str] = None,
diff --git a/geofileops/geoops.py b/geofileops/geoops.py
index c5311d6..6081b80 100644
--- a/geofileops/geoops.py
+++ b/geofileops/geoops.py
@@ -107,9 +107,7 @@ def dissolve_within_distance(
         # First dissolve the input.
         #
         # Note: this reduces the complexity of operations to be executed later on.
-        # Note2: this already applies the gridsize, which needs to be applied anyway to
-        # avoid issues when determining the addedpieces_1neighbour later on.
-        # Note2: don' apply gridsize yet
+        # Note2: don't apply gridsize yet
         logger.info(f"Start, with input file {input_path}")
         step = 1
         logger.info(f"Step {step} of {nb_steps}")
@@ -342,6 +340,7 @@ def dissolve_within_distance(
             output_path=parts_to_add_filtered_path,
             sql_stmt=sql_stmt,
             input2_layer=input_layer,
+            explodecollections=True,
             gridsize=0.0,
             nb_parallel=nb_parallel,
             batchsize=batchsize,
diff --git a/geofileops/util/_geoops_ogr.py b/geofileops/util/_geoops_ogr.py
index ead1e62..98f2b1c 100644
--- a/geofileops/util/_geoops_ogr.py
+++ b/geofileops/util/_geoops_ogr.py
@@ -1,4 +1,5 @@
 import logging
+from collections.abc import Iterable
 from datetime import datetime
 from pathlib import Path
 from typing import Literal, Optional, Union
@@ -27,6 +28,15 @@ def clip_by_geometry(
         geom = wkt.loads(clip_geometry)
         spatial_filter = tuple(geom.bounds)
 
+    force_output_geometrytype = None
+    if not explodecollections:
+        input_layer_info = gfo.get_layerinfo(input_path, input_layer)
+        if input_layer_info.geometrytype is not GeometryType.POINT:
+            # If explodecollections is False and the input type is not point, force the
+            # output type to multi, because clip can cause eg. polygons to be split to
+            # multipolygons.
+            force_output_geometrytype = "PROMOTE_TO_MULTI"
+
     _run_ogr(
         operation="clip_by_geometry",
         input_path=input_path,
@@ -37,6 +47,7 @@ def clip_by_geometry(
         output_layer=output_layer,
         columns=columns,
         explodecollections=explodecollections,
+        force_output_geometrytype=force_output_geometrytype,
         force=force,
     )
 
@@ -112,7 +123,7 @@ def _run_ogr(
     append: bool = False,
     update: bool = False,
     explodecollections: bool = False,
-    force_output_geometrytype: Optional[GeometryType] = None,
+    force_output_geometrytype: Union[GeometryType, str, Iterable[str], None] = None,
     options: dict = {},
     columns: Optional[list[str]] = None,
     warp: Optional[dict] = None,
diff --git a/geofileops/util/_geoops_sql.py b/geofileops/util/_geoops_sql.py
index 9e7f66a..6f57394 100644
--- a/geofileops/util/_geoops_sql.py
+++ b/geofileops/util/_geoops_sql.py
@@ -1023,10 +1023,13 @@ def erase(
     input_layer_info = gfo.get_layerinfo(input_path, input_layer)
     primitivetypeid = input_layer_info.geometrytype.to_primitivetype.value
 
-    # If explodecollections is False and the input type is not point, force the output
-    # type to multi, because erase can cause eg. polygons to be split to multipolygons.
     force_output_geometrytype = input_layer_info.geometrytype
-    if not explodecollections and force_output_geometrytype is not GeometryType.POINT:
+    if explodecollections:
+        force_output_geometrytype = force_output_geometrytype.to_singletype
+    elif force_output_geometrytype is not GeometryType.POINT:
+        # If explodecollections is False and the input type is not point, force the
+        # output type to multi, because erase can cause eg. polygons to be split to
+        # multipolygons.
         force_output_geometrytype = force_output_geometrytype.to_multitype
 
     # Subdivide the erase layer if applicable to speed up further processing.
diff --git a/geofileops/util/_ogr_util.py b/geofileops/util/_ogr_util.py
index 32b83de..eb5ed62 100644
--- a/geofileops/util/_ogr_util.py
+++ b/geofileops/util/_ogr_util.py
@@ -12,7 +12,7 @@ from osgeo import gdal, ogr
 from pygeoops import GeometryType
 
 import geofileops as gfo
-from geofileops import fileops
+from geofileops import _compat, fileops
 
 # Make sure only one instance per process is running
 lock = Lock()
@@ -133,7 +133,7 @@ class VectorTranslateInfo:
         append: bool = False,
         update: bool = False,
         explodecollections: bool = False,
-        force_output_geometrytype: Union[GeometryType, str, None] = None,
+        force_output_geometrytype: Union[GeometryType, str, Iterable[str], None] = None,
         options: dict = {},
         columns: Optional[Iterable[str]] = None,
         warp: Optional[dict] = None,
@@ -208,7 +208,7 @@ def vector_translate(
     append: bool = False,
     update: bool = False,
     explodecollections: bool = False,
-    force_output_geometrytype: Union[GeometryType, str, None] = None,
+    force_output_geometrytype: Union[GeometryType, str, Iterable[str], None] = None,
     options: dict = {},
     columns: Optional[Iterable[str]] = None,
     warp: Optional[dict] = None,
@@ -319,11 +319,27 @@ def vector_translate(
     if force_output_geometrytype is not None:
         if isinstance(force_output_geometrytype, GeometryType):
             output_geometrytypes.append(force_output_geometrytype.name)
-        else:
+        elif isinstance(force_output_geometrytype, str):
             output_geometrytypes.append(force_output_geometrytype)
-    else:
-        if not explodecollections:
-            output_geometrytypes.append("PROMOTE_TO_MULTI")
+        elif isinstance(force_output_geometrytype, Iterable):
+            for geotype in force_output_geometrytype:
+                if isinstance(geotype, GeometryType):
+                    output_geometrytypes.append(geotype.name)
+                elif isinstance(geotype, str):
+                    output_geometrytypes.append(geotype)
+                else:
+                    raise ValueError(f"invalid type in {force_output_geometrytype=}")
+        else:
+            raise ValueError(f"invalid type for {force_output_geometrytype=}")
+    elif (
+        not explodecollections
+        and input_info.driver == "ESRI Shapefile"
+        and output_info.driver != "ESRI Shapefile"
+    ):
+        # Shapefiles are always reported as singlepart type but can also contain
+        # multiparts geometries, so promote to multi
+        output_geometrytypes.append("PROMOTE_TO_MULTI")
+
     if transaction_size is not None:
         args.extend(["-gt", str(transaction_size)])
     if preserve_fid is None:
@@ -380,6 +396,10 @@ def vector_translate(
     # Now we can really get to work
     output_ds = None
     try:
+        # Till gdal 3.10 datetime columns can be interpreted wrongly with arrow.
+        if _compat.GDAL_STE_310:
+            config_options["OGR2OGR_USE_ARROW_API"] = False
+
         # Go!
         with set_config_options(config_options):
             # Open input datasource already
diff --git a/pyproject.toml b/pyproject.toml
index f6282e1..4aae0f5 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -11,6 +11,9 @@ ignore_missing_imports = true
 
 [tool.pytest.ini_options]
 log_level = "DEBUG"
+# Write all logging and warnings immediately. Useful for debugging.
+# log_cli = true
+# addopts = "-p no:warnings"
 
 # Some filters that apply to warnings that can be ignored for any test.
 filterwarnings = [


Test patch:
diff --git a/tests/test_geofile.py b/tests/test_geofile.py
index db051eb..0ffa60b 100644
--- a/tests/test_geofile.py
+++ b/tests/test_geofile.py
@@ -316,9 +316,15 @@ def test_copy_layer_explodecollections(tmp_path, testfile, expected_count):
         ("polygon-parcel", GeometryType.MULTIPOINT),
     ],
 )
+@pytest.mark.filterwarnings(
+    "ignore:.*A geometry of type MULTIPOLYGON is inserted into .*"
+)
 def test_copy_layer_force_output_geometrytype(tmp_path, testfile, force_geometrytype):
     # The conversion is done by ogr, and the "test" is rather written to
-    # explore the behaviour of this ogr functionality
+    # explore the behaviour of this ogr functionality:
+    # Single-part polygons are converted to the destination types, but multipolygons
+    # are kept as they are.
+    # Issue opened for this: https://github.com/OSGeo/gdal/issues/11068
 
     # copy_layer on testfile and force to force_geometrytype
     src = test_helper.get_testfile(testfile)
@@ -954,7 +960,7 @@ def test_read_file_curve():
     # Test
     read_gdf = gfo.read_file(src)
     assert isinstance(read_gdf, gpd.GeoDataFrame)
-    assert isinstance(read_gdf.geometry[0], sh_geom.MultiPolygon)
+    assert isinstance(read_gdf.geometry[0], sh_geom.Polygon)
 
 
 def test_read_file_invalid_params(tmp_path, engine_setter):
diff --git a/tests/test_geofileops_twolayers.py b/tests/test_geofileops_twolayers.py
index 105fbc1..cee2052 100644
--- a/tests/test_geofileops_twolayers.py
+++ b/tests/test_geofileops_twolayers.py
@@ -413,6 +413,7 @@ def test_export_by_location(
         ("within is False", None, 1000, 0, 15),
     ],
 )
+@pytest.mark.filterwarnings("ignore:.*Field format '' not supported.*")
 def test_export_by_location_area(
     tmp_path,
     query,
@@ -1395,6 +1396,7 @@ def test_select_two_layers(tmp_path, suffix, epsg, gridsize):
 
 @pytest.mark.parametrize("suffix", SUFFIXES_GEOOPS)
 @pytest.mark.parametrize("input_nogeom", ["input1", "input2", "both"])
+@pytest.mark.filterwarnings("ignore:.*Field format '' not supported.*")
 def test_select_two_layers_input_without_geom(tmp_path, suffix, input_nogeom):
     # Prepare test file with geom
     input_geom_path = test_helper.get_testfile("polygon-parcel", suffix=suffix)
diff --git a/tests/test_ogr_util.py b/tests/test_ogr_util.py
index ae4da87..fde7b5c 100644
--- a/tests/test_ogr_util.py
+++ b/tests/test_ogr_util.py
@@ -6,6 +6,7 @@ import os
 
 import pytest
 from osgeo import gdal
+from pygeoops import GeometryType
 
 import geofileops as gfo
 from geofileops._compat import GDAL_GTE_38
@@ -148,6 +149,49 @@ def test_set_config_options():
     assert gdal.GetConfigOption(test3_config_envset) == "test3_new_env_value"
 
 
+@pytest.mark.parametrize(
+    "output_geometrytype, exp_geometrytype",
+    [
+        (None, "POLYGON"),
+        ("MULTIPOLYGON", "MULTIPOLYGON"),
+        (GeometryType.MULTIPOLYGON, "MULTIPOLYGON"),
+        (["MULTIPOLYGON"], "MULTIPOLYGON"),
+        ([GeometryType.MULTIPOLYGON], "MULTIPOLYGON"),
+        (["PROMOTE_TO_MULTI"], "MULTIPOLYGON"),
+    ],
+)
+def test_vector_translate_geometrytype(tmp_path, output_geometrytype, exp_geometrytype):
+    input_path = test_helper.get_testfile(
+        "polygon-parcel", dst_dir=tmp_path, explodecollections=True
+    )
+    output_path = tmp_path / f"output{input_path.suffix}"
+
+    _ogr_util.vector_translate(
+        str(input_path), output_path, force_output_geometrytype=output_geometrytype
+    )
+    output_info = gfo.get_layerinfo(output_path)
+    assert output_info.geometrytypename == exp_geometrytype
+
+
+@pytest.mark.parametrize(
+    "output_geometrytype, exp_error",
+    [
+        (1, "invalid type for"),
+        ([1], "invalid type in"),
+    ],
+)
+def test_vector_translate_geometrytype_error(tmp_path, output_geometrytype, exp_error):
+    input_path = test_helper.get_testfile(
+        "polygon-parcel", dst_dir=tmp_path, explodecollections=True
+    )
+    output_path = tmp_path / f"output{input_path.suffix}"
+
+    with pytest.raises(ValueError, match=exp_error):
+        _ogr_util.vector_translate(
+            str(input_path), output_path, force_output_geometrytype=output_geometrytype
+        )
+
+
 def test_vector_translate_input_nolayer(tmp_path):
     input_path = test_helper.get_testfile("polygon-parcel", dst_dir=tmp_path)
     output_path = tmp_path / f"output{input_path.suffix}"