From 48b65fedc627ae9ecef15851d7f72ec84b8a9b3f Mon Sep 17 00:00:00 2001
From: Even Rouault <even.rouault@spatialys.com>
Date: Fri, 17 Jan 2025 01:44:16 +0100
Subject: [PATCH] Add 'gdal vector clip' and 'gdal vector pipeline read ... !
 clip ... ! write ...'

---
 apps/CMakeLists.txt                           |   1 +
 apps/gdalalg_vector.cpp                       |   2 +
 apps/gdalalg_vector_clip.cpp                  | 389 ++++++++++++
 apps/gdalalg_vector_clip.h                    |  64 ++
 apps/gdalalg_vector_pipeline.cpp              |   2 +
 .../utilities/test_gdalalg_vector_clip.py     | 575 ++++++++++++++++++
 doc/source/conf.py                            |   7 +
 .../programs/gdal_options/co_vector.rst       |  18 +
 doc/source/programs/gdal_vector.rst           |   4 +-
 doc/source/programs/gdal_vector_clip.rst      | 132 ++++
 doc/source/programs/gdal_vector_convert.rst   |   2 +-
 doc/source/programs/gdal_vector_pipeline.rst  |  18 +
 doc/source/programs/index.rst                 |   2 +
 13 files changed, 1214 insertions(+), 2 deletions(-)
 create mode 100644 apps/gdalalg_vector_clip.cpp
 create mode 100644 apps/gdalalg_vector_clip.h
 create mode 100755 autotest/utilities/test_gdalalg_vector_clip.py
 create mode 100644 doc/source/programs/gdal_options/co_vector.rst
 create mode 100644 doc/source/programs/gdal_vector_clip.rst

