Repository: geofileops/geofileops
Pull request #595 (2024-12-11 22:48:52)

ENH: rename erase to difference
To use same terminology as most other open source projects...

codecov[bot]: ## [Codecov](https://app.codecov.io/gh/geofileops/geofileops/pull/595?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 `97.77778%` with `1 line` in your changes missing coverage. Please review.
> Project coverage is 90.00%. 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 [(`790278d`)](https://app.codecov.io/gh/geofileops/geofileops/commit/790278d245f0d6f2687354ecadbe8618046d8fa4?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=geofileops).
> Report is 78 commits behind head on main.

| [Files with missing lines](https://app.codecov.io/gh/geofileops/geofileops/pull/595?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/geoops.py](https://app.codecov.io/gh/geofileops/geofileops/pull/595?src=pr&el=tree&filepath=geofileops%2Fgeoops.py&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=geofileops#diff-Z2VvZmlsZW9wcy9nZW9vcHMucHk=) | 90.90% | [1 Missing :warning: ](https://app.codecov.io/gh/geofileops/geofileops/pull/595?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     #595      +/-   ##
==========================================
- Coverage   93.11%   90.00%   -3.12%     
==========================================
  Files          21       21              
  Lines        4113     4322     +209     
==========================================
+ Hits         3830     3890      +60     
- Misses        283      432     +149     
```

</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/geofileops/geofileops/pull/595?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 4dda1db..6f4ad02 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,6 +2,12 @@
 
 ## 0.10 (????-??-??)
 
+### Deprecations and compatibility notes
+
+- `erase` was renamed to `difference`, as most other open source applications/libraries
+  use this terminology. `erase` just keeps existing for backwards compatibility for now,
+  but a warning is shown that it might be removed in the (distant) future. (#595)
+
 ### Improvements
 
 - Add function `apply_vectorized` to apply a vectorized function on a geometry column
diff --git a/README.md b/README.md
index f65098a..a9becce 100644
--- a/README.md
+++ b/README.md
@@ -11,7 +11,7 @@ Geofileops is a python toolbox to process large vector files faster.
 Most typical GIS operations are available: e.g.
 [buffer](https://geofileops.readthedocs.io/en/stable/api/geofileops.apply.html#geofileops.buffer), 
 [dissolve](https://geofileops.readthedocs.io/en/stable/api/geofileops.apply.html#geofileops.dissolve),
-[erase](https://geofileops.readthedocs.io/en/stable/api/geofileops.apply.html#geofileops.erase)/difference, 
+[difference](https://geofileops.readthedocs.io/en/stable/api/geofileops.apply.html#geofileops.difference), 
 [intersection](https://geofileops.readthedocs.io/en/stable/api/geofileops.apply.html#geofileops.intersection), 
 [union](https://geofileops.readthedocs.io/en/stable/api/geofileops.apply.html#geofileops.union),... 
 
diff --git a/docs/reference.rst b/docs/reference.rst
index 8f31405..e6e4181 100644
--- a/docs/reference.rst
+++ b/docs/reference.rst
@@ -39,7 +39,7 @@ Spatial overlays and joins
    :toctree: api/
 
    clip
-   erase
+   difference
    export_by_distance
    export_by_location
    intersection
diff --git a/docs/user_guide.rst b/docs/user_guide.rst
index fd1470a..0796fab 100644
--- a/docs/user_guide.rst
+++ b/docs/user_guide.rst
@@ -107,7 +107,7 @@ Spatial overlays
 ----------------
 
 The standard :ref:`spatial overlays <reference-spatial-overlays-joins>` are
-available: :meth:`~intersection`, :meth:`~erase`, :meth:`~clip`, :meth:`~identity`,
+available: :meth:`~intersection`, :meth:`~difference`, :meth:`~clip`, :meth:`~identity`,
 :meth:`~union`, ...
 
 An example:
diff --git a/geofileops/geoops.py b/geofileops/geoops.py
index 9ce832d..31af3cf 100644
--- a/geofileops/geoops.py
+++ b/geofileops/geoops.py
@@ -200,9 +200,9 @@ def dissolve_within_distance(
         # temporary boundariesstep += 1
         logger.info(f"Step {step} of {nb_steps}")
         parts_to_add_path = tempdir / "200_parts_to_add.gpkg"
-        _geoops_sql.erase(
-            input_path=bufp_diss_bufm_path,
-            erase_path=diss_path,
+        _geoops_sql.difference(
+            input1_path=bufp_diss_bufm_path,
+            input2_path=diss_path,
             output_path=parts_to_add_path,
             overlay_self=False,
             explodecollections=True,
@@ -249,7 +249,7 @@ def dissolve_within_distance(
             #   - if > 1 neighbour, seems OK.
             #
             # For all pieces that don't comply to the above, the following parameters
-            # indicate that they need to be selected to erase them:
+            # indicate that they need to be selected to difference them:
             #   - pieces can be very narrow slivers. E.g. alongside a long boundary with
             #     a small bend, probably due to rounding side effects in the +/- buffer.
             #   - pieces can be spikes. E.g. when a "road" of ~ 'distance' width is not
@@ -1883,13 +1883,13 @@ def clip(
     )
 
 
-def erase(
-    input_path: Union[str, "os.PathLike[Any]"],
-    erase_path: Union[str, "os.PathLike[Any]", None],
+def difference(
+    input1_path: Union[str, "os.PathLike[Any]"],
+    input2_path: Union[str, "os.PathLike[Any]", None],
     output_path: Union[str, "os.PathLike[Any]"],
-    input_layer: Optional[str] = None,
-    input_columns: Optional[list[str]] = None,
-    erase_layer: Optional[str] = None,
+    input1_layer: Optional[str] = None,
+    input1_columns: Optional[list[str]] = None,
+    input2_layer: Optional[str] = None,
     output_layer: Optional[str] = None,
     explodecollections: bool = False,
     gridsize: float = 0.0,
@@ -1899,37 +1899,37 @@ def erase(
     subdivide_coords: int = 2000,
     force: bool = False,
 ):
-    """Erase all features in the erase layer from the features in the input layer.
+    """Calculate the difference of the input1 layer and input2 layer.
 
     Notes:
         - Every row in the input layer will result in maximum one row in the
           output layer.
-        - The output will contain the columns from the input layer, no columns from the
-          erase layer. The attribute values wont't be changed, so columns like area,...
+        - The output will contain the columns from the 1st no columns from the 2nd
+          layer. The attribute values wont't be changed, so columns like area,...
           will have to be recalculated manually.
-        - If ``erase_path`` is None, the 1st input layer is used for both inputs but
+        - If ``input2_path`` is None, the 1st input layer is used for both inputs but
           interactions between the same rows in this layer will be ignored. The output
           will be the (pieces of) features in this layer that don't have any
           intersections with other features in this layer.
 
     Alternative names:
-        - QGIS: difference
+        - ArcMap: erase
 
     Args:
-        input_path (PathLike): The file to erase from.
-        erase_path (PathLike, optional): The file with the geometries to erase with. If
-            None, the 1st input layer is used for both inputs but interactions between
-            the same rows in this layer will be ignored. The output will be the (pieces
-            of) features in this layer that don't have any intersections with other
-            features in this layer.
+        input1_path (PathLike): The file to remove/difference from.
+        input2_path (PathLike, optional): The file with the geometries to remove from
+            input1. If None, the 1st input layer is used for both inputs but interactions
+            between the same rows in this layer will be ignored. The output will be the
+            (pieces of) features in this layer that don't have any intersections with
+            other features in this layer.
         output_path (PathLike): the file to write the result to.
-        input_layer (str, optional): input layer name. If None, ``input_path`` should
+        input1_layer (str, optional): input layer name. If None, ``input1_path`` should
             contain only one layer. Defaults to None.
-        input_columns (List[str], optional): list of columns to retain. If None, all
+        input1_columns (List[str], optional): list of columns to retain. If None, all
             standard columns are retained. In addition to standard columns, it is also
             possible to specify "fid", a unique index available in all input files. Note
             that the "fid" will be aliased eg. to "fid_1". Defaults to None.
-        erase_layer (str, optional): erase layer name. If None, ``erase_path`` should
+        input2_layer (str, optional): input2 layer name. If None, ``input2_path`` should
             contain only one layer. Defaults to None.
         output_layer (str, optional): output layer name. If None, the ``output_path``
             stem is used. Defaults to None.
@@ -1960,26 +1960,67 @@ def erase(
         <a href="https://www.gaia-gis.it/gaia-sins/spatialite-sql-latest.html" target="_blank">spatialite reference</a>
 
     """  # noqa: E501
-    logger = logging.getLogger("geofileops.erase")
-    logger.info(f"Start, on {input_path} with {erase_path} to {output_path}")
+    logger = logging.getLogger("geofileops.difference")
+    logger.info(f"Start, on {input1_path} with {input2_path} to {output_path}")
 
-    # If erase_path is None, we are doing a self-overlay
+    # If input2_path is None, we are doing a self-overlay
     overlay_self = False
-    if erase_path is None:
-        if erase_layer is not None:
-            raise ValueError("erase_layer must be None if erase_path is None")
-        erase_path = input_path
-        erase_layer = input_layer
+    if input2_path is None:
+        if input2_layer is not None:
+            raise ValueError("input2_layer must be None if input2_path is None")
+        input2_path = input1_path
+        input2_layer = input1_layer
         overlay_self = True
 
-    return _geoops_sql.erase(
-        input_path=Path(input_path),
-        erase_path=Path(erase_path),
+    return _geoops_sql.difference(
+        input1_path=Path(input1_path),
+        input2_path=Path(input2_path),
         output_path=Path(output_path),
         overlay_self=overlay_self,
-        input_layer=input_layer,
-        input_columns=input_columns,
-        erase_layer=erase_layer,
+        input_layer=input1_layer,
+        input1_columns=input1_columns,
+        input2_layer=input2_layer,
+        output_layer=output_layer,
+        explodecollections=explodecollections,
+        gridsize=gridsize,
+        where_post=where_post,
+        nb_parallel=nb_parallel,
+        batchsize=batchsize,
+        subdivide_coords=subdivide_coords,
+        force=force,
+    )
+
+
+def erase(
+    input_path: Union[str, "os.PathLike[Any]"],
+    erase_path: Union[str, "os.PathLike[Any]", None],
+    output_path: Union[str, "os.PathLike[Any]"],
+    input_layer: Optional[str] = None,
+    input_columns: Optional[list[str]] = None,
+    erase_layer: Optional[str] = None,
+    output_layer: Optional[str] = None,
+    explodecollections: bool = False,
+    gridsize: float = 0.0,
+    where_post: Optional[str] = None,
+    nb_parallel: int = -1,
+    batchsize: int = -1,
+    subdivide_coords: int = 2000,
+    force: bool = False,
+):
+    """DEPRECATED: please use difference."""
+    warnings.warn(  # pragma: no cover
+        "erase is deprecated because it was renamed to difference. "
+        "Will be removed in a (distant) future version",
+        FutureWarning,
+        stacklevel=2,
+    )
+    return difference(
+        input1_path=input_path,
+        input2_path=erase_path,
+        output_path=output_path,
+        input1_layer=input_layer,
+        input1_columns=input_columns,
+        input2_layer=erase_layer,
         output_layer=output_layer,
         explodecollections=explodecollections,
         gridsize=gridsize,
@@ -2209,7 +2250,7 @@ def identity(
     r"""Calculates the pairwise identity of the two input layers.
 
     The result is the equivalent of the intersection between the two layers + layer 1
-    erased with layer 2.
+    differenced with layer 2.
 
     Notes:
         - The result will contain the attribute columns from both input layers. The
@@ -2372,7 +2413,7 @@ def intersect(
 ):
     """DEPRECATED: please use intersection."""
     warnings.warn(  # pragma: no cover
-        "intersect() is deprecated because it was renamed intersection(). "
+        "intersect is deprecated because it was renamed intersection. "
         "Will be removed in a future version",
         FutureWarning,
         stacklevel=2,
diff --git a/geofileops/util/_geoops_sql.py b/geofileops/util/_geoops_sql.py
index 5595107..ed4da1e 100644
--- a/geofileops/util/_geoops_sql.py
+++ b/geofileops/util/_geoops_sql.py
@@ -898,7 +898,7 @@ def clip(
     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 clip cause eg. polygons to be split to multipolygons.
+    # type to multi, because clip 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:
         force_output_geometrytype = force_output_geometrytype.to_multitype
@@ -979,14 +979,14 @@ def clip(
     )
 
 
-def erase(
-    input_path: Path,
-    erase_path: Path,
+def difference(
+    input1_path: Path,
+    input2_path: Path,
     output_path: Path,
     overlay_self: bool,
     input_layer: Optional[str] = None,
-    input_columns: Optional[list[str]] = None,
-    erase_layer: Optional[str] = None,
+    input1_columns: Optional[list[str]] = None,
+    input2_layer: Optional[str] = None,
     output_layer: Optional[str] = None,
     explodecollections: bool = False,
     gridsize: float = 0.0,
@@ -999,25 +999,25 @@ def erase(
     output_with_spatial_index: Optional[bool] = None,
     operation_prefix: str = "",
 ):
-    # Because there might be extra preparation of the erase layer before going ahead
+    # Because there might be extra preparation of the input2 layer before going ahead
     # with the real calculation, do some additional init + checks here...
     if subdivide_coords < 0:
         raise ValueError("subdivide_coords < 0 is not allowed")
 
-    operation_name = f"{operation_prefix}erase"
+    operation_name = f"{operation_prefix}difference"
     logger = logging.getLogger(f"geofileops.{operation_name}")
 
     # Get layer names
     if input_layer is None:
-        input_layer = gfo.get_only_layer(input_path)
-    if erase_layer is None:
-        erase_layer = gfo.get_only_layer(erase_path)
+        input_layer = gfo.get_only_layer(input1_path)
+    if input2_layer is None:
+        input2_layer = gfo.get_only_layer(input2_path)
 
     if _io_util.output_exists(path=output_path, remove_if_exists=force):
         return
 
     start_time = datetime.now()
-    input_layer_info = gfo.get_layerinfo(input_path, input_layer)
+    input_layer_info = gfo.get_layerinfo(input1_path, input_layer)
     primitivetypeid = input_layer_info.geometrytype.to_primitivetype.value
 
     force_output_geometrytype = input_layer_info.geometrytype
@@ -1025,7 +1025,7 @@ def erase(
         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
+        # output type to multi, because difference can cause eg. polygons to be split to
         # multipolygons.
         force_output_geometrytype = force_output_geometrytype.to_multitype
 
@@ -1033,8 +1033,8 @@ def erase(
     # Save the original fid column in a new fid_1 column, we will need it to filter on
     # it later on.
     tempdir = _io_util.create_tempdir(f"geofileops/{operation_name}")
-    input_subdivided_path = _subdivide_layer(
-        path=input_path,
+    input1_subdivided_path = _subdivide_layer(
+        path=input1_path,
         layer=input_layer,
         output_path=tempdir / "subdivided/input_layer.gpkg",
         subdivide_coords=subdivide_coords,
@@ -1047,22 +1047,22 @@ def erase(
     where_clause_self = "1=1"
     if overlay_self:
         # If we are doing a self overlay, we need to filter out rows with the same rowid
-        if input_subdivided_path is None:
+        if input1_subdivided_path is None:
             where_clause_self = "layer1.rowid <> layer2_sub.rowid"
         else:
             # Filter out the same rowids using the original fids!
             where_clause_self = "layer1_subdiv.fid_1 <> layer2_sub.fid_1"
 
         # For overlay self, both subdivided layers are equal
-        erase_subdivided_path = input_subdivided_path
+        input2_subdivided_path = input1_subdivided_path
 
     else:
-        # Subdivide the erase layer if needed to speed up further processing.
-        # No self overlay, so we don't need the original fid column for the erase layer.
-        erase_subdivided_path = _subdivide_layer(
-            path=erase_path,
-            layer=erase_layer,
-            output_path=tempdir / "subdivided/erase_layer.gpkg",
+        # Subdivide the input2 layer if needed to speed up further processing.
+        # No self overlay, so no original fid column for the input2 layer needed.
+        input2_subdivided_path = _subdivide_layer(
+            path=input2_path,
+            layer=input2_layer,
+            output_path=tempdir / "subdivided/input2_layer.gpkg",
             subdivide_coords=subdivide_coords,
             keep_fid=False,
             nb_parallel=nb_parallel,
@@ -1070,9 +1070,9 @@ def erase(
             operation_prefix=f"{operation_name}/",
         )
 
-    # If the erase layer was subdivided
-    if erase_subdivided_path is not None:
-        erase_path = erase_subdivided_path
+    # If the input2 layer was subdivided
+    if input2_subdivided_path is not None:
+        input2_path = input2_subdivided_path
 
     # Prepare sql template for this operation
     # - WHERE geom IS NOT NULL to avoid rows with a NULL geom, they give issues in
@@ -1096,7 +1096,7 @@ def erase(
     input2_layer_rtree = "rtree_{input2_layer}_{input2_geometrycolumn}"
     input1_subdiv_layer_rtree = "rtree_{input1_layer}_{input1_subdiv_geometrycolumn}"
 
-    if input_subdivided_path is None:
+    if input1_subdivided_path is None:
         # The input layer was not subdivided
         sql_template = f"""
             SELECT * FROM (
@@ -1197,16 +1197,16 @@ def erase(
 
     # Go!
     _two_layer_vector_operation(
-        input1_path=input_path,
-        input1_subdivided_path=input_subdivided_path,
-        input2_path=erase_path,
+        input1_path=input1_path,
+        input1_subdivided_path=input1_subdivided_path,
+        input2_path=input2_path,
         output_path=output_path,
         sql_template=sql_template,
         operation_name=operation_name,
         input1_layer=input_layer,
-        input1_columns=input_columns,
+        input1_columns=input1_columns,
         input1_columns_prefix=input_columns_prefix,
-        input2_layer=erase_layer,
+        input2_layer=input2_layer,
         input2_columns=[],
         input2_columns_prefix="",
         output_layer=output_layer,
@@ -1222,7 +1222,7 @@ def erase(
     )
 
     # Print time taken
-    logger.info(f"Ready, full erase took {datetime.now()-start_time}")
+    logger.info(f"Ready, full difference took {datetime.now()-start_time}")
 
 
 def _subdivide_layer(
@@ -1382,7 +1382,7 @@ def export_by_location(
     subdivide_coords: int = 10000,
     force: bool = False,
 ):
-    # Because there might be extra preparation of the erase layer before going ahead
+    # Because there might be extra preparation of the 2nd layer before going ahead
     # with the real calculation, do some additional init + checks here...
     if subdivide_coords < 0:
         raise ValueError("subdivide_coords < 0 is not allowed")
@@ -1409,7 +1409,7 @@ def export_by_location(
         true_for_disjoint,
     ) = _prepare_filter_by_location_fields(spatial_relations_query)
 
-    # Subdivide the erase layer if applicable to speed up further processing.
+    # Subdivide the 2nd layer if applicable to speed up further processing.
     tmp_dir = _io_util.create_tempdir(f"geofileops/{operation_name}")
     input_to_compare_with_subdivided_path = _subdivide_layer(
         path=input_to_compare_with_path,
@@ -2265,7 +2265,7 @@ def identity(
     force: bool = False,
 ):
     # An identity is the combination of the results of an "intersection" of input1 and
-    # input2 and an erase of input2 with input1.
+    # input2 and an difference of input2 with input1.
 
     # Because the calculations of the intermediate results will be towards temp files,
     # we need to do some additional init + checks here...
@@ -2306,18 +2306,18 @@ def identity(
             operation_prefix="identity/",
         )
 
-        # Now erase input1 from input2 to another temporary output gfo...
-        logger.info("Step 2 of 3: erase")
-        erase_output_path = tempdir / "erase_output.gpkg"
-        erase(
-            input_path=input1_path,
-            erase_path=input2_path,
-            output_path=erase_output_path,
+        # Now difference input1 from input2 to another temporary output gfo...
+        logger.info("Step 2 of 3: difference")
+        difference_output_path = tempdir / "difference_output.gpkg"
+        difference(
+            input1_path=input1_path,
+            input2_path=input2_path,
+            output_path=difference_output_path,
             overlay_self=overlay_self,
             input_layer=input1_layer,
-            input_columns=input1_columns,
+            input1_columns=input1_columns,
             input_columns_prefix=input1_columns_prefix,
-            erase_layer=input2_layer,
+            input2_layer=input2_layer,
             output_layer=output_layer,
             explodecollections=explodecollections,
             gridsize=gridsize,
@@ -2334,7 +2334,7 @@ def identity(
         logger.info("Step 3 of 3: finalize")
         # Note: append will never create an index on an already existing layer.
         _append_to_nolock(
-            src=erase_output_path,
+            src=difference_output_path,
             dst=intersection_output_path,
             src_layer=output_layer,
             dst_layer=output_layer,
@@ -2379,11 +2379,11 @@ def symmetric_difference(
     subdivide_coords: int = 2000,
     force: bool = False,
 ):
-    # A symmetric difference can be simulated by doing an "erase" of input1
-    # and input2 and then append the result of an erase of input2 with
+    # A symmetric difference can be simulated by doing an difference of input1
+    # and input2 and then append the result of an difference of input2 with
     # input1...
 
-    # Because both erase calculations will be towards temp files,
+    # Because both difference calculations will be towards temp files,
     # we need to do some additional init + checks here...
     if subdivide_coords < 0:
         raise ValueError("subdivide_coords < 0 is not allowed")
@@ -2401,18 +2401,18 @@ def symmetric_difference(
 
     tempdir = _io_util.create_tempdir("geofileops/symmdiff")
     try:
-        # First erase input2 from input1 to a temporary output file
-        logger.info("Step 1 of 3: erase 1")
-        erase1_output_path = tempdir / "layer1_erase_layer2_output.gpkg"
-        erase(
-            input_path=input1_path,
-            erase_path=input2_path,
-            output_path=erase1_output_path,
+        # First difference input2 from input1 to a temporary output file
+        logger.info("Step 1 of 3: difference 1")
+        diff1_output_path = tempdir / "layer1_diff_layer2_output.gpkg"
+        difference(
+            input1_path=input1_path,
+            input2_path=input2_path,
+            output_path=diff1_output_path,
             overlay_self=overlay_self,
             input_layer=input1_layer,
-            input_columns=input1_columns,
+            input1_columns=input1_columns,
             input_columns_prefix=input1_columns_prefix,
-            erase_layer=input2_layer,
+            input2_layer=input2_layer,
             output_layer=output_layer,
             explodecollections=explodecollections,
             gridsize=gridsize,
@@ -2432,23 +2432,23 @@ def symmetric_difference(
             )
             for column in columns_to_add:
                 gfo.add_column(
-                    erase1_output_path,
+                    diff1_output_path,
                     name=f"{input2_columns_prefix}{column}",
                     type=input2_info.columns[column].gdal_type,
                 )
 
-        # Now erase input1 from input2 to another temporary output file
-        logger.info("Step 2 of 3: erase 2")
-        erase2_output_path = tempdir / "layer2_erase_layer1_output.gpkg"
-        erase(
-            input_path=input2_path,
-            erase_path=input1_path,
-            output_path=erase2_output_path,
+        # Now difference input1 from input2 to another temporary output file
+        logger.info("Step 2 of 3: difference 2")
+        diff2_output_path = tempdir / "layer2_diff_layer1_output.gpkg"
+        difference(
+            input1_path=input2_path,
+            input2_path=input1_path,
+            output_path=diff2_output_path,
             overlay_self=overlay_self,
             input_layer=input2_layer,
-            input_columns=input2_columns,
+            input1_columns=input2_columns,
             input_columns_prefix=input2_columns_prefix,
-            erase_layer=input1_layer,
+            input2_layer=input1_layer,
             output_layer=output_layer,
             explodecollections=explodecollections,
             gridsize=gridsize,
@@ -2465,18 +2465,18 @@ def symmetric_difference(
         logger.info("Step 3 of 3: finalize")
         # Note: append will never create an index on an already existing layer.
         _append_to_nolock(
-            src=erase2_output_path,
-            dst=erase1_output_path,
+            src=diff2_output_path,
+            dst=diff1_output_path,
             src_layer=output_layer,
             dst_layer=output_layer,
         )
 
         # Convert or add spatial index
-        tmp_output_path = erase1_output_path
-        if erase1_output_path.suffix != output_path.suffix:
+        tmp_output_path = diff1_output_path
+        if diff1_output_path.suffix != output_path.suffix:
             # Output file should be in diffent format, so convert
             tmp_output_path = tempdir / output_path.name
-            gfo.copy_layer(src=erase1_output_path, dst=tmp_output_path)
+            gfo.copy_layer(src=diff1_output_path, dst=tmp_output_path)
         elif GeofileInfo(tmp_output_path).default_spatial_index:
             gfo.create_spatial_index(path=tmp_output_path, layer=output_layer)
 
@@ -2511,7 +2511,8 @@ def union(
     force: bool = False,
 ):
     # A union is the combination of the results of an intersection of input1 and input2,
-    # the result of an erase of input2 with input1 and the erase of input1 with input2.
+    # the result of an difference of input2 with input1 and the difference of input1
+    # with input2.
 
     # Because the calculations of the intermediate results will be towards temp files,
     # we need to do some additional init + checks here...
@@ -2554,18 +2555,18 @@ def union(
             operation_prefix="union/",
         )
 
-        # Erase input1 from input2 to another temporary output gfo.
-        logger.info("Step 2 of 4: erase input 1 from input 2")
-        erase1_output_path = tempdir / "erase_input1_from_input2_output.gpkg"
-        erase(
-            input_path=input2_path,
-            erase_path=input1_path,
-            output_path=erase1_output_path,
+        # Difference input1 from input2 to another temporary output gfo.
+        logger.info("Step 2 of 4: difference of input 1 from input 2")
+        diff1_output_path = tempdir / "diff_input1_from_input2_output.gpkg"
+        difference(
+            input1_path=input2_path,
+            input2_path=input1_path,
+            output_path=diff1_output_path,
             overlay_self=overlay_self,
             input_layer=input2_layer,
-            input_columns=input2_columns,
+            input1_columns=input2_columns,
             input_columns_prefix=input2_columns_prefix,
-            erase_layer=input1_layer,
+            input2_layer=input1_layer,
             output_layer=output_layer,
             explodecollections=explodecollections,
             gridsize=gridsize,
@@ -2579,26 +2580,26 @@ def union(
         )
         # Note: append will never create an index on an already existing layer.
         _append_to_nolock(
-            src=erase1_output_path,
+            src=diff1_output_path,
             dst=intersection_output_path,
             src_layer=output_layer,
             dst_layer=output_layer,
         )
-        gfo.remove(erase1_output_path)
+        gfo.remove(diff1_output_path)
 
-        # Erase input1 from input2 to and add to temporary output file.
-        logger.info("Step 3 of 4: erase input 2 from input 1")
-        erase2_output_path = tempdir / "erase_input2_from_input1_output.gpkg"
+        # Difference input1 from input2 to and add to temporary output file.
+        logger.info("Step 3 of 4: difference input 2 from input 1")
+        diff2_output_path = tempdir / "diff_input2_from_input1_output.gpkg"
 
-        erase(
-            input_path=input1_path,
-            erase_path=input2_path,
-            output_path=erase2_output_path,
+        difference(
+            input1_path=input1_path,
+            input2_path=input2_path,
+            output_path=diff2_output_path,
             overlay_self=overlay_self,
             input_layer=input1_layer,
-            input_columns=input1_columns,
+            input1_columns=input1_columns,
             input_columns_prefix=input1_columns_prefix,
-            erase_layer=input2_layer,
+            input2_layer=input2_layer,
             output_layer=output_layer,
             explodecollections=explodecollections,
             gridsize=gridsize,
@@ -2611,12 +2612,12 @@ def union(
             operation_prefix="union/",
         )
         _append_to_nolock(
-            src=erase2_output_path,
+            src=diff2_output_path,
             dst=intersection_output_path,
             src_layer=output_layer,
             dst_layer=output_layer,
         )
-        gfo.remove(erase2_output_path)
+        gfo.remove(diff2_output_path)
 
         # Convert or add spatial index
         logger.info("Step 4 of 4: finalize")


Test patch:
diff --git a/ci/envs/latest-fiona.yml b/ci/envs/latest-fiona.yml
index f949a19..8911dba 100644
--- a/ci/envs/latest-fiona.yml
+++ b/ci/envs/latest-fiona.yml
@@ -1,4 +1,4 @@
-name: geofileops-latest-nofiona
+name: geofileops-latest-fiona
 channels:
   - conda-forge
 dependencies:
diff --git a/tests/test_geofileops_twolayers.py b/tests/test_geofileops_twolayers.py
index 2760768..b1a8f87 100644
--- a/tests/test_geofileops_twolayers.py
+++ b/tests/test_geofileops_twolayers.py
@@ -116,7 +116,7 @@ def test_clip_resultempty(tmp_path, suffix, clip_empty):
     not GEOPANDAS_GTE_10,
     reason="assert_geodataframe_equal with check_geom_gridsize requires gpd >= 1.0",
 )
-def test_erase(
+def test_difference(
     tmp_path,
     suffix,
     testfile,
@@ -125,25 +125,25 @@ def test_erase(
     subdivide_coords,
     check_geom_gridsize,
 ):
-    input_path = test_helper.get_testfile(testfile, suffix=suffix)
+    input1_path = test_helper.get_testfile(testfile, suffix=suffix)
     if suffix == ".shp":
-        erase_path = test_helper.get_testfile("polygon-zone", suffix=suffix)
-        erase_layer = None
+        input2_path = test_helper.get_testfile("polygon-zone", suffix=suffix)
+        input2_layer = None
     else:
-        erase_path = test_helper.get_testfile("polygon-twolayers", suffix=suffix)
-        erase_layer = "zones"
-    input_layerinfo = gfo.get_layerinfo(input_path)
+        input2_path = test_helper.get_testfile("polygon-twolayers", suffix=suffix)
+        input2_layer = "zones"
+    input_layerinfo = gfo.get_layerinfo(input1_path)
     batchsize = math.ceil(input_layerinfo.featurecount / 2)
-    output_path = tmp_path / f"{input_path.stem}-output{suffix}"
+    output_path = tmp_path / f"{input1_path.stem}-output{suffix}"
 
     kwargs = {}
     if subdivide_coords is not None:
         kwargs["subdivide_coords"] = subdivide_coords
 
-    gfo.erase(
-        input_path=str(input_path),
-        erase_path=str(erase_path),
-        erase_layer=erase_layer,
+    gfo.difference(
+        input1_path=str(input1_path),
+        input2_path=str(input2_path),
+        input2_layer=input2_layer,
         output_path=str(output_path),
         gridsize=gridsize,
         where_post=where_post,
@@ -156,11 +156,11 @@ def test_erase(
     exp_spatial_index = GeofileInfo(output_path).default_spatial_index
     assert gfo.has_spatial_index(output_path) is exp_spatial_index
     output_gdf = gfo.read_file(output_path)
-    input_gdf = gfo.read_file(input_path)
-    erase_gdf = gfo.read_file(erase_path, layer=erase_layer)
+    input1_gdf = gfo.read_file(input1_path)
+    input2_gdf = gfo.read_file(input2_path, layer=input2_layer)
 
     # Prepare expected gdf
-    exp_gdf = gpd.overlay(input_gdf, erase_gdf, how="difference", keep_geom_type=True)
+    exp_gdf = gpd.overlay(input1_gdf, input2_gdf, how="difference", keep_geom_type=True)
     if gridsize != 0.0:
         exp_gdf.geometry = shapely.set_precision(exp_gdf.geometry, grid_size=gridsize)
     if where_post is not None:
@@ -175,7 +175,7 @@ def test_erase(
     exp_gdf = exp_gdf[~exp_gdf.geometry.is_empty]
 
     if test_helper.RUNS_LOCAL:
-        output_exp_path = tmp_path / f"{input_path.stem}-expected{suffix}"
+        output_exp_path = tmp_path / f"{input1_path.stem}-expected{suffix}"
         gfo.to_file(exp_gdf, output_exp_path)
 
     assert_geodataframe_equal(
@@ -192,16 +192,16 @@ def test_erase(
     assert len(output_gdf) > 0
 
 
-def test_erase_explodecollections(tmp_path):
-    input_path = test_helper.get_testfile("polygon-parcel")
-    erase_path = test_helper.get_testfile("polygon-zone")
-    input_layerinfo = gfo.get_layerinfo(input_path)
+def test_difference_explodecollections(tmp_path):
+    input1_path = test_helper.get_testfile("polygon-parcel")
+    input2_path = test_helper.get_testfile("polygon-zone")
+    input_layerinfo = gfo.get_layerinfo(input1_path)
     batchsize = math.ceil(input_layerinfo.featurecount / 2)
 
-    output_path = tmp_path / f"{input_path.stem}-output_exploded{input_path.suffix}"
-    gfo.erase(
-        input_path=input_path,
-        erase_path=erase_path,
+    output_path = tmp_path / f"{input1_path.stem}-output_exploded{input1_path.suffix}"
+    gfo.difference(
+        input1_path=input1_path,
+        input2_path=input2_path,
         output_path=output_path,
         explodecollections=True,
         batchsize=batchsize,
@@ -212,10 +212,10 @@ def test_erase_explodecollections(tmp_path):
     exp_spatial_index = GeofileInfo(output_path).default_spatial_index
     assert gfo.has_spatial_index(output_path) is exp_spatial_index
     output_gdf = gfo.read_file(output_path)
-    input_gdf = gfo.read_file(input_path)
-    erase_gdf = gfo.read_file(erase_path)
+    input1_gdf = gfo.read_file(input1_path)
+    input2_gdf = gfo.read_file(input2_path)
     output_gpd_gdf = gpd.overlay(
-        input_gdf, erase_gdf, how="difference", keep_geom_type=True
+        input1_gdf, input2_gdf, how="difference", keep_geom_type=True
     )
     output_gpd_gdf = output_gpd_gdf.explode(ignore_index=True)
     assert_geodataframe_equal(
@@ -228,26 +228,26 @@ def test_erase_explodecollections(tmp_path):
     )
 
 
-def test_erase_force(tmp_path):
+def test_difference_force(tmp_path):
     # Prepare test data
-    input_path = test_helper.get_testfile("polygon-parcel")
-    erase_path = test_helper.get_testfile("polygon-zone")
-    output_path = tmp_path / f"output{input_path.suffix}"
+    input1_path = test_helper.get_testfile("polygon-parcel")
+    input2_path = test_helper.get_testfile("polygon-zone")
+    output_path = tmp_path / f"output{input1_path.suffix}"
     output_path.touch()
 
     # Test with force False (the default): existing output file should stay the same
     mtime_orig = output_path.stat().st_mtime
-    gfo.erase(
-        input_path=input_path,
-        erase_path=erase_path,
+    gfo.difference(
+        input1_path=input1_path,
+        input2_path=input2_path,
         output_path=output_path,
     )
     assert output_path.stat().st_mtime == mtime_orig
 
     # With force=True
-    gfo.erase(
-        input_path=input_path,
-        erase_path=erase_path,
+    gfo.difference(
+        input1_path=input1_path,
+        input2_path=input2_path,
         output_path=output_path,
         force=True,
     )
@@ -257,33 +257,33 @@ def test_erase_force(tmp_path):
 @pytest.mark.parametrize(
     "kwargs, expected_error",
     [
-        ({"erase_path": None}, "erase_layer must be None if erase_path is None"),
+        ({"input2_path": None}, "input2_layer must be None if input2_path is None"),
         ({"subdivide_coords": -1}, "subdivide_coords < 0 is not allowed"),
     ],
 )
-def test_erase_invalid_params(kwargs, expected_error):
-    if "erase_path" not in kwargs:
-        kwargs["erase_path"] = "erase.gpkg"
+def test_difference_invalid_params(kwargs, expected_error):
+    if "input2_path" not in kwargs:
+        kwargs["input2_path"] = "input2.gpkg"
     with pytest.raises(ValueError, match=expected_error):
-        gfo.erase(
-            input_path="input.gpkg",
+        gfo.difference(
+            input1_path="input.gpkg",
             output_path="output.gpkg",
-            erase_layer="invalid",
+            input2_layer="invalid",
             **kwargs,
         )
 
 
 @pytest.mark.parametrize("subdivide_coords", [2000, 5])
-def test_erase_self(tmp_path, subdivide_coords):
-    input_path = test_helper.get_testfile("polygon-overlappingcircles-all")
-    input_layerinfo = gfo.get_layerinfo(input_path)
+def test_difference_self(tmp_path, subdivide_coords):
+    input1_path = test_helper.get_testfile("polygon-overlappingcircles-all")
+    input_layerinfo = gfo.get_layerinfo(input1_path)
     batchsize = math.ceil(input_layerinfo.featurecount / 2)
 
     # Now run test
-    output_path = tmp_path / f"{input_path.stem}_erase_self.gpkg"
-    gfo.erase(
-        input_path=input_path,
-        erase_path=None,
+    output_path = tmp_path / f"{input1_path.stem}_diff_self.gpkg"
+    gfo.difference(
+        input1_path=input1_path,
+        input2_path=None,
         output_path=output_path,
         subdivide_coords=subdivide_coords,
         nb_parallel=2,
@@ -300,27 +300,28 @@ def test_erase_self(tmp_path, subdivide_coords):
 
 
 @pytest.mark.parametrize("suffix", SUFFIXES_GEOOPS)
-def test_erase_subdivide_multipolygons(tmp_path, suffix):
+def test_difference_subdivide_multipolygons(tmp_path, suffix):
     """
-    Test if erase with subdivide also works if the erase layer contains multipolygons.
+    Test if difference with subdivide also works if the input2 layer contains
+    multipolygons.
 
     It seems spatialite function ST_AsBinary, ST_GeomFromWKB and/or ST_Collect have
     issues processing nested multi-types (e.g. a GeometryCollection containing e.g.
     MultiPolygons).
     """
     # Prepare test data
-    input_path = test_helper.get_testfile("point", suffix=suffix)
-    input_layerinfo = gfo.get_layerinfo(input_path)
-    batchsize = math.ceil(input_layerinfo.featurecount / 2)
+    input1_path = test_helper.get_testfile("point", suffix=suffix)
+    input1_layerinfo = gfo.get_layerinfo(input1_path)
+    batchsize = math.ceil(input1_layerinfo.featurecount / 2)
 
-    # Prepare erase test data: should be multipolygons for good test coverage
+    # Prepare test data: should be multipolygons for good test coverage
     zone_path = test_helper.get_testfile("polygon-zone", suffix=suffix)
     zones_gdf = gfo.read_file(zone_path).explode(ignore_index=True)
 
-    erase_geometries = [
-        {"desc": "erase1", "geometry": zones_gdf.geometry[4]},
+    input2_geometries = [
+        {"desc": "input2_1", "geometry": zones_gdf.geometry[4]},
         {
-            "desc": "erase2",
+            "desc": "input2_2",
             "geometry": sh_geom.MultiPolygon(
                 [
                     zones_gdf.geometry[0],
@@ -331,14 +332,14 @@ def test_erase_subdivide_multipolygons(tmp_path, suffix):
             ),
         },
     ]
-    erase_gdf = gpd.GeoDataFrame(erase_geometries, crs=31370)
-    erase_path = tmp_path / f"{zone_path.stem}_multi{suffix}"
-    gfo.to_file(erase_gdf, erase_path)
+    input2_gdf = gpd.GeoDataFrame(input2_geometries, crs=31370)
+    input2_path = tmp_path / f"{zone_path.stem}_multi{suffix}"
+    gfo.to_file(input2_gdf, input2_path)
 
-    output_path = tmp_path / f"{input_path.stem}-output_exploded{suffix}"
-    gfo.erase(
-        input_path=input_path,
-        erase_path=erase_path,
+    output_path = tmp_path / f"{input1_path.stem}-output_exploded{suffix}"
+    gfo.difference(
+        input1_path=input1_path,
+        input2_path=input2_path,
         output_path=output_path,
         batchsize=batchsize,
         subdivide_coords=10,
@@ -349,9 +350,9 @@ def test_erase_subdivide_multipolygons(tmp_path, suffix):
     exp_spatial_index = GeofileInfo(output_path).default_spatial_index
     assert gfo.has_spatial_index(output_path) is exp_spatial_index
     output_gdf = gfo.read_file(output_path)
-    input_gdf = gfo.read_file(input_path)
+    input1_gdf = gfo.read_file(input1_path)
     output_gpd_gdf = gpd.overlay(
-        input_gdf, erase_gdf, how="difference", keep_geom_type=True
+        input1_gdf, input2_gdf, how="difference", keep_geom_type=True
     )
     output_gpd_gdf = output_gpd_gdf.explode(ignore_index=True)
     assert_geodataframe_equal(
@@ -364,6 +365,30 @@ def test_erase_subdivide_multipolygons(tmp_path, suffix):
     )
 
 
+def test_erase_deprecated(tmp_path):
+    """Minimal test of the deprecated erase function."""
+    input1_path = test_helper.get_testfile("polygon-parcel")
+    input2_path = test_helper.get_testfile("polygon-twolayers")
+    input2_layer = "zones"
+    output_path = tmp_path / f"{input1_path.stem}-output.gpkg"
+
+    with pytest.warns(FutureWarning, match="erase is deprecated"):
+        gfo.erase(
+            input_path=str(input1_path),
+            erase_path=str(input2_path),
+            input_layer="parcels",
+            input_columns=[],
+            erase_layer=input2_layer,
+            output_path=str(output_path),
+            gridsize=0.0,
+            where_post=None,
+            batchsize=-1,
+        )
+
+    # Compare result with geopandas
+    assert output_path.exists()
+
+
 @pytest.mark.parametrize("suffix", SUFFIXES_GEOOPS)
 @pytest.mark.parametrize(
     "columns, gridsize, where_post, subdivide_coords, exp_featurecount",
@@ -736,6 +761,29 @@ def test_identity_self(tmp_path):
     assert output_layerinfo.featurecount == 9
 
 
+def test_intersect_deprecated(tmp_path):
+    """Minimal test of the deprecated intersect function."""
+    input1_path = test_helper.get_testfile("polygon-parcel")
+    input2_path = test_helper.get_testfile("polygon-zone")
+
+    # Now run test
+    output_path = tmp_path / f"{input1_path.stem}_intersection_{input2_path.stem}.gpkg"
+
+    with pytest.warns(FutureWarning, match="intersect is deprecated"):
+        gfo.intersect(
+            input1_path=str(input1_path),
+            input2_path=str(input2_path),
+            output_path=str(output_path),
+            gridsize=0.01,
+            explodecollections=True,
+            nb_parallel=2,
+            batchsize=-1,
+        )
+
+    # Check if the tmp file is correctly created
+    assert output_path.exists()
+
+
 @pytest.mark.parametrize("testfile", ["polygon-parcel"])
 @pytest.mark.parametrize(
     "suffix, epsg, gridsize, explodecollections, nb_parallel",
@@ -1811,59 +1859,24 @@ def test_select_two_layers_select_star_fids_unique(tmp_path, suffix):
     assert len(output_layerinfo.columns) == 1
 
 
-@pytest.mark.filterwarnings(
-    "ignore: split is deprecated because it was renamed to identity"
-)
-def test_split(tmp_path):
-    """Is deprecated, but keep minimal test."""
+def test_split_deprecated(tmp_path):
+    """Minimal tests for the deprecated split function."""
     # Prepare test data
     input1_path = test_helper.get_testfile("polygon-parcel")
     input2_path = test_helper.get_testfile("polygon-zone")
-    input1_layerinfo = gfo.get_layerinfo(input1_path)
-    batchsize = math.ceil(input1_layerinfo.featurecount / 2)
     output_path = tmp_path / f"{input1_path.stem}-output.gpkg"
 
     # Test
-    gfo.split(
-        input1_path=str(input1_path),
-        input2_path=str(input2_path),
-        output_path=str(output_path),
-        batchsize=batchsize,
-    )
+    with pytest.warns(FutureWarning, match="split is deprecated"):
+        gfo.split(
+            input1_path=str(input1_path),
+            input2_path=str(input2_path),
+            output_path=str(output_path),
+            batchsize=-1,
+        )
 
-    # Check if the tmp file is correctly created
+    # Check if the output file is correctly created
     assert output_path.exists()
-    exp_spatial_index = GeofileInfo(output_path).default_spatial_index
-    assert gfo.has_spatial_index(output_path) is exp_spatial_index
-    input2_layerinfo = gfo.get_layerinfo(input2_path)
-    output_layerinfo = gfo.get_layerinfo(output_path)
-    assert output_layerinfo.featurecount == 68
-    assert (len(input1_layerinfo.columns) + len(input2_layerinfo.columns)) == len(
-        output_layerinfo.columns
-    )
-    assert output_layerinfo.geometrytype == GeometryType.MULTIPOLYGON
-
-    # Check the contents of the result file
-    output_gdf = gfo.read_file(output_path)
-    assert output_gdf["geometry"][0] is not None
-    input1_gdf = gfo.read_file(input1_path)
-    input2_gdf = gfo.read_file(input2_path)
-    exp_gdf = input1_gdf.overlay(input2_gdf, how="identity", keep_geom_type=True)
-    renames = dict(zip(exp_gdf.columns, output_gdf.columns))
-    exp_gdf = exp_gdf.rename(columns=renames)
-    # For text columns, gfo gives None rather than np.nan for missing values.
-    for column in exp_gdf.select_dtypes(include="O").columns:
-        exp_gdf[column] = exp_gdf[column].replace({np.nan: None})
-    # OIDN is float vs int? -> check_column_type=False
-    assert_geodataframe_equal(
-        output_gdf,
-        exp_gdf,
-        promote_to_multi=True,
-        sort_values=True,
-        check_less_precise=True,
-        normalize=True,
-        check_dtype=False,
-    )
 
 
 @pytest.mark.parametrize(