From a6555cc3ec691b84e02ed667a6b77e05a4922c62 Mon Sep 17 00:00:00 2001
From: RADHA KRISHNA KAVULURU <50633815+kssvrk@users.noreply.github.com>
Date: Fri, 14 Oct 2022 22:01:16 +0530
Subject: [PATCH] Promoting ogr_layer_algebra to official Python script (#6484)

Fixes #1581
---
 autotest/pyscripts/test_ogr_layer_algebra.py  | 499 ++++++++++++++++++
 doc/source/programs/index.rst                 |   2 +
 doc/source/programs/ogr_layer_algebra.rst     | 153 ++++++
 .../{samples => }/ogr_layer_algebra.py        |   2 +-
 .../gdal-utils/scripts/ogr_layer_algebra.py   |  12 +
 5 files changed, 667 insertions(+), 1 deletion(-)
 create mode 100644 autotest/pyscripts/test_ogr_layer_algebra.py
 create mode 100644 doc/source/programs/ogr_layer_algebra.rst
 rename swig/python/gdal-utils/osgeo_utils/{samples => }/ogr_layer_algebra.py (99%)
 create mode 100644 swig/python/gdal-utils/scripts/ogr_layer_algebra.py

diff --git a/autotest/pyscripts/test_ogr_layer_algebra.py b/autotest/pyscripts/test_ogr_layer_algebra.py
new file mode 100644
index 000000000000..da3f732057f1
--- /dev/null
+++ b/autotest/pyscripts/test_ogr_layer_algebra.py
@@ -0,0 +1,499 @@
+#!/usr/bin/env pytest
+###############################################################################
+# $Id$
+#
+# Project:  GDAL/OGR Test Suite
+# Purpose:  To test the functionality of ogr_layer_algebra script
+# Author:   Radha Krishna Kavuluru <kssvrk@gmail.com>
+#
+###############################################################################
+# Copyright (c) 2022, Radha Krishna Kavuluru <kssvrk@gmail.com>
+#
+# Permission is hereby granted, free of charge, to any person obtaining a
+# copy of this software and associated documentation files (the "Software"),
+# to deal in the Software without restriction, including without limitation
+# the rights to use, copy, modify, merge, publish, distribute, sublicense,
+# and/or sell copies of the Software, and to permit persons to whom the
+# Software is furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included
+# in all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+# DEALINGS IN THE SOFTWARE.
+###############################################################################
+
+import ogrtest
+import pytest
+import test_py_scripts
+
+from osgeo import ogr
+
+###############################################################################
+
+# Skip if no geos
+@pytest.fixture(autouse=True, scope="module")
+def startup_and_cleanup():
+    if not ogrtest.have_geos():
+        pytest.skip()
+
+
+###############################################################################
+
+# Test Intersection
+
+
+def test_ogr_layer_algebra_intersection():
+
+    script_path = test_py_scripts.get_py_script("ogr_layer_algebra")
+    if script_path is None:
+        pytest.skip()
+
+    # Create input,method,output paths for intersection.
+    input_path = "tmp/input_layer.shp"
+    method_path = "tmp/method_layer.shp"
+    output_path = "tmp/output_layer.shp"
+
+    # definition of input,method layers
+    input_layer = ogr.GetDriverByName("ESRI Shapefile").CreateDataSource(input_path)
+    method_layer = ogr.GetDriverByName("ESRI Shapefile").CreateDataSource(method_path)
+
+    A = input_layer.CreateLayer("poly")
+    B = method_layer.CreateLayer("poly")
+
+    a1 = "POLYGON((1 2, 1 3, 3 3, 3 2, 1 2))"
+    a2 = "POLYGON((5 2, 5 3, 7 3, 7 2, 5 2))"
+    b1 = "POLYGON((2 1, 2 4, 6 4, 6 1, 2 1))"
+
+    feat = ogr.Feature(A.GetLayerDefn())
+    feat.SetGeometryDirectly(ogr.Geometry(wkt=a1))
+    A.CreateFeature(feat)
+
+    feat = ogr.Feature(A.GetLayerDefn())
+    feat.SetGeometryDirectly(ogr.Geometry(wkt=a2))
+    A.CreateFeature(feat)
+
+    input_layer = None
+
+    feat = ogr.Feature(B.GetLayerDefn())
+    feat.SetGeometryDirectly(ogr.Geometry(wkt=b1))
+    B.CreateFeature(feat)
+
+    method_layer = None
+
+    # executing script
+    test_py_scripts.run_py_script(
+        script_path,
+        "ogr_layer_algebra",
+        f"Intersection -input_ds {input_path} -output_ds {output_path}  -method_ds {method_path}",
+    )
+
+    driver = ogr.GetDriverByName("ESRI Shapefile")
+    dataSource = driver.Open(output_path, 0)
+    layer = dataSource.GetLayer()
+    featureCount = layer.GetFeatureCount()
+
+    # deleting shapefiles after test
+    driver.DeleteDataSource(input_path)
+    driver.DeleteDataSource(method_path)
+    driver.DeleteDataSource(output_path)
+
+    assert featureCount == 2
+
+
+###############################################################################
+
+# Test Union
+
+
+def test_ogr_layer_algebra_union():
+
+    script_path = test_py_scripts.get_py_script("ogr_layer_algebra")
+    if script_path is None:
+        pytest.skip()
+
+    # Create input,method,output paths for intersection.
+    input_path = "tmp/input_layer.shp"
+    method_path = "tmp/method_layer.shp"
+    output_path = "tmp/output_layer.shp"
+
+    # definition of input,method layers
+    input_layer = ogr.GetDriverByName("ESRI Shapefile").CreateDataSource(input_path)
+    method_layer = ogr.GetDriverByName("ESRI Shapefile").CreateDataSource(method_path)
+
+    A = input_layer.CreateLayer("poly")
+    B = method_layer.CreateLayer("poly")
+
+    a1 = "POLYGON((1 2, 1 3, 3 3, 3 2, 1 2))"
+    a2 = "POLYGON((5 2, 5 3, 7 3, 7 2, 5 2))"
+    b1 = "POLYGON((2 1, 2 4, 6 4, 6 1, 2 1))"
+
+    feat = ogr.Feature(A.GetLayerDefn())
+    feat.SetGeometryDirectly(ogr.Geometry(wkt=a1))
+    A.CreateFeature(feat)
+
+    feat = ogr.Feature(A.GetLayerDefn())
+    feat.SetGeometryDirectly(ogr.Geometry(wkt=a2))
+    A.CreateFeature(feat)
+
+    input_layer = None
+
+    feat = ogr.Feature(B.GetLayerDefn())
+    feat.SetGeometryDirectly(ogr.Geometry(wkt=b1))
+    B.CreateFeature(feat)
+
+    method_layer = None
+
+    # executing script
+    test_py_scripts.run_py_script(
+        script_path,
+        "ogr_layer_algebra",
+        f"Union -input_ds {input_path} -output_ds {output_path}  -method_ds {method_path}",
+    )
+
+    driver = ogr.GetDriverByName("ESRI Shapefile")
+    dataSource = driver.Open(output_path, 0)
+    layer = dataSource.GetLayer()
+    featureCount = layer.GetFeatureCount()
+
+    # deleting shapefiles after test
+    driver.DeleteDataSource(input_path)
+    driver.DeleteDataSource(method_path)
+    driver.DeleteDataSource(output_path)
+
+    assert featureCount == 5
+
+
+###############################################################################
+
+# Test Symmetric Difference
+
+
+def test_ogr_layer_algebra_symdifference():
+
+    script_path = test_py_scripts.get_py_script("ogr_layer_algebra")
+    if script_path is None:
+        pytest.skip()
+
+    # Create input,method,output paths for intersection.
+    input_path = "tmp/input_layer.shp"
+    method_path = "tmp/method_layer.shp"
+    output_path = "tmp/output_layer.shp"
+
+    # definition of input,method layers
+    input_layer = ogr.GetDriverByName("ESRI Shapefile").CreateDataSource(input_path)
+    method_layer = ogr.GetDriverByName("ESRI Shapefile").CreateDataSource(method_path)
+
+    A = input_layer.CreateLayer("poly")
+    B = method_layer.CreateLayer("poly")
+
+    a1 = "POLYGON((1 2, 1 3, 3 3, 3 2, 1 2))"
+    a2 = "POLYGON((5 2, 5 3, 7 3, 7 2, 5 2))"
+    b1 = "POLYGON((2 1, 2 4, 6 4, 6 1, 2 1))"
+    b2 = "POLYGON((2 4, 2 6, 6 6, 6 4, 2 4))"
+
+    feat = ogr.Feature(A.GetLayerDefn())
+    feat.SetGeometryDirectly(ogr.Geometry(wkt=a1))
+    A.CreateFeature(feat)
+
+    feat = ogr.Feature(A.GetLayerDefn())
+    feat.SetGeometryDirectly(ogr.Geometry(wkt=a2))
+    A.CreateFeature(feat)
+
+    input_layer = None
+
+    feat = ogr.Feature(B.GetLayerDefn())
+    feat.SetGeometryDirectly(ogr.Geometry(wkt=b1))
+    B.CreateFeature(feat)
+
+    feat = ogr.Feature(B.GetLayerDefn())
+    feat.SetGeometryDirectly(ogr.Geometry(wkt=b2))
+    B.CreateFeature(feat)
+
+    method_layer = None
+
+    # executing script
+    test_py_scripts.run_py_script(
+        script_path,
+        "ogr_layer_algebra",
+        f"SymDifference -input_ds {input_path} -output_ds {output_path}  -method_ds {method_path}",
+    )
+
+    driver = ogr.GetDriverByName("ESRI Shapefile")
+    dataSource = driver.Open(output_path, 0)
+    layer = dataSource.GetLayer()
+    featureCount = layer.GetFeatureCount()
+
+    # deleting shapefiles after test
+    driver.DeleteDataSource(input_path)
+    driver.DeleteDataSource(method_path)
+    driver.DeleteDataSource(output_path)
+
+    assert featureCount == 4
+
+
+###############################################################################
+
+# Test Identity
+
+
+def test_ogr_layer_algebra_identity():
+
+    script_path = test_py_scripts.get_py_script("ogr_layer_algebra")
+    if script_path is None:
+        pytest.skip()
+
+    # Create input,method,output paths for intersection.
+    input_path = "tmp/input_layer.shp"
+    method_path = "tmp/method_layer.shp"
+    output_path = "tmp/output_layer.shp"
+
+    # definition of input,method layers
+    input_layer = ogr.GetDriverByName("ESRI Shapefile").CreateDataSource(input_path)
+    method_layer = ogr.GetDriverByName("ESRI Shapefile").CreateDataSource(method_path)
+
+    A = input_layer.CreateLayer("poly")
+    B = method_layer.CreateLayer("poly")
+
+    a1 = "POLYGON((1 2, 1 3, 3 3, 3 2, 1 2))"
+    a2 = "POLYGON((5 2, 5 3, 7 3, 7 2, 5 2))"
+    b1 = "POLYGON((2 1, 2 4, 6 4, 6 1, 2 1))"
+    b2 = "POLYGON((2 4, 2 6, 6 6, 6 4, 2 4))"
+
+    feat = ogr.Feature(A.GetLayerDefn())
+    feat.SetGeometryDirectly(ogr.Geometry(wkt=a1))
+    A.CreateFeature(feat)
+
+    feat = ogr.Feature(A.GetLayerDefn())
+    feat.SetGeometryDirectly(ogr.Geometry(wkt=a2))
+    A.CreateFeature(feat)
+
+    input_layer = None
+
+    feat = ogr.Feature(B.GetLayerDefn())
+    feat.SetGeometryDirectly(ogr.Geometry(wkt=b1))
+    B.CreateFeature(feat)
+
+    feat = ogr.Feature(B.GetLayerDefn())
+    feat.SetGeometryDirectly(ogr.Geometry(wkt=b2))
+    B.CreateFeature(feat)
+
+    method_layer = None
+
+    # executing script
+    test_py_scripts.run_py_script(
+        script_path,
+        "ogr_layer_algebra",
+        f"Identity -input_ds {input_path} -output_ds {output_path}  -method_ds {method_path}",
+    )
+
+    driver = ogr.GetDriverByName("ESRI Shapefile")
+    dataSource = driver.Open(output_path, 0)
+    layer = dataSource.GetLayer()
+    featureCount = layer.GetFeatureCount()
+
+    # deleting shapefiles after test
+    driver.DeleteDataSource(input_path)
+    driver.DeleteDataSource(method_path)
+    driver.DeleteDataSource(output_path)
+
+    assert featureCount == 4
+
+
+###############################################################################
+
+# Test Update
+
+
+def test_ogr_layer_algebra_update():
+
+    script_path = test_py_scripts.get_py_script("ogr_layer_algebra")
+    if script_path is None:
+        pytest.skip()
+
+    # Create input,method,output paths for intersection.
+    input_path = "tmp/input_layer.shp"
+    method_path = "tmp/method_layer.shp"
+    output_path = "tmp/output_layer.shp"
+
+    # definition of input,method layers
+    input_layer = ogr.GetDriverByName("ESRI Shapefile").CreateDataSource(input_path)
+    method_layer = ogr.GetDriverByName("ESRI Shapefile").CreateDataSource(method_path)
+
+    A = input_layer.CreateLayer("poly")
+    B = method_layer.CreateLayer("poly")
+
+    a1 = "POLYGON((1 2, 1 3, 3 3, 3 2, 1 2))"
+    a2 = "POLYGON((5 2, 5 3, 7 3, 7 2, 5 2))"
+    b1 = "POLYGON((2 1, 2 4, 6 4, 6 1, 2 1))"
+
+    feat = ogr.Feature(A.GetLayerDefn())
+    feat.SetGeometryDirectly(ogr.Geometry(wkt=a1))
+    A.CreateFeature(feat)
+
+    feat = ogr.Feature(A.GetLayerDefn())
+    feat.SetGeometryDirectly(ogr.Geometry(wkt=a2))
+    A.CreateFeature(feat)
+
+    input_layer = None
+
+    feat = ogr.Feature(B.GetLayerDefn())
+    feat.SetGeometryDirectly(ogr.Geometry(wkt=b1))
+    B.CreateFeature(feat)
+
+    method_layer = None
+
+    # executing script
+    test_py_scripts.run_py_script(
+        script_path,
+        "ogr_layer_algebra",
+        f"Update -input_ds {input_path} -output_ds {output_path}  -method_ds {method_path}",
+    )
+
+    driver = ogr.GetDriverByName("ESRI Shapefile")
+    dataSource = driver.Open(output_path, 0)
+    layer = dataSource.GetLayer()
+    featureCount = layer.GetFeatureCount()
+
+    # deleting shapefiles after test
+    driver.DeleteDataSource(input_path)
+    driver.DeleteDataSource(method_path)
+    driver.DeleteDataSource(output_path)
+
+    assert featureCount == 3
+
+
+###############################################################################
+
+# Test Clip
+
+
+def test_ogr_layer_algebra_clip():
+
+    script_path = test_py_scripts.get_py_script("ogr_layer_algebra")
+    if script_path is None:
+        pytest.skip()
+
+    # Create input,method,output paths for intersection.
+    input_path = "tmp/input_layer.shp"
+    method_path = "tmp/method_layer.shp"
+    output_path = "tmp/output_layer.shp"
+
+    # definition of input,method layers
+    input_layer = ogr.GetDriverByName("ESRI Shapefile").CreateDataSource(input_path)
+    method_layer = ogr.GetDriverByName("ESRI Shapefile").CreateDataSource(method_path)
+
+    A = input_layer.CreateLayer("poly")
+    B = method_layer.CreateLayer("poly")
+
+    a1 = "POLYGON((1 2, 1 3, 3 3, 3 2, 1 2))"
+    a2 = "POLYGON((5 2, 5 3, 7 3, 7 2, 5 2))"
+    b1 = "POLYGON((2 1, 2 4, 6 4, 6 1, 2 1))"
+
+    feat = ogr.Feature(A.GetLayerDefn())
+    feat.SetGeometryDirectly(ogr.Geometry(wkt=a1))
+    A.CreateFeature(feat)
+
+    feat = ogr.Feature(A.GetLayerDefn())
+    feat.SetGeometryDirectly(ogr.Geometry(wkt=a2))
+    A.CreateFeature(feat)
+
+    input_layer = None
+
+    feat = ogr.Feature(B.GetLayerDefn())
+    feat.SetGeometryDirectly(ogr.Geometry(wkt=b1))
+    B.CreateFeature(feat)
+
+    method_layer = None
+
+    # executing script
+    test_py_scripts.run_py_script(
+        script_path,
+        "ogr_layer_algebra",
+        f"Clip -input_ds {input_path} -output_ds {output_path}  -method_ds {method_path}",
+    )
+
+    driver = ogr.GetDriverByName("ESRI Shapefile")
+    dataSource = driver.Open(output_path, 0)
+    layer = dataSource.GetLayer()
+    featureCount = layer.GetFeatureCount()
+
+    # deleting shapefiles after test
+    driver.DeleteDataSource(input_path)
+    driver.DeleteDataSource(method_path)
+    driver.DeleteDataSource(output_path)
+
+    assert featureCount == 2
+
+
+###############################################################################
+
+# Test Erase
+
+
+def test_ogr_layer_algebra_erase():
+
+    script_path = test_py_scripts.get_py_script("ogr_layer_algebra")
+    if script_path is None:
+        pytest.skip()
+
+    # Create input,method,output paths for intersection.
+    input_path = "tmp/input_layer.shp"
+    method_path = "tmp/method_layer.shp"
+    output_path = "tmp/output_layer.shp"
+
+    # definition of input,method layers
+    input_layer = ogr.GetDriverByName("ESRI Shapefile").CreateDataSource(input_path)
+    method_layer = ogr.GetDriverByName("ESRI Shapefile").CreateDataSource(method_path)
+
+    A = input_layer.CreateLayer("poly")
+    B = method_layer.CreateLayer("poly")
+
+    a1 = "POLYGON((1 2, 1 3, 3 3, 3 2, 1 2))"
+    a2 = "POLYGON((5 2, 5 3, 7 3, 7 2, 5 2))"
+    b1 = "POLYGON((2 1, 2 4, 6 4, 6 1, 2 1))"
+    b2 = "POLYGON((2 4, 2 6, 6 6, 6 4, 2 4))"
+
+    feat = ogr.Feature(A.GetLayerDefn())
+    feat.SetGeometryDirectly(ogr.Geometry(wkt=a1))
+    A.CreateFeature(feat)
+
+    feat = ogr.Feature(A.GetLayerDefn())
+    feat.SetGeometryDirectly(ogr.Geometry(wkt=a2))
+    A.CreateFeature(feat)
+
+    input_layer = None
+
+    feat = ogr.Feature(B.GetLayerDefn())
+    feat.SetGeometryDirectly(ogr.Geometry(wkt=b1))
+    B.CreateFeature(feat)
+
+    feat = ogr.Feature(B.GetLayerDefn())
+    feat.SetGeometryDirectly(ogr.Geometry(wkt=b2))
+    B.CreateFeature(feat)
+
+    method_layer = None
+
+    # executing script
+    test_py_scripts.run_py_script(
+        script_path,
+        "ogr_layer_algebra",
+        f"Erase -input_ds {input_path} -output_ds {output_path}  -method_ds {method_path}",
+    )
+
+    driver = ogr.GetDriverByName("ESRI Shapefile")
+    dataSource = driver.Open(output_path, 0)
+    layer = dataSource.GetLayer()
+    featureCount = layer.GetFeatureCount()
+
+    # deleting shapefiles after test
+    driver.DeleteDataSource(input_path)
+    driver.DeleteDataSource(method_path)
+    driver.DeleteDataSource(output_path)
+
+    assert featureCount == 2
diff --git a/doc/source/programs/index.rst b/doc/source/programs/index.rst
index dc00883c3d46..31403002bf89 100644
--- a/doc/source/programs/index.rst
+++ b/doc/source/programs/index.rst
@@ -113,6 +113,7 @@ Vector programs
    ogrtindex
    ogrlineref
    ogrmerge
+   ogr_layer_algebra
 
 .. only:: html
 
@@ -122,6 +123,7 @@ Vector programs
     - :ref:`ogrtindex`: Creates a tileindex.
     - :ref:`ogrlineref`: Create linear reference and provide some calculations using it.
     - :ref:`ogrmerge`: Merge several vector datasets into a single one.
+    - :ref:`ogr_layer_algebra`: Performs various Vector layer algebraic operations.
 
 Geographic network programs
 ---------------------------
diff --git a/doc/source/programs/ogr_layer_algebra.rst b/doc/source/programs/ogr_layer_algebra.rst
new file mode 100644
index 000000000000..92b0977e73b4
--- /dev/null
+++ b/doc/source/programs/ogr_layer_algebra.rst
@@ -0,0 +1,153 @@
+.. _ogr_layer_algebra:
+
+================================================================================
+ogr_layer_algebra.py
+================================================================================
+
+.. versionadded:: 3.6
+
+.. only:: html
+
+    Performs various Vector layer algebraic operations.
+
+.. Index:: ogr_layer_algebra
+
+Synopsis
+--------
+
+.. code-block::
+
+    ogr_layer_algebra.py Union|Intersection|SymDifference|Identity|Update|Clip|Erase
+                        -input_ds name [-input_lyr name]
+                        -method_ds [-method_lyr name]
+                        -output_ds name [-output_lyr name] [-overwrite]
+                        [-opt NAME=VALUE]*
+                        [-f format_name] [-dsco NAME=VALUE]* [-lco NAME=VALUE]*
+                        [-input_fields NONE|ALL|fld1,fl2,...fldN] [-method_fields NONE|ALL|fld1,fl2,...fldN]
+                        [-nlt geom_type] [-a_srs srs_def]
+
+Description
+-----------
+
+The :program:`ogr_layer_algebra.py` provides a command line utility to perform various vector layer algebraic operations. The utility takes a vector 
+input source , a method source and generates the output of the operation in the specified output file
+
+.. program:: ogr_layer_algebra
+
+.. option:: <mode>
+
+    Where <mode> is one of the seven available modes:
+
+    * ``Union``
+
+        A union is a set of features, which represent areas that are in either of the operand layers. 
+
+    * ``Intersection``
+
+        An intersection is a set of features, which represent the common areas of two layers. 
+
+    * ``SymDifference``
+
+        A symmetric difference is a set of features, which represent areas that are in operand layers but which do not intersect. 
+
+    * ``Identity``
+
+        The identity method identifies features in the input layer with features in the method layer possibly splitting features into several features. 
+        By default the result layer has attributes from both operand layers. 
+
+    * ``Update``
+
+        The update method creates a layer, which add features into the input layer from the method layer possibly cutting features in the input layer. 
+        By default the result layer has attributes only from the input layer. 
+
+    * ``Clip``
+
+        The clip method creates a layer, which has features from the input layer clipped to the areas of the features in the method layer. 
+        By default the result layer has attributes of the input layer. 
+
+    * ``Erase``
+
+        The erase method creates a layer, which has features from the input layer whose areas are erased by the features in the method layer. 
+        By default the result layer has attributes of the input layer. 
+        
+.. option:: -input_ds <path>
+
+    Input data set path for the operation to be performed.
+    For operations involving 2 datasets, this is one of the dataset.
+
+.. option:: -input_lyr <name>
+
+    Layer name of the ``input_ds`` for which the operations have to be performed ( Optional )
+
+.. option:: -method_ds <path>
+
+    Method data set path for the operation to be performed. 
+    This is usually the conditional data set supplied to the operation ( ex: clip , erase , update )
+    This is the Second data set in the operation ( ex : Union, Intersection , SymDifference )
+
+.. option:: -method_lyr <name>
+
+    Layer name of the ``method_ds`` for which the operations have to be performed ( Optional )
+
+.. option:: -output_ds <path>
+
+    Output data set path for writing the result of the operations performed by ogr_layer_algebra
+
+.. option:: -output_lyr_name <name>
+
+    Layer name of the ``output_lyr_name`` where the output vector has to be written. ( Optional )
+
+.. option:: -overwrite
+
+    Indicates wether the ``output_ds`` have to be overwritten with the generated result of ogr_layer_algebra
+
+.. option:: -opt <NAME=VALUE>
+
+    Attributes for which the operation has to run on ``input_ds`` and ``method_ds``
+
+.. option:: -f <format_name>
+
+    Select the output format.If not specified,
+    the format is guessed from the extension (previously was ESRI Shapefile).
+    Use the short format name
+
+.. option:: -dsco <NAME=VALUE>
+    
+    Dataset creation option (format specific)
+
+.. option:: -lco <NAME=VALUE>
+    
+    Layer creation option (format specific)
+
+.. option:: -input_fields <NONE|ALL|fld1,fld2,fld3...>
+
+    Comma-delimited list of fields from input layer to copy to the output layer ,
+    if eligible according to the operation
+
+.. option:: -method_fields <NONE|ALL|fld1,fld2,fld3...>
+
+    Comma-delimited list of fields from method layer to copy to the output layer ,
+    if eligible according to the operation
+
+.. option:: -nlt <geom_type>
+
+    Define the geometry type for the created layer. 
+    One of NONE, GEOMETRY, POINT, LINESTRING, POLYGON, GEOMETRYCOLLECTION,
+    MULTIPOINT, MULTIPOLYGON, GEOMETRY25D, POINT25D, LINESTRING25D, POLYGON25D,
+    GEOMETRYCOLLECTION25D, MULTIPOINT25D, MULTIPOLYGON25D.
+
+.. option:: -a_srs <srs_def>
+
+    Assign an output SRS, but without reprojecting
+
+    The coordinate reference systems that can be passed are anything supported by the
+    OGRSpatialReference.SetFromUserInput() call, which includes EPSG Projected,
+    Geographic or Compound CRS (i.e. EPSG:4296), a well known text (WKT) CRS definition, 
+    PROJ.4 declarations, or the name of a .prj file containing a WKT CRS definition.
+
+
+
+
+
+
+
diff --git a/swig/python/gdal-utils/osgeo_utils/samples/ogr_layer_algebra.py b/swig/python/gdal-utils/osgeo_utils/ogr_layer_algebra.py
similarity index 99%
rename from swig/python/gdal-utils/osgeo_utils/samples/ogr_layer_algebra.py
rename to swig/python/gdal-utils/osgeo_utils/ogr_layer_algebra.py
index 221cc811a7f0..0edd1d8887a4 100644
--- a/swig/python/gdal-utils/osgeo_utils/samples/ogr_layer_algebra.py
+++ b/swig/python/gdal-utils/osgeo_utils/ogr_layer_algebra.py
@@ -158,7 +158,7 @@ def CreateLayer(
 ###############################################################################
 
 
-def main(argv=None):
+def main(argv=sys.argv):
     version_num = int(gdal.VersionInfo("VERSION_NUM"))
     if version_num < 1100000:
         print("ERROR: Python bindings of GDAL 1.10 or later required")
diff --git a/swig/python/gdal-utils/scripts/ogr_layer_algebra.py b/swig/python/gdal-utils/scripts/ogr_layer_algebra.py
new file mode 100644
index 000000000000..6acde54ef4b2
--- /dev/null
+++ b/swig/python/gdal-utils/scripts/ogr_layer_algebra.py
@@ -0,0 +1,12 @@
+#!/usr/bin/env python3
+
+import sys
+
+from osgeo.gdal import deprecation_warn
+
+# import osgeo_utils.ogr_layer_algebra as a convenience to use as a script
+from osgeo_utils.ogr_layer_algebra import *  # noqa
+from osgeo_utils.ogr_layer_algebra import main
+
+deprecation_warn("ogr_layer_algebra")
+sys.exit(main(sys.argv))