diff --git a/apps/CMakeLists.txt b/apps/CMakeLists.txt
index a080179aa656..a89f8c456f15 100644
--- a/apps/CMakeLists.txt
+++ b/apps/CMakeLists.txt
@@ -24,6 +24,7 @@ add_library(
   gdalalg_raster_write.cpp
   gdalalg_vector.cpp
   gdalalg_vector_info.cpp
+  gdalalg_vector_clip.cpp
   gdalalg_vector_convert.cpp
   gdalalg_vector_pipeline.cpp
   gdalalg_vector_read.cpp
diff --git a/apps/gdalalg_vector.cpp b/apps/gdalalg_vector.cpp
index 41a15da845ef..75191f2c01ba 100644
--- a/apps/gdalalg_vector.cpp
+++ b/apps/gdalalg_vector.cpp
@@ -13,6 +13,7 @@
 #include "gdalalgorithm.h"
 
 #include "gdalalg_vector_info.h"
+#include "gdalalg_vector_clip.h"
 #include "gdalalg_vector_convert.h"
 #include "gdalalg_vector_pipeline.h"
 #include "gdalalg_vector_filter.h"
@@ -37,6 +38,7 @@ class GDALVectorAlgorithm final : public GDALAlgorithm
     GDALVectorAlgorithm() : GDALAlgorithm(NAME, DESCRIPTION, HELP_URL)
     {
         RegisterSubAlgorithm<GDALVectorInfoAlgorithm>();
+        RegisterSubAlgorithm<GDALVectorClipAlgorithmStandalone>();
         RegisterSubAlgorithm<GDALVectorConvertAlgorithm>();
         RegisterSubAlgorithm<GDALVectorPipelineAlgorithm>();
         RegisterSubAlgorithm<GDALVectorFilterAlgorithmStandalone>();
diff --git a/apps/gdalalg_vector_clip.cpp b/apps/gdalalg_vector_clip.cpp
new file mode 100644
index 000000000000..2d5e15f4af9d
--- /dev/null
+++ b/apps/gdalalg_vector_clip.cpp
@@ -0,0 +1,389 @@
+/******************************************************************************
+ *
+ * Project:  GDAL
+ * Purpose:  "clip" step of "vector pipeline", or "gdal vector clip" standalone
+ * Author:   Even Rouault <even dot rouault at spatialys.com>
+ *
+ ******************************************************************************
+ * Copyright (c) 2025, Even Rouault <even dot rouault at spatialys.com>
+ *
+ * SPDX-License-Identifier: MIT
+ ****************************************************************************/
+
+#include "gdalalg_vector_clip.h"
+
+#include "gdal_priv.h"
+#include "gdal_utils.h"
+#include "ogrsf_frmts.h"
+
+#include <algorithm>
+
+//! @cond Doxygen_Suppress
+
+#ifndef _
+#define _(x) (x)
+#endif
+
+/************************************************************************/
+/*          GDALVectorClipAlgorithm::GDALVectorClipAlgorithm()          */
+/************************************************************************/
+
+GDALVectorClipAlgorithm::GDALVectorClipAlgorithm(bool standaloneStep)
+    : GDALVectorPipelineStepAlgorithm(NAME, DESCRIPTION, HELP_URL,
+                                      standaloneStep)
+{
+    AddBBOXArg(&m_bbox, _("Clipping bounding box as xmin,ymin,xmax,ymax"))
+        .SetMutualExclusionGroup("exclusion-group");
+    AddArg("bbox-crs", 0, _("CRS of clipping bounding box"), &m_bboxCrs)
+        .SetIsCRSArg()
+        .AddHiddenAlias("bbox_srs");
+    AddArg("geometry", 0, _("Clipping geometry (WKT or GeoJSON)"), &m_geometry)
+        .SetMutualExclusionGroup("exclusion-group");
+    AddArg("geometry-crs", 0, _("CRS of clipping geometry"), &m_geometryCrs)
+        .SetIsCRSArg()
+        .AddHiddenAlias("geometry_srs");
+    AddArg("like", 0, _("Dataset to use as a template for bounds"),
+           &m_likeDataset, GDAL_OF_RASTER | GDAL_OF_VECTOR)
+        .SetMetaVar("DATASET")
+        .SetMutualExclusionGroup("exclusion-group");
+}
+
+/************************************************************************/
+/*                   GDALVectorClipAlgorithmDataset                     */
+/************************************************************************/
+
+namespace
+{
+class GDALVectorClipAlgorithmDataset final : public GDALDataset
+{
+    std::vector<std::unique_ptr<OGRLayer>> m_layers{};
+
+  public:
+    GDALVectorClipAlgorithmDataset() = default;
+
+    void AddLayer(std::unique_ptr<OGRLayer> poLayer)
+    {
+        m_layers.push_back(std::move(poLayer));
+    }
+
+    int GetLayerCount() override
+    {
+        return static_cast<int>(m_layers.size());
+    }
+
+    OGRLayer *GetLayer(int idx) override
+    {
+        return idx >= 0 && idx < GetLayerCount() ? m_layers[idx].get()
+                                                 : nullptr;
+    }
+};
+
+class GDALVectorClipAlgorithmLayer final : public OGRLayer
+{
+  public:
+    GDALVectorClipAlgorithmLayer(OGRLayer *poSrcLayer,
+                                 std::unique_ptr<OGRGeometry> poClipGeom)
+        : m_poSrcLayer(poSrcLayer), m_poClipGeom(std::move(poClipGeom)),
+          m_eSrcLayerGeomType(wkbFlatten(m_poSrcLayer->GetGeomType())),
+          m_bSrcLayerGeomTypeIsCollection(
+              OGR_GT_IsSubClassOf(m_eSrcLayerGeomType, wkbGeometryCollection))
+    {
+        SetDescription(poSrcLayer->GetDescription());
+        poSrcLayer->SetSpatialFilter(m_poClipGeom.get());
+    }
+
+    OGRFeatureDefn *GetLayerDefn() override
+    {
+        return m_poSrcLayer->GetLayerDefn();
+    }
+
+    void ResetReading() override
+    {
+        m_poSrcLayer->ResetReading();
+        m_poSrcFeature.reset();
+        m_poCurGeomColl.reset();
+        m_idxInCurGeomColl = 0;
+    }
+
+    OGRFeature *GetNextFeature() override
+    {
+        if (m_poSrcFeature && m_poCurGeomColl)
+        {
+            if (m_idxInCurGeomColl < m_poCurGeomColl->getNumGeometries())
+            {
+                auto poDstFeature =
+                    std::unique_ptr<OGRFeature>(m_poSrcFeature->Clone());
+                poDstFeature->SetGeometryDirectly(
+                    m_poCurGeomColl->getGeometryRef(m_idxInCurGeomColl)
+                        ->clone());
+                ++m_idxInCurGeomColl;
+                return poDstFeature.release();
+            }
+            m_poSrcFeature.reset();
+            m_poCurGeomColl.reset();
+            m_idxInCurGeomColl = 0;
+        }
+
+        while (auto poFeature =
+                   std::unique_ptr<OGRFeature>(m_poSrcLayer->GetNextFeature()))
+        {
+            auto poGeom = poFeature->GetGeometryRef();
+            if (!poGeom)
+                continue;
+
+            auto poIntersection = std::unique_ptr<OGRGeometry>(
+                poGeom->Intersection(m_poClipGeom.get()));
+            if (!poIntersection)
+                continue;
+
+            const auto eFeatGeomType =
+                wkbFlatten(poIntersection->getGeometryType());
+            if (m_eSrcLayerGeomType != wkbUnknown &&
+                m_eSrcLayerGeomType != eFeatGeomType &&
+                !m_bSrcLayerGeomTypeIsCollection &&
+                OGR_GT_IsSubClassOf(eFeatGeomType, wkbGeometryCollection))
+            {
+                m_poSrcFeature = std::move(poFeature);
+                m_poCurGeomColl.reset(
+                    poIntersection.release()->toGeometryCollection());
+                m_idxInCurGeomColl = 0;
+                return GetNextFeature();
+            }
+            else
+            {
+                poFeature->SetGeometryDirectly(poIntersection.release());
+                return poFeature.release();
+            }
+        }
+        return nullptr;
+    }
+
+    int TestCapability(const char *pszCap) override
+    {
+        if (EQUAL(pszCap, OLCStringsAsUTF8) ||
+            EQUAL(pszCap, OLCCurveGeometries) || EQUAL(pszCap, OLCZGeometries))
+            return m_poSrcLayer->TestCapability(pszCap);
+        return false;
+    }
+
+  private:
+    OGRLayer *m_poSrcLayer = nullptr;
+    std::unique_ptr<OGRGeometry> m_poClipGeom{};
+    const OGRwkbGeometryType m_eSrcLayerGeomType;
+    const bool m_bSrcLayerGeomTypeIsCollection;
+    std::unique_ptr<OGRFeature> m_poSrcFeature{};
+    std::unique_ptr<OGRGeometryCollection> m_poCurGeomColl{};
+    int m_idxInCurGeomColl = 0;
+
+    CPL_DISALLOW_COPY_ASSIGN(GDALVectorClipAlgorithmLayer)
+};
+
+}  // namespace
+
+/************************************************************************/
+/*                 GDALVectorClipAlgorithm::RunStep()                   */
+/************************************************************************/
+
+bool GDALVectorClipAlgorithm::RunStep(GDALProgressFunc, void *)
+{
+    CPLAssert(m_inputDataset.GetDatasetRef());
+    CPLAssert(m_outputDataset.GetName().empty());
+    CPLAssert(!m_outputDataset.GetDatasetRef());
+
+    auto poSrcDS = m_inputDataset.GetDatasetRef();
+
+    std::unique_ptr<OGRGeometry> poClipGeom;
+
+    const int nLayerCount = poSrcDS->GetLayerCount();
+    bool bSrcLayerHasSRS = false;
+    for (int i = 0; i < nLayerCount; ++i)
+    {
+        auto poSrcLayer = poSrcDS->GetLayer(i);
+        if (poSrcLayer && poSrcLayer->GetSpatialRef())
+        {
+            bSrcLayerHasSRS = true;
+            break;
+        }
+    }
+
+    if (!m_bbox.empty())
+    {
+        poClipGeom = std::make_unique<OGRPolygon>(m_bbox[0], m_bbox[1],
+                                                  m_bbox[2], m_bbox[3]);
+
+        if (!m_bboxCrs.empty())
+        {
+            auto poSRS = new OGRSpatialReference();
+            poSRS->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
+            CPL_IGNORE_RET_VAL(poSRS->SetFromUserInput(m_bboxCrs.c_str()));
+            poClipGeom->assignSpatialReference(poSRS);
+            poSRS->Release();
+        }
+    }
+    else if (!m_geometry.empty())
+    {
+        {
+            CPLErrorStateBackuper oErrorStateBackuper(CPLQuietErrorHandler);
+            auto [poGeom, eErr] =
+                OGRGeometryFactory::createFromWkt(m_geometry.c_str());
+            if (eErr == OGRERR_NONE)
+            {
+                poClipGeom = std::move(poGeom);
+            }
+            else
+            {
+                poClipGeom.reset(
+                    OGRGeometryFactory::createFromGeoJson(m_geometry.c_str()));
+                if (poClipGeom)
+                {
+                    auto poSRS = new OGRSpatialReference();
+                    poSRS->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
+                    CPL_IGNORE_RET_VAL(poSRS->SetFromUserInput("WGS84"));
+                    poClipGeom->assignSpatialReference(poSRS);
+                    poSRS->Release();
+                }
+            }
+        }
+        if (!poClipGeom)
+        {
+            ReportError(
+                CE_Failure, CPLE_AppDefined,
+                "Clipping geometry is neither a valid WKT or GeoJSON geometry");
+            return false;
+        }
+
+        if (!m_geometryCrs.empty())
+        {
+            auto poSRS = new OGRSpatialReference();
+            poSRS->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
+            CPL_IGNORE_RET_VAL(poSRS->SetFromUserInput(m_geometryCrs.c_str()));
+            poClipGeom->assignSpatialReference(poSRS);
+            poSRS->Release();
+        }
+    }
+    else if (auto poLikeDS = m_likeDataset.GetDatasetRef())
+    {
+        if (poLikeDS->GetLayerCount() > 1)
+        {
+            ReportError(
+                CE_Failure, CPLE_AppDefined,
+                "Only single layer dataset can be specified with --like");
+            return false;
+        }
+        else if (poLikeDS->GetLayerCount() == 1)
+        {
+            auto poLikeLayer = poLikeDS->GetLayer(0);
+            OGREnvelope sEnvelope;
+            if (poLikeLayer->GetExtent(&sEnvelope, true) != OGRERR_NONE)
+            {
+                ReportError(CE_Failure, CPLE_AppDefined,
+                            "Cannot get extent of clip layer");
+                return false;
+            }
+            poClipGeom = std::make_unique<OGRPolygon>(
+                sEnvelope.MinX, sEnvelope.MinY, sEnvelope.MaxX, sEnvelope.MaxY);
+            auto poLikeSRS = poLikeLayer->GetSpatialRef();
+            if (bSrcLayerHasSRS && !poLikeSRS)
+            {
+                ReportError(CE_Warning, CPLE_AppDefined,
+                            "Dataset '%s' has no SRS. Assuming its SRS is the "
+                            "same as the input vector.",
+                            poLikeDS->GetDescription());
+            }
+            poClipGeom->assignSpatialReference(poLikeSRS);
+        }
+        else if (poLikeDS->GetRasterCount() > 0)
+        {
+            double adfGT[6];
+            if (poLikeDS->GetGeoTransform(adfGT) != CE_None)
+            {
+                ReportError(
+                    CE_Failure, CPLE_AppDefined,
+                    "Dataset '%s' has no geotransform matrix. Its bounds "
+                    "cannot be established.",
+                    poLikeDS->GetDescription());
+                return false;
+            }
+            auto poLikeSRS = poLikeDS->GetSpatialRef();
+            if (bSrcLayerHasSRS && !poLikeSRS)
+            {
+                ReportError(CE_Warning, CPLE_AppDefined,
+                            "Dataset '%s' has no SRS. Assuming its SRS is the "
+                            "same as the input vector.",
+                            poLikeDS->GetDescription());
+            }
+            const double dfTLX = adfGT[0];
+            const double dfTLY = adfGT[3];
+            const double dfTRX =
+                adfGT[0] + poLikeDS->GetRasterXSize() * adfGT[1];
+            const double dfTRY =
+                adfGT[3] + poLikeDS->GetRasterXSize() * adfGT[4];
+            const double dfBLX =
+                adfGT[0] + poLikeDS->GetRasterYSize() * adfGT[2];
+            const double dfBLY =
+                adfGT[3] + poLikeDS->GetRasterYSize() * adfGT[5];
+            const double dfBRX = adfGT[0] +
+                                 poLikeDS->GetRasterXSize() * adfGT[1] +
+                                 poLikeDS->GetRasterYSize() * adfGT[2];
+            const double dfBRY = adfGT[3] +
+                                 poLikeDS->GetRasterXSize() * adfGT[4] +
+                                 poLikeDS->GetRasterYSize() * adfGT[5];
+
+            auto poPoly = std::make_unique<OGRPolygon>();
+            auto poLR = std::make_unique<OGRLinearRing>();
+            poLR->addPoint(dfTLX, dfTLY);
+            poLR->addPoint(dfTRX, dfTRY);
+            poLR->addPoint(dfBRX, dfBRY);
+            poLR->addPoint(dfBLX, dfBLY);
+            poLR->addPoint(dfTLX, dfTLY);
+            poPoly->addRingDirectly(poLR.release());
+            poPoly->assignSpatialReference(poLikeSRS);
+            poClipGeom = std::move(poPoly);
+        }
+        else
+        {
+            ReportError(CE_Failure, CPLE_AppDefined,
+                        "Cannot get extent from clip dataset");
+            return false;
+        }
+    }
+    else
+    {
+        ReportError(CE_Failure, CPLE_AppDefined,
+                    "--bbox, --geometry or --like must be specified");
+        return false;
+    }
+
+    auto outDS = std::make_unique<GDALVectorClipAlgorithmDataset>();
+    outDS->SetDescription(poSrcDS->GetDescription());
+
+    bool ret = true;
+    for (int i = 0; ret && i < nLayerCount; ++i)
+    {
+        auto poSrcLayer = poSrcDS->GetLayer(i);
+        ret = (poSrcLayer != nullptr);
+        if (ret)
+        {
+            auto poClipGeomForLayer =
+                std::unique_ptr<OGRGeometry>(poClipGeom->clone());
+            if (poClipGeomForLayer->getSpatialReference() &&
+                poSrcLayer->GetSpatialRef())
+            {
+                ret = poClipGeomForLayer->transformTo(
+                          poSrcLayer->GetSpatialRef()) == OGRERR_NONE;
+            }
+            if (ret)
+            {
+                outDS->AddLayer(std::make_unique<GDALVectorClipAlgorithmLayer>(
+                    poSrcLayer, std::move(poClipGeomForLayer)));
+            }
+        }
+    }
+
+    if (ret)
+        m_outputDataset.Set(std::move(outDS));
+
+    return ret;
+}
+
+//! @endcond
diff --git a/apps/gdalalg_vector_clip.h b/apps/gdalalg_vector_clip.h
new file mode 100644
index 000000000000..8eb75ca6b3d8
--- /dev/null
+++ b/apps/gdalalg_vector_clip.h
@@ -0,0 +1,64 @@
+/******************************************************************************
+ *
+ * Project:  GDAL
+ * Purpose:  "clip" step of "vector pipeline", or "gdal vector clip" standalone
+ * Author:   Even Rouault <even dot rouault at spatialys.com>
+ *
+ ******************************************************************************
+ * Copyright (c) 2025, Even Rouault <even dot rouault at spatialys.com>
+ *
+ * SPDX-License-Identifier: MIT
+ ****************************************************************************/
+
+#ifndef GDALALG_VECTOR_CLIP_INCLUDED
+#define GDALALG_VECTOR_CLIP_INCLUDED
+
+#include "gdalalg_vector_pipeline.h"
+
+//! @cond Doxygen_Suppress
+
+/************************************************************************/
+/*                       GDALVectorClipAlgorithm                        */
+/************************************************************************/
+
+class GDALVectorClipAlgorithm /* non final */
+    : public GDALVectorPipelineStepAlgorithm
+{
+  public:
+    static constexpr const char *NAME = "clip";
+    static constexpr const char *DESCRIPTION = "Clip a vector dataset.";
+    static constexpr const char *HELP_URL = "/programs/gdal_vector_clip.html";
+
+    static std::vector<std::string> GetAliases()
+    {
+        return {};
+    }
+
+    explicit GDALVectorClipAlgorithm(bool standaloneStep = false);
+
+  private:
+    bool RunStep(GDALProgressFunc pfnProgress, void *pProgressData) override;
+
+    std::vector<double> m_bbox{};
+    std::string m_bboxCrs{};
+    std::string m_geometry{};
+    std::string m_geometryCrs{};
+    GDALArgDatasetValue m_likeDataset{};
+};
+
+/************************************************************************/
+/*                   GDALVectorClipAlgorithmStandalone                  */
+/************************************************************************/
+
+class GDALVectorClipAlgorithmStandalone final : public GDALVectorClipAlgorithm
+{
+  public:
+    GDALVectorClipAlgorithmStandalone()
+        : GDALVectorClipAlgorithm(/* standaloneStep = */ true)
+    {
+    }
+};
+
+//! @endcond
+
+#endif /* GDALALG_VECTOR_CLIP_INCLUDED */
diff --git a/apps/gdalalg_vector_pipeline.cpp b/apps/gdalalg_vector_pipeline.cpp
index 7bee6af9a19a..eba4c0ae10ab 100644
--- a/apps/gdalalg_vector_pipeline.cpp
+++ b/apps/gdalalg_vector_pipeline.cpp
@@ -12,6 +12,7 @@
 
 #include "gdalalg_vector_pipeline.h"
 #include "gdalalg_vector_read.h"
+#include "gdalalg_vector_clip.h"
 #include "gdalalg_vector_filter.h"
 #include "gdalalg_vector_reproject.h"
 #include "gdalalg_vector_write.h"
@@ -174,6 +175,7 @@ GDALVectorPipelineAlgorithm::GDALVectorPipelineAlgorithm()
 
     m_stepRegistry.Register<GDALVectorReadAlgorithm>();
     m_stepRegistry.Register<GDALVectorWriteAlgorithm>();
+    m_stepRegistry.Register<GDALVectorClipAlgorithm>();
     m_stepRegistry.Register<GDALVectorReprojectAlgorithm>();
     m_stepRegistry.Register<GDALVectorFilterAlgorithm>();
 }
diff --git a/autotest/utilities/test_gdalalg_vector_clip.py b/autotest/utilities/test_gdalalg_vector_clip.py
new file mode 100755
index 000000000000..c8d4fa15c62a
--- /dev/null
+++ b/autotest/utilities/test_gdalalg_vector_clip.py
@@ -0,0 +1,575 @@
+#!/usr/bin/env pytest
+# -*- coding: utf-8 -*-
+###############################################################################
+# Project:  GDAL/OGR Test Suite
+# Purpose:  'gdal vector clip' testing
+# Author:   Even Rouault <even dot rouault @ spatialys.com>
+#
+###############################################################################
+# Copyright (c) 2025, Even Rouault <even dot rouault at spatialys.com>
+#
+# SPDX-License-Identifier: MIT
+###############################################################################
+
+import ogrtest
+import pytest
+
+from osgeo import gdal, ogr, osr
+
+pytestmark = pytest.mark.require_geos
+
+
+def get_clip_alg():
+    reg = gdal.GetGlobalAlgorithmRegistry()
+    vector = reg.InstantiateAlg("vector")
+    return vector.InstantiateSubAlgorithm("clip")
+
+
+@pytest.mark.require_driver("GPKG")
+def test_gdalalg_vector_clip_general_behavior(tmp_vsimem):
+
+    out_filename = str(tmp_vsimem / "out.gpkg")
+
+    clip = get_clip_alg()
+    assert clip.ParseRunAndFinalize(
+        ["--bbox", "-1e10,-1e10,1e10,1e10", "../ogr/data/poly.shp", out_filename]
+    )
+
+    with gdal.OpenEx(out_filename, gdal.OF_UPDATE) as ds:
+        assert ds.GetLayer(0).GetFeatureCount() == 10
+        for i in range(10):
+            ds.GetLayer(0).DeleteFeature(i + 1)
+
+    clip = get_clip_alg()
+    with pytest.raises(
+        Exception, match="already exists. Specify the --overwrite option"
+    ):
+        clip.ParseRunAndFinalize(
+            ["--bbox", "-1e10,-1e10,1e10,1e10", "../ogr/data/poly.shp", out_filename]
+        )
+
+    with gdal.OpenEx(out_filename) as ds:
+        assert ds.GetLayer(0).GetFeatureCount() == 0
+
+    clip = get_clip_alg()
+    assert clip.ParseRunAndFinalize(
+        [
+            "--bbox",
+            "-1e10,-1e10,1e10,1e10",
+            "--overwrite",
+            "../ogr/data/poly.shp",
+            out_filename,
+        ]
+    )
+
+    with gdal.OpenEx(out_filename) as ds:
+        assert ds.GetLayer(0).GetFeatureCount() == 10
+
+    clip = get_clip_alg()
+    assert clip.ParseRunAndFinalize(
+        [
+            "--bbox",
+            "-1e10,-1e10,1e10,1e10",
+            "--append",
+            "../ogr/data/poly.shp",
+            out_filename,
+        ]
+    )
+
+    with gdal.OpenEx(out_filename) as ds:
+        assert ds.GetLayer(0).GetFeatureCount() == 20
+
+    clip = get_clip_alg()
+    assert clip.ParseRunAndFinalize(
+        [
+            "--bbox",
+            "-1e10,-1e10,1e10,1e10",
+            "--update",
+            "--nln",
+            "layer2",
+            "../ogr/data/poly.shp",
+            out_filename,
+        ]
+    )
+
+    with gdal.OpenEx(out_filename) as ds:
+        assert ds.GetLayerByName("poly").GetFeatureCount() == 20
+        assert ds.GetLayerByName("layer2").GetFeatureCount() == 10
+
+    clip = get_clip_alg()
+    assert clip.ParseRunAndFinalize(
+        [
+            "--bbox",
+            "-1e10,-1e10,1e10,1e10",
+            "--of",
+            "GPKG",
+            "--overwrite-layer",
+            "--nln",
+            "poly",
+            "../ogr/data/poly.shp",
+            out_filename,
+        ]
+    )
+
+    with gdal.OpenEx(out_filename) as ds:
+        assert ds.GetLayerByName("poly").GetFeatureCount() == 10
+        assert ds.GetLayerByName("layer2").GetFeatureCount() == 10
+
+
+def test_gdalalg_vector_clip_bbox():
+
+    src_ds = gdal.GetDriverByName("Memory").Create("", 0, 0, 0, gdal.GDT_Unknown)
+    src_lyr = src_ds.CreateLayer("test")
+    src_lyr.CreateField(ogr.FieldDefn("foo"))
+
+    f = ogr.Feature(src_lyr.GetLayerDefn())
+    f["foo"] = "bar"
+    f.SetGeometry(ogr.CreateGeometryFromWkt("POLYGON ((0 0,0 -1,-1 -1,-1 0,0 0))"))
+    src_lyr.CreateFeature(f)
+
+    f = ogr.Feature(src_lyr.GetLayerDefn())
+    f["foo"] = "bar"
+    f.SetGeometry(ogr.CreateGeometryFromWkt("POLYGON ((0 0,0 1,1 1,1 0,0 0))"))
+    src_lyr.CreateFeature(f)
+
+    clip = get_clip_alg()
+    clip.GetArg("input").Get().SetDataset(src_ds)
+
+    assert clip.ParseCommandLineArguments(
+        ["--bbox", "0.2,0.3,0.7,0.8", "--of", "Memory", "--output", "memory_ds"]
+    )
+    assert clip.Run()
+
+    out_ds = clip.GetArg("output").Get().GetDataset()
+    out_lyr = out_ds.GetLayer(0)
+    out_f = out_lyr.GetNextFeature()
+    assert out_f["foo"] == "bar"
+    ogrtest.check_feature_geometry(
+        out_f, "POLYGON ((0.2 0.8,0.7 0.8,0.7 0.3,0.2 0.3,0.2 0.8))"
+    )
+
+    assert out_lyr.GetNextFeature() is None
+
+
+def test_gdalalg_vector_clip_bbox_srs():
+
+    src_ds = gdal.GetDriverByName("Memory").Create("", 0, 0, 0, gdal.GDT_Unknown)
+    srs_wgs84 = osr.SpatialReference()
+    srs_wgs84.SetFromUserInput("WGS84")
+    srs_wgs84.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
+    src_lyr = src_ds.CreateLayer("test", srs=srs_wgs84)
+    src_lyr.CreateField(ogr.FieldDefn("foo"))
+
+    f = ogr.Feature(src_lyr.GetLayerDefn())
+    f["foo"] = "bar"
+    f.SetGeometry(ogr.CreateGeometryFromWkt("POLYGON ((0 0,0 1,1 1,1 0,0 0))"))
+    src_lyr.CreateFeature(f)
+
+    clip = get_clip_alg()
+    clip.GetArg("input").Get().SetDataset(src_ds)
+
+    assert clip.ParseCommandLineArguments(
+        [
+            "--bbox",
+            "22263.8981586547,33395.9998333802,77923.6435552915,89058.4864167096",
+            "--bbox-crs",
+            "EPSG:3857",
+            "--of",
+            "Memory",
+            "--output",
+            "memory_ds",
+        ]
+    )
+    assert clip.Run()
+
+    out_ds = clip.GetArg("output").Get().GetDataset()
+    out_lyr = out_ds.GetLayer(0)
+    assert out_lyr.GetSpatialRef().GetAuthorityCode(None) == "4326"
+    out_f = out_lyr.GetNextFeature()
+    assert out_f["foo"] == "bar"
+    ogrtest.check_feature_geometry(
+        out_f, "POLYGON ((0.2 0.8,0.7 0.8,0.7 0.3,0.2 0.3,0.2 0.8))"
+    )
+
+    assert out_lyr.GetNextFeature() is None
+
+
+def test_gdalalg_vector_clip_split_multipart():
+
+    src_ds = gdal.GetDriverByName("Memory").Create("", 0, 0, 0, gdal.GDT_Unknown)
+    src_lyr = src_ds.CreateLayer("test", geom_type=ogr.wkbPolygon)
+    src_lyr.CreateField(ogr.FieldDefn("foo"))
+
+    f = ogr.Feature(src_lyr.GetLayerDefn())
+    f["foo"] = "bar"
+    f.SetGeometry(
+        ogr.CreateGeometryFromWkt(
+            "POLYGON ((0 0,0 1,1 1,1 0,0.8 0,0.8 0.8,0.2 0.8,0.2 0,0 0))"
+        )
+    )
+    src_lyr.CreateFeature(f)
+
+    f = ogr.Feature(src_lyr.GetLayerDefn())
+    f["foo"] = "baz"
+    f.SetGeometry(
+        ogr.CreateGeometryFromWkt("POLYGON ((0.4 0.4,0.4 0.6,0.6 0.6,0.6 0.4,0.4 0.4))")
+    )
+    src_lyr.CreateFeature(f)
+
+    clip = get_clip_alg()
+    clip.GetArg("input").Get().SetDataset(src_ds)
+
+    assert clip.ParseCommandLineArguments(
+        ["--bbox", "-0.1,0.3,1.1,0.7", "--of", "Memory", "--output", "memory_ds"]
+    )
+    assert clip.Run()
+
+    out_ds = clip.GetArg("output").Get().GetDataset()
+    out_lyr = out_ds.GetLayer(0)
+
+    ref1_geom = ogr.CreateGeometryFromWkt(
+        "POLYGON ((0.0 0.7,0.2 0.7,0.2 0.3,0.0 0.3,0.0 0.7))"
+    )
+    ref2_geom = ogr.CreateGeometryFromWkt(
+        "POLYGON ((1.0 0.3,0.8 0.3,0.8 0.7,1.0 0.7,1.0 0.3))"
+    )
+
+    out_f = out_lyr.GetNextFeature()
+    assert out_f["foo"] == "bar"
+    g = out_f.GetGeometryRef()
+    was_one = g.Within(ref1_geom) and ref1_geom.Within(g)
+    assert was_one or (g.Within(ref2_geom) and ref2_geom.Within(g))
+
+    out_f = out_lyr.GetNextFeature()
+    assert out_f["foo"] == "bar"
+    g = out_f.GetGeometryRef()
+    if was_one:
+        assert g.Within(ref2_geom) and ref2_geom.Within(g)
+    else:
+        assert g.Within(ref1_geom) and ref1_geom.Within(g)
+
+    out_f = out_lyr.GetNextFeature()
+    assert out_f["foo"] == "baz"
+    ogrtest.check_feature_geometry(
+        out_f, "POLYGON ((0.4 0.4,0.4 0.6,0.6 0.6,0.6 0.4,0.4 0.4))"
+    )
+
+    assert out_lyr.GetNextFeature() is None
+
+
+@pytest.mark.parametrize(
+    "clip_geom",
+    [
+        "POLYGON ((0.2 0.8,0.7 0.8,0.7 0.3,0.2 0.3,0.2 0.8))",
+        '{"type":"Polygon","coordinates":[[[0.2,0.8],[0.7,0.8],[0.7,0.3],[0.2,0.3],[0.2,0.8]]]}',
+    ],
+)
+def test_gdalalg_vector_clip_geom(clip_geom):
+
+    src_ds = gdal.GetDriverByName("Memory").Create("", 0, 0, 0, gdal.GDT_Unknown)
+    src_lyr = src_ds.CreateLayer("test")
+    src_lyr.CreateField(ogr.FieldDefn("foo"))
+
+    f = ogr.Feature(src_lyr.GetLayerDefn())
+    f["foo"] = "bar"
+    f.SetGeometry(ogr.CreateGeometryFromWkt("POLYGON ((0 0,0 -1,-1 -1,-1 0,0 0))"))
+    src_lyr.CreateFeature(f)
+
+    f = ogr.Feature(src_lyr.GetLayerDefn())
+    f["foo"] = "bar"
+    f.SetGeometry(ogr.CreateGeometryFromWkt("POLYGON ((0 0,0 1,1 1,1 0,0 0))"))
+    src_lyr.CreateFeature(f)
+
+    clip = get_clip_alg()
+    clip.GetArg("input").Get().SetDataset(src_ds)
+
+    assert clip.ParseCommandLineArguments(
+        ["--geometry", clip_geom, "--of", "Memory", "--output", "memory_ds"]
+    )
+    assert clip.Run()
+
+    out_ds = clip.GetArg("output").Get().GetDataset()
+    out_lyr = out_ds.GetLayer(0)
+    out_f = out_lyr.GetNextFeature()
+    assert out_f["foo"] == "bar"
+    ogrtest.check_feature_geometry(
+        out_f, "POLYGON ((0.2 0.8,0.7 0.8,0.7 0.3,0.2 0.3,0.2 0.8))"
+    )
+
+    assert out_lyr.GetNextFeature() is None
+
+
+@pytest.mark.parametrize(
+    "clip_geom",
+    [
+        "POLYGON ((0.2 0.8,0.7 0.8,0.7 0.3,0.2 0.3,0.2 0.8))",
+        '{"type":"Polygon","coordinates":[[[0.2,0.8],[0.7,0.8],[0.7,0.3],[0.2,0.3],[0.2,0.8]]]}',
+    ],
+)
+def test_gdalalg_vector_clip_geom_srs(clip_geom):
+
+    src_ds = gdal.GetDriverByName("Memory").Create("", 0, 0, 0, gdal.GDT_Unknown)
+    srs_wgs84 = osr.SpatialReference()
+    srs_wgs84.SetFromUserInput("WGS84")
+    srs_wgs84.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
+    src_lyr = src_ds.CreateLayer("test", srs=srs_wgs84)
+    src_lyr.CreateField(ogr.FieldDefn("foo"))
+
+    f = ogr.Feature(src_lyr.GetLayerDefn())
+    f["foo"] = "bar"
+    f.SetGeometry(ogr.CreateGeometryFromWkt("POLYGON ((0 0,0 1,1 1,1 0,0 0))"))
+    src_lyr.CreateFeature(f)
+
+    clip = get_clip_alg()
+    clip.GetArg("input").Get().SetDataset(src_ds)
+    clip_geom = clip_geom.replace("0.2", "22263.8981586547")
+    clip_geom = clip_geom.replace("0.3", "33395.9998333802")
+    clip_geom = clip_geom.replace("0.7", "77923.6435552915")
+    clip_geom = clip_geom.replace("0.8", "89058.4864167096")
+    assert clip.ParseCommandLineArguments(
+        [
+            "--geometry",
+            clip_geom,
+            "--geometry-crs",
+            "EPSG:3857",
+            "--of",
+            "Memory",
+            "--output",
+            "memory_ds",
+        ]
+    )
+    assert clip.Run()
+
+    out_ds = clip.GetArg("output").Get().GetDataset()
+    out_lyr = out_ds.GetLayer(0)
+    out_f = out_lyr.GetNextFeature()
+    assert out_f["foo"] == "bar"
+    ogrtest.check_feature_geometry(
+        ogr.CreateGeometryFromWkt(out_f.GetGeometryRef().ExportToWkt()),
+        "POLYGON ((0.2 0.8,0.7 0.8,0.7 0.3,0.2 0.3,0.2 0.8))",
+    )
+
+    assert out_lyr.GetNextFeature() is None
+
+
+def test_gdalalg_vector_clip_like_vector():
+
+    src_ds = gdal.GetDriverByName("Memory").Create("", 0, 0, 0, gdal.GDT_Unknown)
+    src_lyr = src_ds.CreateLayer("test")
+    src_lyr.CreateField(ogr.FieldDefn("foo"))
+    f = ogr.Feature(src_lyr.GetLayerDefn())
+    f["foo"] = "bar"
+    f.SetGeometry(ogr.CreateGeometryFromWkt("POLYGON ((0 0,0 1,1 1,1 0,0 0))"))
+    src_lyr.CreateFeature(f)
+
+    # Create "like" dataset
+    like_ds = gdal.GetDriverByName("Memory").Create("", 0, 0, 0, gdal.GDT_Unknown)
+    like_lyr = like_ds.CreateLayer("test")
+    f = ogr.Feature(like_lyr.GetLayerDefn())
+    f.SetGeometry(
+        ogr.CreateGeometryFromWkt("POLYGON ((0.2 0.8,0.7 0.8,0.7 0.3,0.2 0.3,0.2 0.8))")
+    )
+    like_lyr.CreateFeature(f)
+
+    clip = get_clip_alg()
+    clip.GetArg("input").Get().SetDataset(src_ds)
+    clip.GetArg("like").Get().SetDataset(like_ds)
+
+    assert clip.ParseCommandLineArguments(["--of", "Memory", "--output", "memory_ds"])
+    assert clip.Run()
+
+    out_ds = clip.GetArg("output").Get().GetDataset()
+    out_lyr = out_ds.GetLayer(0)
+    out_f = out_lyr.GetNextFeature()
+    assert out_f["foo"] == "bar"
+    ogrtest.check_feature_geometry(
+        out_f, "POLYGON ((0.2 0.8,0.7 0.8,0.7 0.3,0.2 0.3,0.2 0.8))"
+    )
+
+    assert out_lyr.GetNextFeature() is None
+
+
+def test_gdalalg_vector_clip_like_vector_srs():
+
+    src_ds = gdal.GetDriverByName("Memory").Create("", 0, 0, 0, gdal.GDT_Unknown)
+    srs_wgs84 = osr.SpatialReference()
+    srs_wgs84.SetFromUserInput("WGS84")
+    srs_wgs84.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
+    src_lyr = src_ds.CreateLayer("test", srs=srs_wgs84)
+    f = ogr.Feature(src_lyr.GetLayerDefn())
+    f.SetGeometry(ogr.CreateGeometryFromWkt("POLYGON ((0 0,0 1,1 1,1 0,0 0))"))
+    src_lyr.CreateFeature(f)
+
+    # Create "like" dataset
+    like_ds = gdal.GetDriverByName("Memory").Create("", 0, 0, 0, gdal.GDT_Unknown)
+    srs_3857 = osr.SpatialReference()
+    srs_3857.SetFromUserInput("EPSG:3857")
+    srs_3857.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
+    like_lyr = like_ds.CreateLayer("test", srs=srs_3857)
+    f = ogr.Feature(like_lyr.GetLayerDefn())
+    wkt = "POLYGON ((0.2 0.8,0.7 0.8,0.7 0.3,0.2 0.3,0.2 0.8))"
+    wkt = wkt.replace("0.2", "22263.8981586547")
+    wkt = wkt.replace("0.3", "33395.9998333802")
+    wkt = wkt.replace("0.7", "77923.6435552915")
+    wkt = wkt.replace("0.8", "89058.4864167096")
+    f.SetGeometry(ogr.CreateGeometryFromWkt(wkt))
+    like_lyr.CreateFeature(f)
+
+    clip = get_clip_alg()
+    clip.GetArg("input").Get().SetDataset(src_ds)
+    clip.GetArg("like").Get().SetDataset(like_ds)
+
+    assert clip.ParseCommandLineArguments(["--of", "Memory", "--output", "memory_ds"])
+    assert clip.Run()
+
+    out_ds = clip.GetArg("output").Get().GetDataset()
+    out_lyr = out_ds.GetLayer(0)
+    out_f = out_lyr.GetNextFeature()
+    ogrtest.check_feature_geometry(
+        out_f, "POLYGON ((0.2 0.8,0.7 0.8,0.7 0.3,0.2 0.3,0.2 0.8))"
+    )
+
+    assert out_lyr.GetNextFeature() is None
+
+
+def test_gdalalg_vector_clip_like_raster():
+
+    src_ds = gdal.GetDriverByName("Memory").Create("", 0, 0, 0, gdal.GDT_Unknown)
+    src_lyr = src_ds.CreateLayer("test")
+    src_lyr.CreateField(ogr.FieldDefn("foo"))
+    f = ogr.Feature(src_lyr.GetLayerDefn())
+    f["foo"] = "bar"
+    f.SetGeometry(ogr.CreateGeometryFromWkt("POLYGON ((0 0,0 1,1 1,1 0,0 0))"))
+    src_lyr.CreateFeature(f)
+
+    # Create "like" dataset
+    like_ds = gdal.GetDriverByName("MEM").Create("", 1, 1, 1)
+    like_ds.SetGeoTransform([0.2, 0.5, 0, 0.3, 0, 0.5])
+
+    clip = get_clip_alg()
+    clip.GetArg("input").Get().SetDataset(src_ds)
+    clip.GetArg("like").Get().SetDataset(like_ds)
+
+    assert clip.ParseCommandLineArguments(["--of", "Memory", "--output", "memory_ds"])
+    assert clip.Run()
+
+    out_ds = clip.GetArg("output").Get().GetDataset()
+    out_lyr = out_ds.GetLayer(0)
+    out_f = out_lyr.GetNextFeature()
+    assert out_f["foo"] == "bar"
+    ogrtest.check_feature_geometry(
+        out_f, "POLYGON ((0.2 0.8,0.7 0.8,0.7 0.3,0.2 0.3,0.2 0.8))"
+    )
+
+    assert out_lyr.GetNextFeature() is None
+
+
+def test_gdalalg_vector_clip_like_raster_srs():
+
+    src_ds = gdal.GetDriverByName("Memory").Create("", 0, 0, 0, gdal.GDT_Unknown)
+    srs_wgs84 = osr.SpatialReference()
+    srs_wgs84.SetFromUserInput("WGS84")
+    srs_wgs84.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
+    src_lyr = src_ds.CreateLayer("test", srs=srs_wgs84)
+    src_lyr.CreateField(ogr.FieldDefn("foo"))
+    f = ogr.Feature(src_lyr.GetLayerDefn())
+    f["foo"] = "bar"
+    f.SetGeometry(ogr.CreateGeometryFromWkt("POLYGON ((0 0,0 1,1 1,1 0,0 0))"))
+    src_lyr.CreateFeature(f)
+
+    # Create "like" dataset
+    like_ds = gdal.GetDriverByName("MEM").Create("", 1, 1, 1)
+    like_ds.SetGeoTransform(
+        [22263.8981586547, 55660.4518654215, 0, 33395.9998333802, 0, 55659.7453966368]
+    )
+    srs_3857 = osr.SpatialReference()
+    srs_3857.SetFromUserInput("EPSG:3857")
+    srs_3857.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
+    like_ds.SetSpatialRef(srs_3857)
+
+    clip = get_clip_alg()
+    clip.GetArg("input").Get().SetDataset(src_ds)
+    clip.GetArg("like").Get().SetDataset(like_ds)
+
+    assert clip.ParseCommandLineArguments(["--of", "Memory", "--output", "memory_ds"])
+    assert clip.Run()
+
+    out_ds = clip.GetArg("output").Get().GetDataset()
+    out_lyr = out_ds.GetLayer(0)
+    out_f = out_lyr.GetNextFeature()
+    assert out_f["foo"] == "bar"
+    ogrtest.check_feature_geometry(
+        out_f, "POLYGON ((0.2 0.8,0.7 0.8,0.7 0.3,0.2 0.3,0.2 0.8))"
+    )
+
+    assert out_lyr.GetNextFeature() is None
+
+
+def test_gdalalg_vector_clip_missing_arg(tmp_vsimem):
+
+    out_filename = str(tmp_vsimem / "out.shp")
+
+    clip = get_clip_alg()
+    with pytest.raises(
+        Exception, match="clip: --bbox, --geometry or --like must be specified"
+    ):
+        clip.ParseRunAndFinalize(["../ogr/data/poly.shp", out_filename])
+
+
+def test_gdalalg_vector_clip_geometry_invalid():
+
+    src_ds = gdal.GetDriverByName("Memory").Create("", 0, 0, 0, gdal.GDT_Unknown)
+    src_ds.CreateLayer("test")
+
+    clip = get_clip_alg()
+    clip.GetArg("input").Get().SetDataset(src_ds)
+
+    assert clip.ParseCommandLineArguments(
+        ["--geometry", "invalid", "--of", "Memory", "--output", "memory_ds"]
+    )
+    with pytest.raises(
+        Exception,
+        match="clip: Clipping geometry is neither a valid WKT or GeoJSON geometry",
+    ):
+        clip.Run()
+
+
+def test_gdalalg_vector_clip_like_vector_too_many_layers():
+
+    src_ds = gdal.GetDriverByName("Memory").Create("", 0, 0, 0, gdal.GDT_Unknown)
+    src_ds.CreateLayer("test")
+
+    # Create "like" dataset
+    like_ds = gdal.GetDriverByName("MEM").Create("", 1, 1, 1)
+
+    clip = get_clip_alg()
+    clip.GetArg("input").Get().SetDataset(src_ds)
+    clip.GetArg("like").Get().SetDataset(like_ds)
+
+    assert clip.ParseCommandLineArguments(["--of", "Memory", "--output", "memory_ds"])
+    with pytest.raises(
+        Exception,
+        match="clip: Dataset '' has no geotransform matrix. Its bounds cannot be established",
+    ):
+        clip.Run()
+
+
+def test_gdalalg_vector_clip_like_raster_no_geotransform():
+
+    src_ds = gdal.GetDriverByName("Memory").Create("", 0, 0, 0, gdal.GDT_Unknown)
+    src_ds.CreateLayer("test")
+
+    # Create "like" dataset
+    like_ds = gdal.GetDriverByName("Memory").Create("", 0, 0, 0, gdal.GDT_Unknown)
+    like_ds.CreateLayer("test1")
+    like_ds.CreateLayer("test2")
+
+    clip = get_clip_alg()
+    clip.GetArg("input").Get().SetDataset(src_ds)
+    clip.GetArg("like").Get().SetDataset(like_ds)
+
+    assert clip.ParseCommandLineArguments(["--of", "Memory", "--output", "memory_ds"])
+    with pytest.raises(
+        Exception, match="clip: Only single layer dataset can be specified with --like"
+    ):
+        clip.Run()
diff --git a/doc/source/conf.py b/doc/source/conf.py
index cb9096e2d997..51be7d65f45c 100644
--- a/doc/source/conf.py
+++ b/doc/source/conf.py
@@ -292,6 +292,13 @@
         [author_evenr],
         1,
     ),
