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(); + RegisterSubAlgorithm(); RegisterSubAlgorithm(); RegisterSubAlgorithm(); RegisterSubAlgorithm(); diff --git a/apps/gdalalg_vector_clip.cpp b/apps/gdalalg_vector_clip.cpp new file mode 100644 index 000000000000..8cd6edd5d123 --- /dev/null +++ b/apps/gdalalg_vector_clip.cpp @@ -0,0 +1,466 @@ +/****************************************************************************** + * + * Project: GDAL + * Purpose: "clip" step of "vector pipeline", or "gdal vector clip" standalone + * Author: Even Rouault + * + ****************************************************************************** + * Copyright (c) 2025, Even Rouault + * + * SPDX-License-Identifier: MIT + ****************************************************************************/ + +#include "gdalalg_vector_clip.h" + +#include "gdal_priv.h" +#include "gdal_utils.h" +#include "ogrsf_frmts.h" + +#include + +//! @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("bbox-geometry-like"); + 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("bbox-geometry-like"); + 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("bbox-geometry-like"); + AddArg("like-sql", 0, ("SELECT statement to run on the 'like' dataset"), + &m_likeSQL) + .SetMetaVar("SELECT-STATEMENT") + .SetMutualExclusionGroup("sql-where"); + AddArg("like-layer", 0, ("Name of the layer of the 'like' dataset"), + &m_likeLayer) + .SetMetaVar("LAYER-NAME"); + AddArg("like-where", 0, ("WHERE SQL clause to run on the 'like' dataset"), + &m_likeWhere) + .SetMetaVar("WHERE-EXPRESSION") + .SetMutualExclusionGroup("sql-where"); +} + +/************************************************************************/ +/* GDALVectorClipAlgorithmDataset */ +/************************************************************************/ + +namespace +{ +class GDALVectorClipAlgorithmDataset final : public GDALDataset +{ + std::vector> m_layers{}; + + public: + GDALVectorClipAlgorithmDataset() = default; + + void AddLayer(std::unique_ptr poLayer) + { + m_layers.push_back(std::move(poLayer)); + } + + int GetLayerCount() override + { + return static_cast(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 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(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(m_poSrcLayer->GetNextFeature())) + { + auto poGeom = poFeature->GetGeometryRef(); + if (!poGeom) + continue; + + auto poIntersection = std::unique_ptr( + 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 m_poClipGeom{}; + const OGRwkbGeometryType m_eSrcLayerGeomType; + const bool m_bSrcLayerGeomTypeIsCollection; + std::unique_ptr m_poSrcFeature{}; + std::unique_ptr m_poCurGeomColl{}; + int m_idxInCurGeomColl = 0; + + CPL_DISALLOW_COPY_ASSIGN(GDALVectorClipAlgorithmLayer) +}; + +} // namespace + +/************************************************************************/ +/* LoadGeometry() */ +/************************************************************************/ + +static std::unique_ptr LoadGeometry(GDALDataset *poDS, + const std::string &osSQL, + const std::string &osLyr, + const std::string &osWhere) +{ + OGRLayer *poLyr = nullptr; + if (!osSQL.empty()) + poLyr = poDS->ExecuteSQL(osSQL.c_str(), nullptr, nullptr); + else if (!osLyr.empty()) + poLyr = poDS->GetLayerByName(osLyr.c_str()); + else + poLyr = poDS->GetLayer(0); + + if (poLyr == nullptr) + { + CPLError(CE_Failure, CPLE_AppDefined, + "Failed to identify source layer from clipping dataset."); + return nullptr; + } + + if (!osWhere.empty()) + poLyr->SetAttributeFilter(osWhere.c_str()); + + OGRGeometryCollection oGC; + + const auto poSRSSrc = poLyr->GetSpatialRef(); + if (poSRSSrc) + { + auto poSRSClone = poSRSSrc->Clone(); + oGC.assignSpatialReference(poSRSClone); + poSRSClone->Release(); + } + + for (auto &poFeat : poLyr) + { + auto poSrcGeom = std::unique_ptr(poFeat->StealGeometry()); + if (poSrcGeom) + { + // Only take into account areal geometries. + if (poSrcGeom->getDimension() == 2) + { + if (!poSrcGeom->IsValid()) + { + CPLError(CE_Failure, CPLE_AppDefined, + "Geometry of feature " CPL_FRMT_GIB " of %s " + "is invalid.", + poFeat->GetFID(), poDS->GetDescription()); + return nullptr; + } + else + { + oGC.addGeometry(std::move(poSrcGeom)); + } + } + } + } + + if (!osSQL.empty()) + poDS->ReleaseResultSet(poLyr); + + if (oGC.IsEmpty()) + { + CPLError(CE_Failure, CPLE_AppDefined, "No clipping geometry found"); + return nullptr; + } + + return std::unique_ptr(oGC.UnaryUnion()); +} + +/************************************************************************/ +/* 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 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(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 && m_likeLayer.empty() && + m_likeSQL.empty()) + { + ReportError( + CE_Failure, CPLE_AppDefined, + "Only single layer dataset can be specified with --like when " + "neither --like-layer or --like-sql have been specified"); + return false; + } + else if (poLikeDS->GetLayerCount() > 0) + { + poClipGeom = + LoadGeometry(poLikeDS, m_likeSQL, m_likeLayer, m_likeWhere); + if (!poClipGeom) + return false; + } + 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(); + auto poLR = std::make_unique(); + 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; + } + + if (!poClipGeom) + { + m_outputDataset.Set(m_inputDataset.GetDatasetRef()); + return true; + } + + auto outDS = std::make_unique(); + 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(poClipGeom->clone()); + if (poClipGeomForLayer->getSpatialReference() && + poSrcLayer->GetSpatialRef()) + { + ret = poClipGeomForLayer->transformTo( + poSrcLayer->GetSpatialRef()) == OGRERR_NONE; + } + if (ret) + { + outDS->AddLayer(std::make_unique( + 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..edc30c0e5fd1 --- /dev/null +++ b/apps/gdalalg_vector_clip.h @@ -0,0 +1,67 @@ +/****************************************************************************** + * + * Project: GDAL + * Purpose: "clip" step of "vector pipeline", or "gdal vector clip" standalone + * Author: Even Rouault + * + ****************************************************************************** + * Copyright (c) 2025, Even Rouault + * + * 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 GetAliases() + { + return {}; + } + + explicit GDALVectorClipAlgorithm(bool standaloneStep = false); + + private: + bool RunStep(GDALProgressFunc pfnProgress, void *pProgressData) override; + + std::vector m_bbox{}; + std::string m_bboxCrs{}; + std::string m_geometry{}; + std::string m_geometryCrs{}; + GDALArgDatasetValue m_likeDataset{}; + std::string m_likeLayer{}; + std::string m_likeSQL{}; + std::string m_likeWhere{}; +}; + +/************************************************************************/ +/* 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(); m_stepRegistry.Register(); + m_stepRegistry.Register(); m_stepRegistry.Register(); m_stepRegistry.Register(); } diff --git a/autotest/utilities/test_gdalalg_vector_clip.py b/autotest/utilities/test_gdalalg_vector_clip.py new file mode 100755 index 000000000000..bd5a1398c8b2 --- /dev/null +++ b/autotest/utilities/test_gdalalg_vector_clip.py @@ -0,0 +1,839 @@ +#!/usr/bin/env pytest +# -*- coding: utf-8 -*- +############################################################################### +# Project: GDAL/OGR Test Suite +# Purpose: 'gdal vector clip' testing +# Author: Even Rouault +# +############################################################################### +# Copyright (c) 2025, Even Rouault +# +# 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()) + 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) + + 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_geom_not_rectangle(): + + 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" + # Bounding box of this polygon intersects the bounding box of the clipping + # geometry, but their full intersection is empty + f.SetGeometry( + ogr.CreateGeometryFromWkt("POLYGON ((0.4 0.4,0.4 0.5,0.5 0.5,0.4 0.4))") + ) + src_lyr.CreateFeature(f) + + f = ogr.Feature(src_lyr.GetLayerDefn()) + f["foo"] = "bar" + f.SetGeometry(ogr.CreateGeometryFromWkt("POLYGON ((-1 -1,-1 2,2 2,2 -1,-1 -1))")) + src_lyr.CreateFeature(f) + + clip = get_clip_alg() + clip.GetArg("input").Get().SetDataset(src_ds) + + assert clip.ParseCommandLineArguments( + [ + "--geometry", + "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))", + "--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 0,0 1,1 1,1 0,0.8 0,0.8 0.8,0.2 0.8,0.2 0,0 0))" + ) + + 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( + 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_invalid_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) + + # Create "like" dataset + like_ds = gdal.GetDriverByName("Memory").Create( + "clip_ds", 0, 0, 0, gdal.GDT_Unknown + ) + like_lyr = like_ds.CreateLayer("test") + f = ogr.Feature(like_lyr.GetLayerDefn()) + f.SetGeometry(ogr.CreateGeometryFromWkt("POLYGON ((0 0,1 1,0 1,1 0,0 0))")) + 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"]) + with gdal.quiet_errors(), pytest.raises( + Exception, match="Geometry of feature 0 of clip_ds is invalid" + ): + clip.Run() + + +def test_gdalalg_vector_clip_like_vector_like_layer(): + + 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_ds.CreateLayer("unused") + like_lyr = like_ds.CreateLayer("my_layer") + 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( + ["--like-layer", "my_layer", "--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_like_layer_invalid(): + + src_ds = gdal.GetDriverByName("Memory").Create("", 0, 0, 0, gdal.GDT_Unknown) + + # Create "like" dataset + like_ds = gdal.GetDriverByName("Memory").Create("", 0, 0, 0, gdal.GDT_Unknown) + like_ds.CreateLayer("test") + + clip = get_clip_alg() + clip.GetArg("input").Get().SetDataset(src_ds) + clip.GetArg("like").Get().SetDataset(like_ds) + + assert clip.ParseCommandLineArguments( + ["--like-layer", "invalid", "--of", "Memory", "--output", "memory_ds"] + ) + with pytest.raises( + Exception, match="Failed to identify source layer from clipping dataset." + ): + clip.Run() + + +def test_gdalalg_vector_clip_like_vector_like_sql(): + + 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_ds.CreateLayer("unused") + like_lyr = like_ds.CreateLayer("my_layer") + 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( + [ + "--like-sql", + "SELECT * FROM my_layer", + "--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_like_where(): + + 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("my_layer") + like_lyr.CreateField(ogr.FieldDefn("id", ogr.OFTInteger)) + f = ogr.Feature(like_lyr.GetLayerDefn()) + f["id"] = 0 + f.SetGeometry(ogr.CreateGeometryFromWkt("POLYGON ((0 0,0 1,1 1,1 0,0 0))")) + like_lyr.CreateFeature(f) + + f = ogr.Feature(like_lyr.GetLayerDefn()) + f["id"] = 1 + 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( + ["--like-where", "id=1", "--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_like_where_empty(): + + 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("my_layer") + like_lyr.CreateField(ogr.FieldDefn("id", ogr.OFTInteger)) + + clip = get_clip_alg() + clip.GetArg("input").Get().SetDataset(src_ds) + clip.GetArg("like").Get().SetDataset(like_ds) + + assert clip.ParseCommandLineArguments( + ["--like-where", "id=1", "--of", "Memory", "--output", "memory_ds"] + ) + with pytest.raises(Exception, match="No clipping geometry found"): + clip.Run() + + +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( + 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_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 = + + 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 ` + 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 where 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..7a7196a7e096 --- /dev/null +++ b/doc/source/programs/gdal_vector_clip.rst @@ -0,0 +1,160 @@ +.. _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] + + Clip a vector dataset. + + Positional arguments: + -i, --input Input vector dataset [required] + -o, --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 = Configuration option [may be repeated] + --progress Display progress bar + + Options: + -l, --layer, --input-layer Input layer name(s) [may be repeated] + -f, --of, --format, --output-format Output format + --co, --creation-option = Creation option [may be repeated] + --lco, --layer-creation-option = 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 name + --bbox Clipping bounding box as xmin,ymin,xmax,ymax + Mutually exclusive with --geometry, --like + --bbox-crs CRS of clipping bounding box + --geometry Clipping geometry (WKT or GeoJSON) + Mutually exclusive with --bbox, --like + --geometry-crs CRS of clipping geometry + --like Dataset to use as a template for bounds + Mutually exclusive with --bbox, --geometry + --like-sql SELECT statement to run on the 'like' dataset + --like-layer Name of the layer of the 'like' dataset + --like-where Expression for a WHERE SQL clause to run on the 'like' dataset + + + Advanced Options: + --if, --input-format Input formats [may be repeated] + --oo, --open-option 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 ,,, + + 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. + Mutually exclusive with :option:`--like` and :option:`--geometry`. + +.. option:: --bbox-crs + + CRS in which the ,,, 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 + + 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. + Mutually exclusive with :option:`--bbox` and :option:`--like`. + +.. option:: --geometry-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 + + Vector or raster dataset to use as a template for bounds. + If the specified dataset is a raster, its rectangular bounds are use as + the clipping geometry. + If the specified dataset is a vector dataset, its polygonal geometries + are unioned together to form the clipping geometry. If several layers are present, + :option:`--like-sql` or :option:`--like-layer` must be specified. + Mutually exclusive with :option:`--bbox` and :option:`--geometry`. + +.. option:: --like-sql + + Select desired geometries from the vector clip dataset using an SQL query. + e.g ``SELECT geom FROM my_layer WHERE country = 'France'``. + Mutually exclusive with :option:`--like-layer` and :option:`--like-where` + +.. option:: --like-layer + + Select the named layer from the vector clip dataset. + Mutually exclusive with :option:`--like-sql` + +.. option:: --like-where + + Restrict desired geometries from vector clip dataset layer based on an attribute query. + e.g ``country = 'France'``. + +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..20d6b756fccd 100644 --- a/doc/source/programs/gdal_vector_pipeline.rst +++ b/doc/source/programs/gdal_vector_pipeline.rst @@ -52,6 +52,28 @@ Potential steps are: --if, --input-format Input formats [may be repeated] --oo, --open-option Open options [may be repeated] +* clip [OPTIONS] + +.. code-block:: + + Clip a vector dataset. + + Options: + --bbox Clipping bounding box as xmin,ymin,xmax,ymax + Mutually exclusive with --geometry, --like + --bbox-crs CRS of clipping bounding box + --geometry Clipping geometry (WKT or GeoJSON) + Mutually exclusive with --bbox, --like + --geometry-crs CRS of clipping geometry + --like Dataset to use as a template for bounds + Mutually exclusive with --bbox, --geometry + --like-sql SELECT statement to run on the 'like' dataset + --like-layer Name of the layer of the 'like' dataset + --like-where Expression for a WHERE SQL clause to run on the 'like' dataset + + +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