Repository: geofileops/geofileops
Pull request #673 (2025-03-25 09:05:09)

ENH: add basic support for .gpkg.zip and .shp.zip files in general file/layer operations
Following issues need to be solved upstream:
- [x] For to_file there seems to be an issue in pyogrio, PR: https://github.com/geopandas/pyogrio/pull/527

ref #674 

codecov[bot]: ## [Codecov](https://app.codecov.io/gh/geofileops/geofileops/pull/673?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 `91.46341%` with `7 lines` in your changes missing coverage. Please review.
> Project coverage is 92.96%. 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 [(`5dfae34`)](https://app.codecov.io/gh/geofileops/geofileops/commit/5dfae348981bfd29440a3e4d6a081e93a2b4dfa2?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=geofileops).
> Report is 136 commits behind head on main.

| [Files with missing lines](https://app.codecov.io/gh/geofileops/geofileops/pull/673?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/fileops.py](https://app.codecov.io/gh/geofileops/geofileops/pull/673?src=pr&el=tree&filepath=geofileops%2Ffileops.py&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=geofileops#diff-Z2VvZmlsZW9wcy9maWxlb3BzLnB5) | 75.00% | [7 Missing :warning: ](https://app.codecov.io/gh/geofileops/geofileops/pull/673?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     #673      +/-   ##
==========================================
- Coverage   93.11%   92.96%   -0.16%     
==========================================
  Files          21       23       +2     
  Lines        4113     4664     +551     
==========================================
+ Hits         3830     4336     +506     
- Misses        283      328      +45     
```

</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/geofileops/geofileops/pull/673?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).

<details><summary> :rocket: New features to boost your workflow: </summary>

- :snowflake: [Test Analytics](https://docs.codecov.com/docs/test-analytics): Detect flaky tests, report on failures, and find test suite problems.
</details>

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 0b4828c..8215622 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -24,7 +24,9 @@
 - Improve performance of `export_by_location` and `join_by_location` for simple queries
   (#548)
 - Add support for `query=""` in `export_by_location` (#597)
-- Support GDAL vsi handlers in paths in most general file/layer operations (#669)
+- Add support for GDAL vsi handlers in paths in most general file/layer operations (#669)
+- Add basic support for ".gpkg.zip" and ".shp.zip" files in most general file/layer
+  operations (#673)
 - Add support to rename columns and layers with only a difference in casing (#549, #593)
 - Use `ST_Equals` and add priority feature to `delete_duplicate_geometries` (#638)
 - Avoid integer overflow when gpkg written by geofileops is read from .NET (#612)
diff --git a/geofileops/fileops.py b/geofileops/fileops.py
index 214e382..2897936 100644
--- a/geofileops/fileops.py
+++ b/geofileops/fileops.py
@@ -3,21 +3,18 @@
 import enum
 import filecmp
 import logging
+import os
 import pprint
 import shutil
 import string
 import tempfile
 import time
 import warnings
+import zipfile
 from collections.abc import Iterable
 from datetime import date, datetime
 from pathlib import Path
-from typing import (
-    TYPE_CHECKING,
-    Any,
-    Literal,
-    Union,
-)
+from typing import Any, Literal, Union
 
 import geopandas as gpd
 import numpy as np
@@ -32,15 +29,13 @@ from pygeoops import GeometryType, PrimitiveType  # noqa: F401
 from geofileops.helpers._configoptions_helper import ConfigOptions
 from geofileops.util import (
     _geofileinfo,
+    _geopath_util,
     _geoseries_util,
     _io_util,
     _ogr_sql_util,
     _ogr_util,
 )
 
-if TYPE_CHECKING:  # pragma: no cover
-    import os
-
 logger = logging.getLogger(__name__)
 
 # Enable exceptions for GDAL
@@ -111,11 +106,10 @@ def listlayers(
         <a href="https://gdal.org/en/stable/user/virtual_file_systems.html" target="_blank">GDAL vsi</a>
 
     """  # noqa: E501
-    path_p = Path(path)
-    if path_p.suffix.lower() == ".shp":
+    if str(path).lower().endswith((".shp", ".shp.zip")):
         if not _vsi_exists(path):
             raise FileNotFoundError(f"File not found: {path}")
-        return [path_p.stem]
+        return [_geopath_util.stem(path)]
 
     datasource = None
     try:
@@ -534,7 +528,7 @@ def get_default_layer(path: Union[str, "os.PathLike[Any]"]) -> str:
     Returns:
         str: The default layer name.
     """
-    return Path(path).stem
+    return _geopath_util.stem(path)
 
 
 def execute_sql(
@@ -609,6 +603,25 @@ def create_spatial_index(
         layer = get_layerinfo(path, layer, raise_on_nogeom=not no_geom_ok)
     if no_geom_ok and layer.geometrycolumn is None:
         return
+    if exist_ok and force_rebuild:
+        raise ValueError("exist_ok and force_rebuild can't both be True")
+
+    # .gpkg.zip files don't support update, so use has_spatial_index up-front to check
+    # if there is an index.
+    if str(path).lower().endswith(".gpkg.zip"):
+        if has_spatial_index(path, layer):
+            if force_rebuild:
+                pass
+            elif exist_ok:
+                return
+            else:
+                raise RuntimeError(
+                    f"spatial index already exists on {path}#{layer.name}"
+                )
+
+        raise RuntimeError(
+            f"create_spatial_index not supported for .gpkg.zip files: {path}"
+        )
 
     # Add index
     path_info = _geofileinfo.get_geofileinfo(path)
@@ -767,7 +780,9 @@ def remove_spatial_index(
 
     # Now really remove index
     try:
-        if path_info.is_spatialite_based:
+        if path_info.is_spatialite_based and not str(path).lower().endswith(
+            ".gpkg.zip"
+        ):
             datasource = gdal.OpenEx(str(path), nOpenFlags=gdal.OF_UPDATE)
             result = datasource.ExecuteSQL(
                 "SELECT DisableSpatialIndex("
@@ -2712,6 +2727,39 @@ def copy_layer(
     _ogr_util.vector_translate_by_info(info=translate_info)
 
 
+def _zip(src: Union[str, "os.PathLike[Any]"], dst: Union[str, "os.PathLike[Any]"]):
+    """Zip a file or directory.
+
+    Args:
+        src (PathLike): the file or directory to zip.
+        dst (PathLike): the destination zip file.
+    """
+    # Init input parameters
+    src = Path(src)
+
+    # Zip the file or directory
+    with zipfile.ZipFile(dst, "w") as zipf:
+        if src.is_file():
+            zipf.write(src, src.name)
+        elif src.is_dir():
+            for root, _, files in os.walk(src):
+                for file in files:
+                    file_path = Path(root) / file
+                    zipf.write(file_path, file_path.relative_to(src))
+
+
+def _unzip(src: Union[str, "os.PathLike[Any]"], dst: Union[str, "os.PathLike[Any]"]):
+    """Unzip a zip file.
+
+    Args:
+        src (PathLike): the zip file to unzip.
+        dst (PathLike): the destination directory.
+    """
+    # Unzip the file
+    with zipfile.ZipFile(src, "r") as zipf:
+        zipf.extractall(dst)
+
+
 def _determine_access_mode(
     dst: Union[str, "os.PathLike[Any]"],
     dst_layer: str | None,
diff --git a/geofileops/util/_geopath_util.py b/geofileops/util/_geopath_util.py
new file mode 100644
index 0000000..0876814
--- /dev/null
+++ b/geofileops/util/_geopath_util.py
@@ -0,0 +1,105 @@
+"""Helper functions for working with paths to geo files."""
+
+import os
+from pathlib import Path, PurePath
+from typing import Any, Union
+
+GEO_MULTI_SUFFIXES = [".gpkg.zip", ".shp.zip"]
+
+
+def stem(path: Union[str, "os.PathLike[Any]"]) -> str:
+    """Return the stem of a path, supporting multi-suffixes for geo formats.
+
+    Args:
+        path (PathLike): The path to get the stem of.
+
+    Returns:
+        str: the stem of the path.
+    """
+    # If there are multiple "."s in the name, check if it is a geo multi-suffix
+    path = PurePath(path)
+    if path.name.count(".") > 1:
+        name_lower = path.name.lower()
+        is_geo_multi_suffix = False
+        for suffixes in GEO_MULTI_SUFFIXES:
+            if name_lower.endswith(suffixes):
+                is_geo_multi_suffix = True
+                break
+
+        # A geo multi-suffix has been found
+        if is_geo_multi_suffix:
+            return path.name[: -len(suffixes)]
+
+    return path.stem
+
+
+def suffixes(path: Union[str, "os.PathLike[Any]"]) -> str:
+    """Return the suffixes of a path, supporting multi-suffixes for geo formats.
+
+    Args:
+        path (PathLike): The path to get the suffixes of.
+
+    Returns:
+        str: The suffixes of the path.
+    """
+    # If there are multiple "."s in the name, check if it is a geo multi-suffix
+    path = PurePath(path)
+    if path.name.count(".") > 1:
+        name_lower = path.name.lower()
+        for suffixes in GEO_MULTI_SUFFIXES:
+            if name_lower.endswith(suffixes):
+                return path.name[-len(suffixes) :]
+
+    return path.suffix
+
+
+def with_stem(path: Union[str, "os.PathLike[Any]"], new_stem: str) -> str | Path:
+    """Return a Path with a new stem, supporting multi-suffixes for geo formats.
+
+    If the input is a string, the output will be a string. If the input is a Path, the
+    output will be a Path.
+
+    Args:
+        path (PathLike): The path to change the stem of.
+        new_stem (str): The new stem to set.
+
+    Returns:
+        str | Path: The path with the new stem.
+
+    """
+    # If there are multiple "."s in the name, check if it is a geo multi-suffix
+    path_p = PurePath(path)
+
+    new_path_str = None
+    if path_p.name.count(".") > 1:
+        # Determine if the filename has a geo multi-suffix (case insensitive)
+        name_lower = path_p.name.lower()
+        is_geo_multi_suffix = False
+        for suffixes in GEO_MULTI_SUFFIXES:
+            if name_lower.endswith(suffixes):
+                is_geo_multi_suffix = True
+                break
+
+        # The file name has a geo multi-suffix
+        if is_geo_multi_suffix:
+            # Recover the suffixes from the file name to keep the casing
+            suffixes = path_p.name[-len(suffixes) :]
+            new_name = f"{new_stem}{suffixes}"
+            # Use os.path to retain forward/backward slashes in the path
+            parent, _ = os.path.split(path)
+            new_path_str = f"{parent}/{new_name}"
+
+    # If the input was a Path, return a Path
+    if isinstance(path, Path):
+        if new_path_str is not None:
+            return Path(new_path_str)
+        else:
+            return Path(path_p).with_stem(new_stem)
+
+    # The input was not a Path, return a string
+    if new_path_str is not None:
+        return new_path_str
+    else:
+        # Retain forward/backward slashes in the path by using os.path.
+        parent, tail = os.path.split(path)
+        return f"{parent}/{new_stem}{path_p.suffix}"
diff --git a/geofileops/util/_io_util.py b/geofileops/util/_io_util.py
index e0ae7e9..e8223d8 100644
--- a/geofileops/util/_io_util.py
+++ b/geofileops/util/_io_util.py
@@ -161,12 +161,6 @@ def create_file_atomic(path: Union[str, "os.PathLike[Any]"]) -> bool:
             raise RuntimeError(f"Error creating file {path}") from ex
 
 
-def with_stem(path: Path, new_stem) -> Path:
-    # Remark: from python 3.9 this is available on any Path, but to avoid
-    # having to require 3.9 for this, this hack...
-    return path.parent / f"{new_stem}{path.suffix}"
-
-
 def output_exists(path: Path, remove_if_exists: bool) -> bool:
     """Checks and returns whether the file specified exists.
 
diff --git a/geofileops/util/_ogr_util.py b/geofileops/util/_ogr_util.py
index 550e353..c8d1bc5 100644
--- a/geofileops/util/_ogr_util.py
+++ b/geofileops/util/_ogr_util.py
@@ -14,6 +14,7 @@ from pygeoops import GeometryType
 
 import geofileops as gfo
 from geofileops import _compat, fileops
+from geofileops.util import _geopath_util
 from geofileops.util._general_util import MissingRuntimeDependencyError
 
 # Make sure only one instance per process is running
@@ -395,7 +396,7 @@ def vector_translate(
 
     # Shapefiles only can have one layer, and the layer name == the stem of the file
     if output_info.driver == "ESRI Shapefile":
-        output_layer = Path(output_path).stem
+        output_layer = _geopath_util.stem(output_path)
 
     # SRS
     if output_srs is not None and isinstance(output_srs, int):
diff --git a/geofileops/util/_sqlite_util.py b/geofileops/util/_sqlite_util.py
index ac706be..c1aab1c 100644
--- a/geofileops/util/_sqlite_util.py
+++ b/geofileops/util/_sqlite_util.py
@@ -16,7 +16,7 @@ from pyproj import CRS, Transformer
 
 import geofileops as gfo
 from geofileops.helpers._configoptions_helper import ConfigOptions
-from geofileops.util import _sqlite_userdefined as sqlite_userdefined
+from geofileops.util import _geopath_util, _sqlite_userdefined
 from geofileops.util._general_util import MissingRuntimeDependencyError
 
 if TYPE_CHECKING:  # pragma: no cover
@@ -220,17 +220,22 @@ def get_columns(
     start = time.perf_counter()
     tmp_dir = None
 
+    def get_filetype(path: Path) -> str:
+        if _geopath_util.suffixes(path).lower() in (".gpkg", ".gpkg.zip"):
+            return "gpkg"
+        return path.suffix.lstrip(".")
+
     # Connect to/create sqlite main database
     if "main" in input_databases:
         # If an input database is main, use it as the main database
         main_db_path = input_databases["main"]
-        filetype = main_db_path.suffix.lstrip(".")
+        filetype = get_filetype(main_db_path)
         conn = sqlite3.connect(main_db_path, detect_types=sqlite3.PARSE_DECLTYPES)
         new_db = False
     else:
         # Create temp output db to be sure the output DB is writable, even though we
         # only create a temporary table.
-        filetype = next(iter(input_databases.values())).suffix.lstrip(".")
+        filetype = get_filetype(next(iter(input_databases.values())))
         conn = create_new_spatialdb(":memory:", filetype=filetype)
         new_db = True
 
@@ -831,7 +836,7 @@ def load_spatialite(conn):
     conn.create_function(
         "GFO_ReducePrecision",
         -1,
-        sqlite_userdefined.gfo_reduceprecision,
+        _sqlite_userdefined.gfo_reduceprecision,
         deterministic=True,
     )
 


Test patch:
diff --git a/tests/test_geofile.py b/tests/test_geofile.py
index 7b75662..179925b 100644
--- a/tests/test_geofile.py
+++ b/tests/test_geofile.py
@@ -5,6 +5,7 @@ Tests for functionalities in geofileops.general.
 import os
 import shutil
 from itertools import product
+from pathlib import Path
 
 import geopandas as gpd
 import pandas as pd
@@ -16,10 +17,14 @@ from pygeoops import GeometryType
 
 import geofileops as gfo
 from geofileops import fileops
-from geofileops.util import _geofileinfo, _geoseries_util
+from geofileops.util import _geofileinfo, _geopath_util, _geoseries_util
 from geofileops.util._geofileinfo import GeofileInfo
 from tests import test_helper
-from tests.test_helper import SUFFIXES_FILEOPS, assert_geodataframe_equal
+from tests.test_helper import (
+    SUFFIXES_FILEOPS,
+    SUFFIXES_FILEOPS_EXT,
+    assert_geodataframe_equal,
+)
 
 try:
     import fiona  # noqa: F401
@@ -172,7 +177,7 @@ def test_append_to(tmp_path):
     assert info.featurecount == 96
 
 
-@pytest.mark.parametrize("suffix", SUFFIXES_FILEOPS)
+@pytest.mark.parametrize("suffix", SUFFIXES_FILEOPS_EXT)
 def test_cmp(tmp_path, suffix):
     src = test_helper.get_testfile("polygon-parcel", suffix=suffix)
     src2 = test_helper.get_testfile("polygon-invalid", suffix=suffix)
@@ -208,22 +213,25 @@ def test_convert(tmp_path):
 
 
 @pytest.mark.parametrize(
-    "testfile, suffix_input, suffix_output, dimensions",
+    "testfile, suffix_input, suffix_output",
     [
-        *product(["polygon-parcel"], SUFFIXES_FILEOPS, SUFFIXES_FILEOPS, [None, "XYZ"]),
-        ["curvepolygon", ".gpkg", ".gpkg", None],
-        ["polygon-parcel", ".gpkg", ".gpkg.zip", None],
+        *product(["polygon-parcel"], SUFFIXES_FILEOPS_EXT, SUFFIXES_FILEOPS_EXT),
+        ["curvepolygon", ".gpkg", ".gpkg"],
     ],
 )
-def test_copy_layer(tmp_path, testfile, dimensions, suffix_input, suffix_output):
+def test_copy_layer(tmp_path, testfile, suffix_input, suffix_output):
+    if suffix_input == ".shp.zip" and suffix_output == ".shp":
+        # GDAL < 3.10 determines the layer name wrong for .shp.zip leading to this error
+        pytest.xfail("Copy of .shp.zip gives issues in GDAL <= 3.10")
+
     # Prepare test data
-    src = test_helper.get_testfile(testfile, suffix=suffix_input, dimensions=dimensions)
+    src = test_helper.get_testfile(testfile, suffix=suffix_input)
     if suffix_input == ".csv" or suffix_output == ".csv":
         raise_on_nogeom = False
     else:
         raise_on_nogeom = True
 
-    if suffix_input == ".csv" and suffix_output == ".shp":
+    if suffix_input == ".csv" and suffix_output in (".shp", ".shp.zip"):
         # If no geometry column, there will only be a .dbf output file
         dst = tmp_path / f"{src.stem}-output.dbf"
     else:
@@ -235,7 +243,7 @@ def test_copy_layer(tmp_path, testfile, dimensions, suffix_input, suffix_output)
     # Now compare source and dst file
     src_layerinfo = gfo.get_layerinfo(src, raise_on_nogeom=raise_on_nogeom)
     dst_layerinfo = gfo.get_layerinfo(dst, raise_on_nogeom=raise_on_nogeom)
-    assert dst_layerinfo.name == dst.stem
+    assert dst_layerinfo.name == gfo.get_default_layer(dst)
     assert src_layerinfo.featurecount == dst_layerinfo.featurecount
     assert len(src_layerinfo.columns) == len(dst_layerinfo.columns)
     if not (
@@ -245,6 +253,31 @@ def test_copy_layer(tmp_path, testfile, dimensions, suffix_input, suffix_output)
         assert src_layerinfo.geometrytypename == dst_layerinfo.geometrytypename
 
 
+@pytest.mark.parametrize(
+    "testfile, suffix_input, suffix_output",
+    [
+        *product(["polygon-parcel"], [".gpkg", ".shp"], [".gpkg", ".shp"]),
+        ["curvepolygon", ".gpkg", ".gpkg"],
+    ],
+)
+def test_copy_layer_dimensions_xyz(tmp_path, testfile, suffix_input, suffix_output):
+    # Prepare test data
+    src = test_helper.get_testfile(testfile, suffix=suffix_input, dimensions="XYZ")
+    dst = tmp_path / f"{src.stem}-output{suffix_output}"
+
+    # Test
+    gfo.copy_layer(str(src), str(dst))
+
+    # Now compare source and dst file
+    src_layerinfo = gfo.get_layerinfo(src)
+    dst_layerinfo = gfo.get_layerinfo(dst)
+    assert dst_layerinfo.name == dst.stem
+    assert src_layerinfo.featurecount == dst_layerinfo.featurecount
+    assert len(src_layerinfo.columns) == len(dst_layerinfo.columns)
+    if not (suffix_input == ".shp" and suffix_output == ".gpkg"):
+        assert src_layerinfo.geometrytypename == dst_layerinfo.geometrytypename
+
+
 def test_copy_layer_add_layer_gpkg(tmp_path):
     # Prepare test data
     src = test_helper.get_testfile("polygon-parcel")
@@ -777,6 +810,24 @@ def test_copy_layer_sql_placeholders(tmp_path, suffix):
     assert_geodataframe_equal(input_gdf, copy_gdf)
 
 
+def test_copy_layer_to_gpkg_zip(tmp_path):
+    # Prepare test data
+    src = test_helper.get_testfile("polygon-parcel")
+    dst = tmp_path / "output.gpkg.zip"
+
+    # copy_layer
+    gfo.copy_layer(src, dst)
+
+    # Now compare source and dst file
+    src_layerinfo = gfo.get_layerinfo(src)
+    dst_layerinfo = gfo.get_layerinfo(dst)
+    assert src_layerinfo.featurecount == dst_layerinfo.featurecount
+
+    src_gdf = gfo.read_file(src)
+    dst_gdf = gfo.read_file(dst)
+    assert_geodataframe_equal(src_gdf, dst_gdf)
+
+
 def test_copy_layer_vsi(tmp_path):
     # Prepare test data
     src = f"/vsizip//vsicurl/{test_helper.data_url}/poly_shp.zip/poly.shp"
@@ -864,7 +915,7 @@ def test_copy_layer_write_mode_force(tmp_path, write_mode, force):
         assert dst.stat().st_mtime == mtime_orig
 
 
-@pytest.mark.parametrize("suffix", SUFFIXES_FILEOPS)
+@pytest.mark.parametrize("suffix", SUFFIXES_FILEOPS_EXT)
 def test_copy(tmp_path, suffix):
     src = test_helper.get_testfile("polygon-parcel", suffix=suffix)
 
@@ -952,21 +1003,37 @@ def test_get_crs_vsi():
     assert crs.to_epsg() == 27700
 
 
+@pytest.mark.parametrize(
+    "path, exp_default_layer",
+    [
+        ("/tmp/polygons.gpkg", "polygons"),
+        (Path("/tmp/polygons.gpkg"), "polygons"),
+        ("/tmp/polygons.gpkg.zip", "polygons"),
+        (Path("/tmp/polygons.gpkg.zip"), "polygons"),
+        (r"C:\tmp\polygons.gpkg.zip", "polygons"),
+        (Path(r"C:\tmp\polygons.gpkg.zip"), "polygons"),
+        ("/tmp/polygons.shp", "polygons"),
+        ("/tmp/polygons.csv", "polygons"),
+        ("/tmp/polygons.csv", "polygons"),
+        ("/vsizip//vsicurl/poly_shp.zip/poly.shp", "poly"),
+    ],
+)
+def test_get_default_layer(path, exp_default_layer):
+    if os.name != "nt" and str(path).lower().startswith("c:"):
+        pytest.skip("Test only valid on Windows")
+
+    layer = gfo.get_default_layer(path)
+    assert layer == exp_default_layer
+
+
 @pytest.mark.parametrize("suffix", SUFFIXES_FILEOPS)
-def test_get_default_layer(suffix):
+def test_get_default_layer_files(suffix):
     # Prepare test data + test
     src = test_helper.get_testfile("polygon-parcel", suffix=suffix)
     layer = gfo.get_default_layer(str(src))
     assert layer == src.stem
 
 
-def test_get_default_layer_vsi():
-    # Prepare test data + test
-    src = f"/vsizip//vsicurl/{test_helper.data_url}/poly_shp.zip/poly.shp"
-    layer = gfo.get_default_layer(src)
-    assert layer == "poly"
-
-
 @pytest.mark.parametrize("suffix", [s for s in SUFFIXES_FILEOPS if s != ".csv"])
 def test_get_layer_geometrytypes(suffix):
     # Prepare test data + test
@@ -1006,9 +1073,14 @@ def test_get_layer_geometrytypes_vsi(tmp_path):
     assert geometrytypes == ["POLYGON"]
 
 
-@pytest.mark.parametrize("suffix", [".gpkg", ".shp"])
+@pytest.mark.parametrize("suffix", [".gpkg", ".gpkg.zip", ".shp", ".shp.zip"])
 @pytest.mark.parametrize("dimensions", [None, "XYZ"])
 def test_get_layerinfo(suffix, dimensions):
+    if dimensions == "XYZ" and suffix == ".gpkg.zip":
+        pytest.skip(
+            "get_testfile for dim=XYZ requires updating: so skip .gpkg.zip + XYZ"
+        )
+
     src = test_helper.get_testfile(
         "polygon-parcel", suffix=suffix, dimensions=dimensions
     )
@@ -1189,7 +1261,9 @@ def test_listlayers_errors(suffix):
     "suffix, only_spatial_layers, expected",
     [
         (".gpkg", True, ["parcels"]),
+        (".gpkg.zip", True, ["parcels"]),
         (".shp", True, ["{src_stem}"]),
+        (".shp.zip", True, ["{src_stem}"]),
         (".csv", True, []),
         (".csv", False, ["{src_stem}"]),
     ],
@@ -1199,7 +1273,7 @@ def test_listlayers_one_layer(suffix, only_spatial_layers, expected):
     src = test_helper.get_testfile("polygon-parcel", suffix=suffix)
     layers = gfo.listlayers(src, only_spatial_layers=only_spatial_layers)
 
-    expected = [exp.format(src_stem=src.stem) for exp in expected]
+    expected = [exp.format(src_stem=_geopath_util.stem(src)) for exp in expected]
     assert layers == expected
 
 
@@ -1218,7 +1292,7 @@ def test_listlayers_vsi():
     assert "poly" in layers
 
 
-@pytest.mark.parametrize("suffix", SUFFIXES_FILEOPS)
+@pytest.mark.parametrize("suffix", SUFFIXES_FILEOPS_EXT)
 def test_move(tmp_path, suffix):
     src = test_helper.get_testfile("polygon-parcel", dst_dir=tmp_path, suffix=suffix)
 
@@ -1295,11 +1369,16 @@ def test_update_column_error(tmp_path):
         gfo.update_column(test_path, name="OPPERVL", expression="invalid_expression")
 
 
-@pytest.mark.parametrize("suffix", SUFFIXES_FILEOPS)
+@pytest.mark.parametrize("suffix", SUFFIXES_FILEOPS_EXT)
 @pytest.mark.parametrize("dimensions", [None, "XYZ"])
 def test_read_file(suffix, dimensions, engine_setter):
     # Remark: it seems like Z dimensions aren't read in geopandas.
     # Prepare and validate test data
+    if dimensions == "XYZ" and suffix == ".gpkg.zip":
+        pytest.skip(
+            "get_testfile for dim=XYZ requires updating: so skip .gpkg.zip + XYZ"
+        )
+
     src = test_helper.get_testfile(
         "polygon-parcel", suffix=suffix, dimensions=dimensions
     )
@@ -1720,9 +1799,7 @@ def test_spatial_index(tmp_path, suffix):
     assert has_spatial_index is True
 
     # Spatial index if it already exists by default gives error
-    with pytest.raises(
-        Exception, match="create_spatial_index error: spatial index already exists"
-    ):
+    with pytest.raises(Exception, match="spatial index already exists"):
         gfo.create_spatial_index(path=test_path, layer=layer)
     gfo.create_spatial_index(path=test_path, layer=layer, exist_ok=True)
 
@@ -1736,6 +1813,55 @@ def test_spatial_index(tmp_path, suffix):
         assert qix_path.stat().st_mtime > qix_modified_time_orig
 
 
+def test_spatial_index_gpkg_zip(tmp_path):
+    """Spatial index tests are specific for .gpkg.zip as it is read-only."""
+    # Test cases where the input file has an index
+    # --------------------------------------------
+    test_path = test_helper.get_testfile(
+        "polygon-parcel", dst_dir=tmp_path, suffix=".gpkg.zip"
+    )
+    layer = gfo.get_only_layer(test_path)
+
+    # Check if spatial index present
+    has_spatial_index = gfo.has_spatial_index(path=test_path, layer=layer)
+    assert has_spatial_index is True
+
+    # Removing spatial index in not supported as .gpkg.zip is read-only
+    with pytest.raises(ValueError, match="remove_spatial_index not supported for"):
+        gfo.remove_spatial_index(path=test_path, layer=layer)
+
+    # Create spatial index is supported with `exist_ok=True`
+    gfo.create_spatial_index(path=test_path, layer=layer, exist_ok=True)
+
+    # If no `exist_ok=True` is specified, error
+    with pytest.raises(RuntimeError, match="spatial index already exists"):
+        gfo.create_spatial_index(path=test_path, layer=layer)
+
+    # If  `force_rebuild=True` is specified, error
+    with pytest.raises(
+        RuntimeError, match="create_spatial_index not supported for .gpkg.zip files"
+    ):
+        gfo.create_spatial_index(path=test_path, layer=layer, force_rebuild=True)
+
+    # Test cases where the input file does not have an index
+    # ------------------------------------------------------
+    # Prepare .gpkg file without spatial index
+    test_path = test_helper.get_testfile(
+        "polygon-parcel", dst_dir=tmp_path, suffix=".gpkg"
+    )
+    layer = gfo.get_only_layer(test_path)
+    gfo.remove_spatial_index(path=test_path, layer=layer)
+    assert not gfo.has_spatial_index(path=test_path, layer=layer)
+    test_zip_path = tmp_path / f"{test_path.name}.zip"
+    fileops._zip(test_path, test_zip_path)
+    assert not gfo.has_spatial_index(path=test_path, layer=layer)
+
+    with pytest.raises(
+        RuntimeError, match="create_spatial_index not supported for .gpkg.zip files"
+    ):
+        gfo.create_spatial_index(path=test_zip_path, layer=layer)
+
+
 def test_spatial_index_unsupported(tmp_path):
     # Prepare test data
     suffix = ".geojson"
@@ -1769,7 +1895,17 @@ def test_spatial_index_unsupported(tmp_path):
         _ = gfo.remove_spatial_index(path, path.stem)
 
 
-@pytest.mark.parametrize("suffix", SUFFIXES_FILEOPS)
+def test_create_spatial_index_invalid_params():
+    path = test_helper.get_testfile("polygon-parcel")
+
+    # Test invalid parameters
+    with pytest.raises(
+        ValueError, match="exist_ok and force_rebuild can't both be True"
+    ):
+        _ = gfo.create_spatial_index(path=path, exist_ok=True, force_rebuild=True)
+
+
+@pytest.mark.parametrize("suffix", SUFFIXES_FILEOPS_EXT)
 @pytest.mark.parametrize("dimensions", [None])
 def test_to_file(tmp_path, suffix, dimensions, engine_setter):
     # Remark: geopandas doesn't seem seem to read the Z dimension, so writing can't be
@@ -1778,7 +1914,7 @@ def test_to_file(tmp_path, suffix, dimensions, engine_setter):
     src = test_helper.get_testfile(
         "polygon-parcel", suffix=suffix, dimensions=dimensions
     )
-    output_path = tmp_path / f"{src.stem}-output{suffix}"
+    output_path = tmp_path / f"{_geopath_util.stem(src)}-output{suffix}"
     uidn = str(2318781) if suffix == ".csv" else 2318781
     encoding = "utf-8" if suffix == ".csv" else None
 
@@ -1788,6 +1924,9 @@ def test_to_file(tmp_path, suffix, dimensions, engine_setter):
     # Validate if string (encoding) is correct for data read.
     assert read_gdf.loc[read_gdf["UIDN"] == uidn]["LBLHFDTLT"].item() == "Silomaïs"
 
+    if suffix in (".gpkg.zip", ".shp.zip"):
+        pytest.xfail("writing a dataframe to gpkg.zip or .shp.zip has issue")
+
     gfo.to_file(read_gdf, str(output_path))
     written_gdf = gfo.read_file(output_path)
 
@@ -1811,10 +1950,21 @@ def test_to_file(tmp_path, suffix, dimensions, engine_setter):
 
         assert_geodataframe_equal(written_gdf, read_gdf)
 
-    # Append the file again to tmppath
-    gfo.to_file(read_gdf, output_path, append=True)
-    written_gdf = gfo.read_file(output_path)
-    assert 2 * len(read_gdf) == len(written_gdf)
+
+@pytest.mark.parametrize("suffix", SUFFIXES_FILEOPS)
+def test_to_file_append(tmp_path, suffix, engine_setter):
+    test_path = test_helper.get_testfile(
+        "polygon-parcel", dst_dir=tmp_path, suffix=suffix
+    )
+    raise_on_nogeom = False if suffix == ".csv" else True
+
+    test_gdf = gfo.read_file(test_path)
+    gfo.to_file(test_gdf, path=test_path, append=True)
+
+    # Check result
+    assert test_path.exists()
+    dst_info = gfo.get_layerinfo(test_path, raise_on_nogeom=raise_on_nogeom)
+    assert dst_info.featurecount == len(test_gdf) * 2
 
 
 @pytest.mark.parametrize("suffix", SUFFIXES_FILEOPS)
@@ -2327,3 +2477,37 @@ def test_launder_columns():
         NotImplementedError, match="Not supported to launder > 99 columns starting with"
     ):
         laundered = fileops._launder_column_names(columns)
+
+
+def test_zip_unzip(tmp_path):
+    # Prepare test data
+    src = test_helper.get_testfile("polygon-parcel")
+    zip_path = tmp_path / "zipped.zip"
+    fileops._zip(src, zip_path)
+
+    # Unzip and check result
+    dst_dir = tmp_path / "unzipped"
+    fileops._unzip(zip_path, dst_dir)
+    assert len(list(dst_dir.iterdir())) == 1
+    assert (dst_dir / src.name).exists()
+
+
+def test_zip_unzip_dir(tmp_path):
+    # Prepare test data
+    src = test_helper.get_testfile("polygon-parcel")
+    zip_dir = tmp_path / "dir_to_zip"
+    zip_dir.mkdir()
+    file1 = zip_dir / f"{src.stem}_1{src.suffix}"
+    file2 = zip_dir / f"{src.stem}_2{src.suffix}"
+    gfo.copy(src, file1)
+    gfo.copy(src, file2)
+    zip_path = tmp_path / "zipped.zip"
+    fileops._zip(zip_dir, zip_path)
+
+    # Unzip and check result
+    dst_dir = tmp_path / "unzipped"
+    fileops._unzip(zip_path, dst_dir)
+    assert dst_dir.exists()
+    assert len(list(dst_dir.iterdir())) == 2
+    assert (dst_dir / file1.name).exists()
+    assert (dst_dir / file2.name).exists()
diff --git a/tests/test_geofileops_singlelayer.py b/tests/test_geofileops_singlelayer.py
index 5c8f42b..2d0a172 100644
--- a/tests/test_geofileops_singlelayer.py
+++ b/tests/test_geofileops_singlelayer.py
@@ -14,14 +14,14 @@ from shapely import MultiPolygon, Polygon
 
 from geofileops import GeometryType, fileops, geoops
 from geofileops._compat import SPATIALITE_GTE_51
-from geofileops.util import _general_util, _geofileinfo, _geoops_sql
-from geofileops.util import _io_util as io_util
+from geofileops.util import _general_util, _geofileinfo, _geoops_sql, _geopath_util
 from geofileops.util._geofileinfo import GeofileInfo
 from tests import test_helper
 from tests.test_helper import (
     EPSGS,
     GRIDSIZE_DEFAULT,
     SUFFIXES_GEOOPS,
+    SUFFIXES_GEOOPS_INPUT,
     TESTFILES,
     WHERE_AREA_GT_400,
     assert_geodataframe_equal,
@@ -149,12 +149,40 @@ def basic_combinations_to_test(
     return result
 
 
+@pytest.mark.parametrize("suffix_input", SUFFIXES_GEOOPS_INPUT)
+@pytest.mark.parametrize("geoops_module", GEOOPS_MODULES)
+def test_buffer(tmp_path, suffix_input, geoops_module):
+    """Buffer minimal test."""
+    # Prepare test data
+    input_path = test_helper.get_testfile("polygon-parcel", suffix=suffix_input)
+    input_layerinfo = fileops.get_layerinfo(input_path)
+    batchsize = math.ceil(input_layerinfo.featurecount / 2)
+
+    # Now run test
+    output_path = tmp_path / "output.gpkg"
+    set_geoops_module(geoops_module)
+    geoops.buffer(
+        input_path=input_path,
+        output_path=output_path,
+        distance=1,
+        nb_parallel=2,
+        keep_empty_geoms=True,
+        batchsize=batchsize,
+    )
+
+    # Now check if the output file is correctly created
+    assert output_path.exists()
+    output_layerinfo = fileops.get_layerinfo(output_path)
+    assert len(output_layerinfo.columns) == len(input_layerinfo.columns)
+    assert output_layerinfo.featurecount == input_layerinfo.featurecount
+
+
 @pytest.mark.parametrize(
     "suffix, epsg, geoops_module, testfile, empty_input, gridsize, keep_empty_geoms, "
     "where_post, dimensions",
     basic_combinations_to_test(),
 )
-def test_buffer(
+def test_buffer_basic(
     tmp_path,
     suffix,
     epsg,
@@ -1114,7 +1142,9 @@ def test_simplify(
     kwargs = {}
     if keep_empty_geoms is not None:
         kwargs["keep_empty_geoms"] = keep_empty_geoms
-    output_path = io_util.with_stem(input_path, output_path)
+    output_path = _geopath_util.with_stem(
+        input_path, f"{output_path.stem}_{keep_empty_geoms}"
+    )
     geoops.simplify(
         input_path=input_path,
         output_path=output_path,
diff --git a/tests/test_geopath.py b/tests/test_geopath.py
new file mode 100644
index 0000000..9ce2188
--- /dev/null
+++ b/tests/test_geopath.py
@@ -0,0 +1,65 @@
+"""Test for functionalities in _geopath_util."""
+
+from pathlib import Path
+
+import pytest
+
+from geofileops.util import _geopath_util
+
+
+@pytest.mark.parametrize(
+    "path, exp_stem",
+    [
+        (Path("/tmp/testje.gpkg"), "testje"),
+        ("/tmp/testje.gpkg.zip", "testje"),
+        ("/tmp/testje.shp", "testje"),
+        ("/tmp/testje.shp.zip", "testje"),
+        ("/tmp/testje.txt", "testje"),
+        ("/tmp/testje", "testje"),
+        ("/tmp/testje.tar.gz", "testje.tar"),
+        ("/tmp/t.estj.e.gpkg", "t.estj.e"),
+        ("/tmp/t.estj.e.gpkg.zip", "t.estj.e"),
+    ],
+)
+def test_stem(path, exp_stem):
+    assert _geopath_util.stem(path) == exp_stem
+
+
+@pytest.mark.parametrize(
+    "path, exp_suffixes",
+    [
+        (Path("/tmp/testje.gpkg"), ".gpkg"),
+        ("/tmp/testje.gpkg.zip", ".gpkg.zip"),
+        ("/tmp/testje.shp", ".shp"),
+        ("/tmp/testje.shp.zip", ".shp.zip"),
+        ("/tmp/testje.txt", ".txt"),
+        ("/tmp/testje", ""),
+        ("/tmp/testje.tar.gz", ".gz"),
+        ("/tmp/t.estj.e.gpkg", ".gpkg"),
+        ("/tmp/t.estj.e.gpkg.zip", ".gpkg.zip"),
+    ],
+)
+def test_suffixes(path, exp_suffixes):
+    assert _geopath_util.suffixes(path) == exp_suffixes
+
+
+@pytest.mark.parametrize(
+    "path, new_stem, exp_path",
+    [
+        ("/tmp/testje.gpkg", "testje_2", "/tmp/testje_2.gpkg"),
+        ("/tmp/testje.gpkg.zip", "testje_2", "/tmp/testje_2.gpkg.zip"),
+        ("/tmp/testje.shp", "testje_2", "/tmp/testje_2.shp"),
+        ("/tmp/testje.shp.zip", "testje_2", "/tmp/testje_2.shp.zip"),
+        ("/tmp/testje.txt", "testje_2", "/tmp/testje_2.txt"),
+        ("/tmp/testje", "testje_2", "/tmp/testje_2"),
+        ("/tmp/testje.tar.gz", "testje_2", "/tmp/testje_2.gz"),
+        ("/tmp/t.estj.e.gpkg", "t.estj.e_2", "/tmp/t.estj.e_2.gpkg"),
+        ("/tmp/t.estj.e.gpkg.zip", "t.estj.e_2", "/tmp/t.estj.e_2.gpkg.zip"),
+    ],
+)
+def test_with_stem(path, new_stem, exp_path):
+    # If input is a string, output should be a string
+    assert _geopath_util.with_stem(path, new_stem) == exp_path
+
+    # If input is a Path, output should be a Path
+    assert _geopath_util.with_stem(Path(path), new_stem) == Path(exp_path)
diff --git a/tests/test_helper.py b/tests/test_helper.py
index 9eb93e1..3b97d72 100644
--- a/tests/test_helper.py
+++ b/tests/test_helper.py
@@ -12,7 +12,13 @@ import shapely
 import shapely.geometry as sh_geom
 
 import geofileops as gfo
-from geofileops.util import _geofileinfo, _geoseries_util, _io_util, geodataframe_util
+from geofileops.util import (
+    _geofileinfo,
+    _geopath_util,
+    _geoseries_util,
+    _io_util,
+    geodataframe_util,
+)
 
 data_dir = Path(__file__).parent.resolve() / "data"
 data_url = "https://raw.githubusercontent.com/geofileops/geofileops/main/tests/data"
@@ -20,7 +26,9 @@ data_url = "https://raw.githubusercontent.com/geofileops/geofileops/main/tests/d
 EPSGS = [31370, 4326]
 GRIDSIZE_DEFAULT = 0.0
 SUFFIXES_FILEOPS = [".gpkg", ".shp", ".csv"]
+SUFFIXES_FILEOPS_EXT = [".gpkg", ".gpkg.zip", ".shp.zip", ".shp", ".csv"]
 SUFFIXES_GEOOPS = [".gpkg", ".shp"]
+SUFFIXES_GEOOPS_INPUT = [".gpkg", ".gpkg.zip", ".shp.zip", ".shp"]
 TESTFILES = ["polygon-parcel", "linestring-row-trees", "point"]
 WHERE_AREA_GT_400 = "ST_Area({geometrycolumn}) > 400"
 WHERE_AREA_GT_5000 = "ST_Area({geometrycolumn}) > 5000"
@@ -173,7 +181,8 @@ def get_testfile(
 
         # Prepare the file in a tmp file so the file is not visible to other
         # processes until it is completely ready.
-        tmp_path = prepared_path.with_stem(f"{prepared_path.stem}_tmp")
+        tmp_stem = f"{_geopath_util.stem(prepared_path)}_tmp"
+        tmp_path = _geopath_util.with_stem(prepared_path, tmp_stem)
         layers = gfo.listlayers(testfile_path)
         dst_info = _geofileinfo.get_geofileinfo(tmp_path)
         if len(layers) > 1 and dst_info.is_singlelayer:
@@ -185,7 +194,13 @@ def get_testfile(
         # Convert all layers found
         for src_layer in layers:
             # Single layer files have stem as layername
-            dst_layer = tmp_path.stem if dst_info.is_singlelayer else src_layer
+            assert isinstance(tmp_path, Path)
+            if dst_info.is_singlelayer:
+                dst_layer = tmp_stem
+                preserve_fid = False
+            else:
+                dst_layer = src_layer
+                preserve_fid = not explodecollections
 
             gfo.copy_layer(
                 testfile_path,
@@ -195,7 +210,7 @@ def get_testfile(
                 write_mode="add_layer",
                 dst_crs=epsg,
                 reproject=True,
-                preserve_fid=not explodecollections,
+                preserve_fid=preserve_fid,
                 dst_dimensions=dimensions,
                 explodecollections=explodecollections,
             )
diff --git a/tests/test_sqlite_userdefined.py b/tests/test_sqlite_userdefined.py
index ac4d937..e74f0bb 100644
--- a/tests/test_sqlite_userdefined.py
+++ b/tests/test_sqlite_userdefined.py
@@ -4,7 +4,7 @@ import shapely.ops
 from shapely import MultiPolygon, Point, Polygon
 from shapely.testing import assert_geometries_equal
 
-from geofileops.util import _sqlite_userdefined as sqlite_userdefined
+from geofileops.util import _sqlite_userdefined
 
 
 def test_gfo_difference_collection():
@@ -12,12 +12,12 @@ def test_gfo_difference_collection():
     box0_5 = shapely.box(0, 0, 5, 5)
     box0_10 = shapely.box(0, 0, 10, 10)
     assert (
-        sqlite_userdefined.gfo_difference_collection(box0_10.wkb, box0_5.wkb)
+        _sqlite_userdefined.gfo_difference_collection(box0_10.wkb, box0_5.wkb)
         == shapely.difference(box0_10, box0_5).wkb
     )
     # Difference of polygons gives a linestring to test -> No testcase at hand
     assert (
-        sqlite_userdefined.gfo_difference_collection(
+        _sqlite_userdefined.gfo_difference_collection(
             box0_10.wkb, box0_5.wkb, keep_geom_type=1
         )
         == shapely.difference(box0_10, box0_5).wkb
@@ -25,7 +25,7 @@ def test_gfo_difference_collection():
 
     # Result of the difference is None.
     assert (
-        sqlite_userdefined.gfo_difference_collection(
+        _sqlite_userdefined.gfo_difference_collection(
             shapely.Point(1, 1).wkb, box0_10.wkb
         )
         is None
@@ -47,7 +47,7 @@ def test_gfo_difference_collection():
 def test_gfo_difference_collection_empty_geoms(
     test_id, geom, geoms_to_subtract, exp_result
 ):
-    result = sqlite_userdefined.gfo_difference_collection(geom, geoms_to_subtract)
+    result = _sqlite_userdefined.gfo_difference_collection(geom, geoms_to_subtract)
     if exp_result is None:
         assert result is None, f"Issue with test {test_id}"
     else:
@@ -57,22 +57,22 @@ def test_gfo_difference_collection_empty_geoms(
 def test_gfo_difference_collection_invalid_params():
     # geom_to_subtract is not a wkb
     with pytest.raises(TypeError, match="Expected bytes or string, got Point"):
-        sqlite_userdefined.gfo_difference_collection(
+        _sqlite_userdefined.gfo_difference_collection(
             shapely.Point(1, 1).wkb, shapely.Point(1, 1)
         )
     # geom is not a wkb
     with pytest.raises(TypeError, match="Expected bytes or string, got Point"):
-        sqlite_userdefined.gfo_difference_collection(
+        _sqlite_userdefined.gfo_difference_collection(
             shapely.Point(1, 1), shapely.Point(1, 1).wkb
         )
 
     # keep_geom_type should be int (0 or 1)
     with pytest.raises(TypeError, match="keep_geom_type must be int"):
-        sqlite_userdefined.gfo_difference_collection(
+        _sqlite_userdefined.gfo_difference_collection(
             shapely.Point(1, 1).wkb, shapely.Point(1, 1).wkb, keep_geom_type="TRUE"
         )
     with pytest.raises(ValueError, match="keep_geom_type has invalid value"):
-        sqlite_userdefined.gfo_difference_collection(
+        _sqlite_userdefined.gfo_difference_collection(
             shapely.Point(1, 1).wkb, shapely.Point(1, 1).wkb, keep_geom_type=5
         )
 
@@ -101,7 +101,7 @@ def test_gfo_reduceprecision(test_descr, geom, exp_result):
     geom_wkb = None if geom is None else geom.wkb
 
     # Run test
-    result_wkb = sqlite_userdefined.gfo_reduceprecision(geom_wkb, gridsize=1)
+    result_wkb = _sqlite_userdefined.gfo_reduceprecision(geom_wkb, gridsize=1)
 
     # Check result
     result = shapely.from_wkb(result_wkb)
@@ -125,7 +125,7 @@ def test_gfo_split(test):
 
     if test == "poly":
         # Split of polygon with linestring blade gives collection of polygons
-        result = shapely.from_wkb(sqlite_userdefined.gfo_split(box1_4.wkb, blade.wkb))
+        result = shapely.from_wkb(_sqlite_userdefined.gfo_split(box1_4.wkb, blade.wkb))
         assert result == shapely.ops.split(box1_4, blade)
 
     elif test == "multipoly":
@@ -133,7 +133,7 @@ def test_gfo_split(test):
         box5_9 = shapely.box(5, 0, 9, 4)
         multipoly = MultiPolygon([box1_4, box5_9])
         result = shapely.from_wkb(
-            sqlite_userdefined.gfo_split(multipoly.wkb, blade.wkb)
+            _sqlite_userdefined.gfo_split(multipoly.wkb, blade.wkb)
         )
         expected_result = """
             GEOMETRYCOLLECTION (POLYGON ((4 2, 4 0, 1 0, 1 2, 4 2)),
@@ -146,7 +146,7 @@ def test_gfo_split(test):
 
     elif test == "None":
         # Result of splitting None is None.
-        assert sqlite_userdefined.gfo_split(None, blade) is None
+        assert _sqlite_userdefined.gfo_split(None, blade) is None
 
     else:
         raise ValueError(f"test not implemented: {test}")
@@ -165,7 +165,7 @@ def test_gfo_split(test):
     ],
 )
 def test_gfo_split_empty_geoms(test_id, geom, blade, exp_result):
-    result = sqlite_userdefined.gfo_split(geom, blade)
+    result = _sqlite_userdefined.gfo_split(geom, blade)
     if exp_result is None:
         assert result is None, f"Issue with test {test_id}"
     else:
@@ -202,7 +202,7 @@ def test_gfo_subdivide(test_descr, geom, subdivide_coords, exp_result):
     geom_wkb = None if geom is None else geom.wkb
 
     # Test
-    result_wkb = sqlite_userdefined.gfo_subdivide(geom_wkb, coords=subdivide_coords)
+    result_wkb = _sqlite_userdefined.gfo_subdivide(geom_wkb, coords=subdivide_coords)
 
     # Check result
     result = shapely.from_wkb(result_wkb)
diff --git a/tests/two_layer_operations/test_geofileops_twolayers.py b/tests/two_layer_operations/test_geofileops_twolayers.py
index 7b71997..1525587 100644
--- a/tests/two_layer_operations/test_geofileops_twolayers.py
+++ b/tests/two_layer_operations/test_geofileops_twolayers.py
@@ -615,24 +615,31 @@ def test_intersect_deprecated(tmp_path):
 
 
 @pytest.mark.parametrize(
-    "suffix, epsg, gridsize, explodecollections, nb_parallel",
+    "suffix_in, suffix_out, epsg, gridsize, explodecollections, nb_parallel",
     [
-        (".gpkg", 31370, 0.0, True, 1),
-        (".gpkg", 31370, 0.01, True, 1),
-        (".gpkg", 31370, 0.0, False, 2),
-        (".gpkg", 4326, 0.0, True, 2),
-        (".shp", 31370, 0.0, True, 1),
-        (".shp", 31370, 0.0, False, 2),
+        (".gpkg.zip", ".gpkg", 31370, 0.0, True, 1),
+        (".gpkg", ".gpkg", 31370, 0.01, True, 1),
+        (".gpkg", ".gpkg", 31370, 0.0, False, 2),
+        (".gpkg", ".gpkg", 4326, 0.0, True, 2),
+        (".shp", ".shp", 31370, 0.0, True, 1),
+        (".shp", ".shp", 31370, 0.0, False, 2),
     ],
 )
 def test_intersection(
-    tmp_path, suffix, epsg, explodecollections, gridsize, nb_parallel
+    tmp_path, suffix_in, suffix_out, epsg, explodecollections, gridsize, nb_parallel
 ):
+    if suffix_in == ".gpkg.zip":
+        pytest.xfail(
+            "Two layer operations use sqlite directly, so .gpkg.zip does not work"
+        )
+
     # Prepare test data/parameters
-    input1_path = test_helper.get_testfile("polygon-parcel", suffix=suffix, epsg=epsg)
-    input2_path = test_helper.get_testfile("polygon-zone", suffix=suffix, epsg=epsg)
+    input1_path = test_helper.get_testfile(
+        "polygon-parcel", suffix=suffix_in, epsg=epsg
+    )
+    input2_path = test_helper.get_testfile("polygon-zone", suffix=suffix_in, epsg=epsg)
 
-    output_path = tmp_path / f"{input1_path.stem}_inters_{input2_path.stem}{suffix}"
+    output_path = tmp_path / f"{input1_path.stem}_inters_{input2_path.stem}{suffix_out}"
     batchsize = -1
     input1_layerinfo = gfo.get_layerinfo(input1_path)
     if nb_parallel > 1:
@@ -659,7 +666,7 @@ def test_intersection(
         len(input1_layerinfo.columns) + len(input2_layerinfo.columns)
     )
 
-    if explodecollections and suffix != ".shp":
+    if explodecollections and suffix_out != ".shp":
         assert output_layerinfo.geometrytype == GeometryType.POLYGON
     else:
         assert output_layerinfo.geometrytype == GeometryType.MULTIPOLYGON