+    (
+        "programs/gdal_vector_clip",
+        "gdal-vector-clip",
+        "Clip a vector dataset",
+        [author_evenr],
+        1,
+    ),
     (
         "programs/gdal_vector_convert",
         "gdal-vector-convert",
diff --git a/doc/source/programs/gdal_options/co_vector.rst b/doc/source/programs/gdal_options/co_vector.rst
new file mode 100644
index 000000000000..329bb104623a
--- /dev/null
+++ b/doc/source/programs/gdal_options/co_vector.rst
@@ -0,0 +1,18 @@
+.. option:: --co <NAME>=<VALUE>
+
+    Many formats have one or more optional dataset creation options that can be
+    used to control particulars about the file created. For instance,
+    the GeoPackage driver supports creation options to control the version.
+
+    May be repeated.
+
+    The dataset creation options available vary by format driver, and some
+    simple formats have no creation options at all. A list of options
+    supported for a format can be listed with the
+    :ref:`--formats <vector_common_options>`
+    command line option but the documentation for the format is the
+    definitive source of information on driver creation options.
+    See :ref:`vector_drivers` format
+    specific documentation for legal creation options for each format.
+
+    Note that dataset creation options are different from layer creation options.
diff --git a/doc/source/programs/gdal_vector.rst b/doc/source/programs/gdal_vector.rst
index 49d1a8ae48a7..9b503db16b7a 100644
--- a/doc/source/programs/gdal_vector.rst
+++ b/doc/source/programs/gdal_vector.rst
@@ -19,6 +19,7 @@ Synopsis
 
     Usage: gdal vector <SUBCOMMAND>
     where <SUBCOMMAND> is one of:
+      - clip:      Clip a vector dataset.
       - convert:   Convert a vector dataset.
       - filter:    Filter a vector dataset.
       - info:      Return information on a vector dataset.
@@ -28,8 +29,9 @@ Synopsis
 Available sub-commands
 ----------------------
 
-- :ref:`gdal_vector_info_subcommand`
+- :ref:`gdal_vector_clip_subcommand`
 - :ref:`gdal_vector_convert_subcommand`
+- :ref:`gdal_vector_info_subcommand`
 - :ref:`gdal_vector_pipeline_subcommand`
 
 Examples
diff --git a/doc/source/programs/gdal_vector_clip.rst b/doc/source/programs/gdal_vector_clip.rst
new file mode 100644
index 000000000000..71c857ad469d
--- /dev/null
+++ b/doc/source/programs/gdal_vector_clip.rst
@@ -0,0 +1,132 @@
+.. _gdal_vector_clip_subcommand:
+
+================================================================================
+"gdal vector clip" sub-command
+================================================================================
+
+.. versionadded:: 3.11
+
+.. only:: html
+
+    Clip a vector dataset.
+
+.. Index:: gdal vector clip
+
+Synopsis
+--------
+
+.. code-block::
+
+    Usage: gdal vector clip [OPTIONS] <INPUT> <OUTPUT>
+
+    Clip a vector dataset.
+
+    Positional arguments:
+      -i, --input <INPUT>                                  Input vector dataset [required]
+      -o, --output <OUTPUT>                                Output vector dataset [required]
+
+    Common Options:
+      -h, --help                                           Display help message and exit
+      --version                                            Display GDAL version and exit
+      --json-usage                                         Display usage as JSON document and exit
+      --drivers                                            Display driver list as JSON document and exit
+      --config <KEY>=<VALUE>                               Configuration option [may be repeated]
+      --progress                                           Display progress bar
+
+    Options:
+      -l, --layer, --input-layer <INPUT-LAYER>             Input layer name(s) [may be repeated]
+      -f, --of, --format, --output-format <OUTPUT-FORMAT>  Output format
+      --co, --creation-option <KEY>=<VALUE>                Creation option [may be repeated]
+      --lco, --layer-creation-option <KEY>=<VALUE>         Layer creation option [may be repeated]
+      --overwrite                                          Whether overwriting existing output is allowed
+      --update                                             Whether to open existing dataset in update mode
+      --overwrite-layer                                    Whether overwriting existing layer is allowed
+      --append                                             Whether appending to existing layer is allowed
+      --output-layer <OUTPUT-LAYER>                        Output layer name
+      --bbox <BBOX>                                        Clipping bounding box as xmin,ymin,xmax,ymax
+                                                           Mutually exclusive with --geometry, --like
+      --bbox-crs <BBOX-CRS>                                CRS of clipping bounding box
+      --geometry <GEOMETRY>                                Clipping geometry (WKT or GeoJSON)
+                                                           Mutually exclusive with --bbox, --like
+      --geometry-crs <GEOMETRY-CRS>                        CRS of clipping geometry
+      --like <DATASET>                                     Dataset to use as a template for bounds
+                                                           Mutually exclusive with --bbox, --geometry
+
+    Advanced Options:
+      --if, --input-format <INPUT-FORMAT>                  Input formats [may be repeated]
+      --oo, --open-option <KEY=VALUE>                      Open options [may be repeated]
+
+
+
+Description
+-----------
+
+:program:`gdal vector clip` can be used to clip a vector dataset using
+georeferenced coordinates.
+
+Either :option:`--bbox`, :option:`--geometry` or :option:`--like` must be specified.
+
+``clip`` can also be used as a step of :ref:`gdal_vector_pipeline_subcommand`.
+
+Standard options
+++++++++++++++++
+
+.. include:: gdal_options/of_vector.rst
+
+.. include:: gdal_options/co_vector.rst
+
+.. include:: gdal_options/overwrite.rst
+
+.. option:: --bbox <xmin>,<ymin>,<xmax>,<ymax>
+
+    Bounds to which to clip the dataset. They are assumed to be in the CRS of
+    the input dataset, unless :option:`--bbox-crs` is specified.
+    The X and Y axis are the "GIS friendly ones", that is X is longitude or easting,
+    and Y is latitude or northing.
+
+.. option:: --bbox-crs <CRS>
+
+    CRS in which the <xmin>,<ymin>,<xmax>,<ymax> values of :option:`--bbox`
+    are expressed. If not specified, it is assumed to be the CRS of the input
+    dataset.
+    Not that specifying --bbox-crs does not involve doing vector reprojection.
+    Instead, the bounds are reprojected from the bbox-crs to the CRS of the
+    input dataset.
+
+.. option:: --geometry <WKT_or_GeoJSON>
+
+    Geometry as a WKT or GeoJSON string to which to clip the dataset.
+    It is assumed to be in the CRS of
+    the input dataset, unless :option:`--geometry-crs` is specified.
+    The X and Y axis are the "GIS friendly ones", that is X is longitude or easting,
+    and Y is latitude or northing.
+
+.. option:: --geometry-crs <CRS>
+
+    CRS in which the coordinates values of :option:`--geometry`
+    are expressed. If not specified, it is assumed to be the CRS of the input
+    dataset.
+    Not that specifying --geometry-crs does not involve doing vector reprojection.
+    Instead, the bounds are reprojected from the geometry-crs to the CRS of the
+    input dataset.
+
+.. option:: --like <DATASET>
+
+    Vector or raster dataset to use as a template for bounds.
+
+Advanced options
+++++++++++++++++
+
+.. include:: gdal_options/oo.rst
+
+.. include:: gdal_options/if.rst
+
+Examples
+--------
+
+.. example::
+   :title: Clip a GeoPackage file to the bounding box from longitude 2, latitude 49, to longitude 3, latitude 50 in WGS 84
+
+   .. code-block:: bash
+
+        $ gdal vector clip --bbox=2,49,3,50 --bbox-crs=EPSG:4326 in.gpkg out.gpkg --overwrite
diff --git a/doc/source/programs/gdal_vector_convert.rst b/doc/source/programs/gdal_vector_convert.rst
index fdd736101b28..f475ecccc818 100644
--- a/doc/source/programs/gdal_vector_convert.rst
+++ b/doc/source/programs/gdal_vector_convert.rst
@@ -59,7 +59,7 @@ Standard options
 
 .. include:: gdal_options/of_vector.rst
 
-.. include:: gdal_options/co.rst
+.. include:: gdal_options/co_vector.rst
 
 .. include:: gdal_options/overwrite.rst
 
diff --git a/doc/source/programs/gdal_vector_pipeline.rst b/doc/source/programs/gdal_vector_pipeline.rst
index 394ff0d092ae..1b5fb2d82d71 100644
--- a/doc/source/programs/gdal_vector_pipeline.rst
+++ b/doc/source/programs/gdal_vector_pipeline.rst
@@ -52,6 +52,24 @@ Potential steps are:
       --if, --input-format <INPUT-FORMAT>                  Input formats [may be repeated]
       --oo, --open-option <KEY=VALUE>                      Open options [may be repeated]
 
+* clip [OPTIONS]
+
+.. code-block::
+
+    Clip a vector dataset.
+
+    Options:
+      --bbox <BBOX>                                        Clipping bounding box as xmin,ymin,xmax,ymax
+                                                           Mutually exclusive with --geometry, --like
+      --bbox-crs <BBOX-CRS>                                CRS of clipping bounding box
+      --geometry <GEOMETRY>                                Clipping geometry (WKT or GeoJSON)
+                                                           Mutually exclusive with --bbox, --like
+      --geometry-crs <GEOMETRY-CRS>                        CRS of clipping geometry
+      --like <DATASET>                                     Dataset to use as a template for bounds
+                                                           Mutually exclusive with --bbox, --geometry
+
+Details for options can be found in :ref:`gdal_vector_clip_subcommand`.
+
 * filter [OPTIONS]
 
 .. code-block::
diff --git a/doc/source/programs/index.rst b/doc/source/programs/index.rst
index 8fd8d3ab2945..aebb8a2c34b6 100644
--- a/doc/source/programs/index.rst
+++ b/doc/source/programs/index.rst
@@ -41,6 +41,7 @@ single :program:`gdal` program that accepts commands and subcommands.
    gdal_raster_reproject
    gdal_vector
    gdal_vector_info
+   gdal_vector_clip
    gdal_vector_convert
    gdal_vector_pipeline
 
@@ -62,6 +63,7 @@ single :program:`gdal` program that accepts commands and subcommands.
     - :ref:`gdal_raster_reproject_subcommand`: Reproject a raster dataset
     - :ref:`gdal_vector_command`: Entry point for vector commands
     - :ref:`gdal_vector_info_subcommand`: Get information on a vector dataset
+    - :ref:`gdal_vector_clip_subcommand`: Clip a vector dataset
     - :ref:`gdal_vector_convert_subcommand`: Convert a vector dataset
     - :ref:`gdal_vector_pipeline_subcommand`: Process a vector dataset