/src/gdal/apps/gdalalg_raster_clip.cpp

Source
/******************************************************************************
 *
 * Project:  GDAL
 * Purpose:  "clip" step of "raster pipeline", or "gdal raster 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_raster_clip.h"

#include "gdal_priv.h"
#include "gdal_utils.h"

#include <algorithm>
#include <cmath>

//! @cond Doxygen_Suppress

#ifndef _
#define _(x) (x)
#endif

/************************************************************************/
/*          GDALRasterClipAlgorithm::GDALRasterClipAlgorithm()          */
/************************************************************************/

GDALRasterClipAlgorithm::GDALRasterClipAlgorithm(bool standaloneStep)
    : GDALRasterPipelineStepAlgorithm(NAME, DESCRIPTION, HELP_URL,
                                      standaloneStep)
{
    constexpr const char *EXCLUSION_GROUP = "bbox-window-geometry-like";
    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("window", 0, _("Raster window as col,line,width,height in pixels"),
           &m_window)
        .SetRepeatedArgAllowed(false)
        .SetMinCount(4)
        .SetMaxCount(4)
        .SetDisplayHintAboutRepetition(false)
        .SetMutualExclusionGroup(EXCLUSION_GROUP)
        .AddValidationAction(
            [this]()
            {
                CPLAssert(m_window.size() == 4);
                if (m_window[2] <= 0 || m_window[3] <= 0)
                {
                    CPLError(CE_Failure, CPLE_AppDefined,
                             "Value of 'window' should be "
                             "col,line,width,height with "
                             "width > 0 and height > 0");
                    return false;
                }
                return true;
            });

    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);
    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");
    AddArg("only-bbox", 0,
           _("For 'geometry' and 'like', only consider their bounding box"),
           &m_onlyBBOX);
    AddArg("allow-bbox-outside-source", 0,
           _("Allow clipping box to include pixels outside input dataset"),
           &m_allowExtentOutsideSource);
    // Note: this would be a use case for multiple exclusion groups because it is not
    //       compatible with "window"
    AddArg("add-alpha", 0,
           _("Adds an alpha mask band to the destination when the source "
             "raster have none."),
           &m_addAlpha);
}

/************************************************************************/
/*                  GDALRasterClipAlgorithm::RunStep()                  */
/************************************************************************/

bool GDALRasterClipAlgorithm::RunStep(GDALPipelineStepRunContext &)
{
    auto poSrcDS = m_inputDataset[0].GetDatasetRef();
    CPLAssert(poSrcDS);
    CPLAssert(m_outputDataset.GetName().empty());
    CPLAssert(!m_outputDataset.GetDatasetRef());

    if (!m_window.empty())
    {
        if (m_addAlpha)
        {
            ReportError(CE_Failure, CPLE_NotSupported,
                        "'alpha' argument is not supported with 'window'");
            return false;
        }

        CPLStringList aosOptions;
        aosOptions.AddString("-of");
        aosOptions.AddString("VRT");

        aosOptions.AddString("-srcwin");
        aosOptions.AddString(CPLSPrintf("%d", m_window[0]));
        aosOptions.AddString(CPLSPrintf("%d", m_window[1]));
        aosOptions.AddString(CPLSPrintf("%d", m_window[2]));
        aosOptions.AddString(CPLSPrintf("%d", m_window[3]));

        if (!m_allowExtentOutsideSource)
        {
            // Unless we've specifically allowed the bounding box to extend beyond
            // the source raster, raise an error.
            aosOptions.AddString("-epo");
        }

        GDALTranslateOptions *psOptions =
            GDALTranslateOptionsNew(aosOptions.List(), nullptr);

        GDALDatasetH hSrcDS = GDALDataset::ToHandle(poSrcDS);
        auto poRetDS = GDALDataset::FromHandle(
            GDALTranslate("", hSrcDS, psOptions, nullptr));
        GDALTranslateOptionsFree(psOptions);

        const bool bOK = poRetDS != nullptr;
        if (bOK)
        {
            m_outputDataset.Set(std::unique_ptr<GDALDataset>(poRetDS));
        }

        return bOK;
    }

    GDALGeoTransform gt;
    if (poSrcDS->GetGeoTransform(gt) != CE_None)
    {
        ReportError(
            CE_Failure, CPLE_NotSupported,
            "Clipping is not supported on a raster without a geotransform");
        return false;
    }
    if (!gt.IsAxisAligned())
    {
        ReportError(CE_Failure, CPLE_NotSupported,
                    "Clipping is not supported on a raster whose geotransform "
                    "has rotation terms");
        return false;
    }

    auto [poClipGeom, errMsg] = GetClipGeometry();
    if (!poClipGeom)
    {
        ReportError(CE_Failure, CPLE_AppDefined, "%s", errMsg.c_str());
        return false;
    }

    auto poLikeDS = m_likeDataset.GetDatasetRef();
    if (!poClipGeom->getSpatialReference() && poLikeDS &&
        poLikeDS->GetLayerCount() == 0)
    {
        ReportError(CE_Failure, CPLE_AppDefined,
                    "Dataset '%s' has no CRS. Its bounds cannot be used.",
                    poLikeDS->GetDescription());
        return false;
    }

    CPLStringList aosOptions;
    aosOptions.AddString("-of");
    aosOptions.AddString("VRT");

    OGREnvelope env;
    poClipGeom->getEnvelope(&env);

    if (m_onlyBBOX)
    {
        auto poPoly = std::make_unique<OGRPolygon>(env);
        poPoly->assignSpatialReference(poClipGeom->getSpatialReference());
        poClipGeom = std::move(poPoly);
    }

    const bool bBottomUpRaster = gt.yscale > 0;

    if (poClipGeom->IsRectangle() && !m_addAlpha && !bBottomUpRaster)
    {
        aosOptions.AddString("-projwin");
        aosOptions.AddString(CPLSPrintf("%.17g", env.MinX));
        aosOptions.AddString(CPLSPrintf("%.17g", env.MaxY));
        aosOptions.AddString(CPLSPrintf("%.17g", env.MaxX));
        aosOptions.AddString(CPLSPrintf("%.17g", env.MinY));

        auto poClipGeomSRS = poClipGeom->getSpatialReference();
        if (poClipGeomSRS)
        {
            const char *const apszOptions[] = {"FORMAT=WKT2", nullptr};
            const std::string osWKT = poClipGeomSRS->exportToWkt(apszOptions);
            aosOptions.AddString("-projwin_srs");
            aosOptions.AddString(osWKT.c_str());
        }

        if (m_allowExtentOutsideSource)
        {
            aosOptions.AddString("--no-warn-about-outside-window");
        }
        else
        {
            // Unless we've specifically allowed the bounding box to extend beyond
            // the source raster, raise an error.
            aosOptions.AddString("-epo");
        }

        GDALTranslateOptions *psOptions =
            GDALTranslateOptionsNew(aosOptions.List(), nullptr);

        GDALDatasetH hSrcDS = GDALDataset::ToHandle(poSrcDS);
        auto poRetDS = GDALDataset::FromHandle(
            GDALTranslate("", hSrcDS, psOptions, nullptr));
        GDALTranslateOptionsFree(psOptions);

        const bool bOK = poRetDS != nullptr;
        if (bOK)
        {
            m_outputDataset.Set(std::unique_ptr<GDALDataset>(poRetDS));
        }

        return bOK;
    }
    else
    {
        if (bBottomUpRaster)
        {
            gt.yorig += gt.yscale * poSrcDS->GetRasterYSize();
            gt.yscale = -gt.yscale;
        }

        {
            auto poClipGeomInSrcSRS =
                std::unique_ptr<OGRGeometry>(poClipGeom->clone());
            if (poClipGeom->getSpatialReference() && poSrcDS->GetSpatialRef())
                poClipGeomInSrcSRS->transformTo(poSrcDS->GetSpatialRef());
            poClipGeomInSrcSRS->getEnvelope(&env);
        }

        OGREnvelope rasterEnv;
        poSrcDS->GetExtent(&rasterEnv, nullptr);
        if (!m_allowExtentOutsideSource && !rasterEnv.Contains(env))
        {
            ReportError(CE_Failure, CPLE_AppDefined,
                        "Clipping geometry is partially or totally outside the "
                        "extent of the raster. You can set the "
                        "'allow-bbox-outside-source' argument to proceed.");
            return false;
        }

        if (m_addAlpha)
        {
            aosOptions.AddString("-dstalpha");
        }

        aosOptions.AddString("-cutline");
        aosOptions.AddString(poClipGeom->exportToWkt());

        aosOptions.AddString("-wo");
        aosOptions.AddString("CUTLINE_ALL_TOUCHED=YES");

        auto poClipGeomSRS = poClipGeom->getSpatialReference();
        if (poClipGeomSRS)
        {
            const char *const apszOptions[] = {"FORMAT=WKT2", nullptr};
            const std::string osWKT = poClipGeomSRS->exportToWkt(apszOptions);
            aosOptions.AddString("-cutline_srs");
            aosOptions.AddString(osWKT.c_str());
        }

        constexpr double REL_EPS_PIXEL = 1e-3;
        const double dfMinX =
            gt.xorig +
            floor((env.MinX - gt.xorig) / gt.xscale + REL_EPS_PIXEL) *
                gt.xscale;
        const double dfMinY =
            gt.yorig +
            ceil((env.MinY - gt.yorig) / gt.yscale - REL_EPS_PIXEL) * gt.yscale;
        const double dfMaxX =
            gt.xorig +
            ceil((env.MaxX - gt.xorig) / gt.xscale - REL_EPS_PIXEL) * gt.xscale;
        const double dfMaxY =
            gt.yorig +
            floor((env.MaxY - gt.yorig) / gt.yscale + REL_EPS_PIXEL) *
                gt.yscale;

        aosOptions.AddString("-te");
        aosOptions.AddString(CPLSPrintf("%.17g", dfMinX));
        aosOptions.AddString(
            CPLSPrintf("%.17g", bBottomUpRaster ? dfMaxY : dfMinY));
        aosOptions.AddString(CPLSPrintf("%.17g", dfMaxX));
        aosOptions.AddString(
            CPLSPrintf("%.17g", bBottomUpRaster ? dfMinY : dfMaxY));

        aosOptions.AddString("-tr");
        aosOptions.AddString(CPLSPrintf("%.17g", gt.xscale));
        aosOptions.AddString(CPLSPrintf("%.17g", std::fabs(gt.yscale)));

        GDALWarpAppOptions *psOptions =
            GDALWarpAppOptionsNew(aosOptions.List(), nullptr);

        GDALDatasetH hSrcDS = GDALDataset::ToHandle(poSrcDS);
        auto poRetDS = GDALDataset::FromHandle(
            GDALWarp("", nullptr, 1, &hSrcDS, psOptions, nullptr));
        GDALWarpAppOptionsFree(psOptions);

        const bool bOK = poRetDS != nullptr;
        if (bOK)
        {
            m_outputDataset.Set(std::unique_ptr<GDALDataset>(poRetDS));
        }

        return bOK;
    }
}

GDALRasterClipAlgorithmStandalone::~GDALRasterClipAlgorithmStandalone() =
    default;

//! @endcond

Coverage Report

Created: 2026-07-14 06:21

Line	Count	Source
1		/******************************************************************************
2		*
3		* Project: GDAL
4		* Purpose: "clip" step of "raster pipeline", or "gdal raster clip" standalone
5		* Author: Even Rouault <even dot rouault at spatialys.com>
6		*
7		******************************************************************************
8		* Copyright (c) 2025, Even Rouault <even dot rouault at spatialys.com>
9		*
10		* SPDX-License-Identifier: MIT
11		****************************************************************************/
12
13		#include "gdalalg_raster_clip.h"
14
15		#include "gdal_priv.h"
16		#include "gdal_utils.h"
17
18		#include <algorithm>
19		#include <cmath>
20
21		//! @cond Doxygen_Suppress
22
23		#ifndef _
24	0	#define _(x) (x)
25		#endif
26
27		/************************************************************************/
28		/* GDALRasterClipAlgorithm::GDALRasterClipAlgorithm() */
29		/************************************************************************/
30
31		GDALRasterClipAlgorithm::GDALRasterClipAlgorithm(bool standaloneStep)
32	0	: GDALRasterPipelineStepAlgorithm(NAME, DESCRIPTION, HELP_URL,
33	0	standaloneStep)
34	0	{
35	0	constexpr const char *EXCLUSION_GROUP = "bbox-window-geometry-like";
36	0	AddBBOXArg(&m_bbox, _("Clipping bounding box as xmin,ymin,xmax,ymax"))
37	0	.SetMutualExclusionGroup(EXCLUSION_GROUP);
38	0	AddArg("bbox-crs", 0, _("CRS of clipping bounding box"), &m_bboxCrs)
39	0	.SetIsCRSArg()
40	0	.AddHiddenAlias("bbox_srs");
41
42	0	AddArg("window", 0, _("Raster window as col,line,width,height in pixels"),
43	0	&m_window)
44	0	.SetRepeatedArgAllowed(false)
45	0	.SetMinCount(4)
46	0	.SetMaxCount(4)
47	0	.SetDisplayHintAboutRepetition(false)
48	0	.SetMutualExclusionGroup(EXCLUSION_GROUP)
49	0	.AddValidationAction(
50	0	[this]()
51	0	{
52	0	CPLAssert(m_window.size() == 4);
53	0	if (m_window[2] <= 0 \|\| m_window[3] <= 0)
54	0	{
55	0	CPLError(CE_Failure, CPLE_AppDefined,
56	0	"Value of 'window' should be "
57	0	"col,line,width,height with "
58	0	"width > 0 and height > 0");
59	0	return false;
60	0	}
61	0	return true;
62	0	});
63
64	0	AddArg("geometry", 0, _("Clipping geometry (WKT or GeoJSON)"), &m_geometry)
65	0	.SetMutualExclusionGroup(EXCLUSION_GROUP);
66	0	AddArg("geometry-crs", 0, _("CRS of clipping geometry"), &m_geometryCrs)
67	0	.SetIsCRSArg()
68	0	.AddHiddenAlias("geometry_srs");
69	0	AddArg("like", 0, _("Dataset to use as a template for bounds"),
70	0	&m_likeDataset, GDAL_OF_RASTER \| GDAL_OF_VECTOR)
71	0	.SetMetaVar("DATASET")
72	0	.SetMutualExclusionGroup(EXCLUSION_GROUP);
73	0	AddArg("like-sql", 0, ("SELECT statement to run on the 'like' dataset"),
74	0	&m_likeSQL)
75	0	.SetMetaVar("SELECT-STATEMENT")
76	0	.SetMutualExclusionGroup("sql-where");
77	0	AddArg("like-layer", 0, ("Name of the layer of the 'like' dataset"),
78	0	&m_likeLayer)
79	0	.SetMetaVar("LAYER-NAME");
80	0	AddArg("like-where", 0, ("WHERE SQL clause to run on the 'like' dataset"),
81	0	&m_likeWhere)
82	0	.SetMetaVar("WHERE-EXPRESSION")
83	0	.SetMutualExclusionGroup("sql-where");
84	0	AddArg("only-bbox", 0,
85	0	_("For 'geometry' and 'like', only consider their bounding box"),
86	0	&m_onlyBBOX);
87	0	AddArg("allow-bbox-outside-source", 0,
88	0	_("Allow clipping box to include pixels outside input dataset"),
89	0	&m_allowExtentOutsideSource);
90		// Note: this would be a use case for multiple exclusion groups because it is not
91		// compatible with "window"
92	0	AddArg("add-alpha", 0,
93	0	_("Adds an alpha mask band to the destination when the source "
94	0	"raster have none."),
95	0	&m_addAlpha);
96	0	}
97
98		/************************************************************************/
99		/* GDALRasterClipAlgorithm::RunStep() */
100		/************************************************************************/
101
102		bool GDALRasterClipAlgorithm::RunStep(GDALPipelineStepRunContext &)
103	0	{
104	0	auto poSrcDS = m_inputDataset[0].GetDatasetRef();
105	0	CPLAssert(poSrcDS);
106	0	CPLAssert(m_outputDataset.GetName().empty());
107	0	CPLAssert(!m_outputDataset.GetDatasetRef());
108
109	0	if (!m_window.empty())
110	0	{
111	0	if (m_addAlpha)
112	0	{
113	0	ReportError(CE_Failure, CPLE_NotSupported,
114	0	"'alpha' argument is not supported with 'window'");
115	0	return false;
116	0	}
117
118	0	CPLStringList aosOptions;
119	0	aosOptions.AddString("-of");
120	0	aosOptions.AddString("VRT");
121
122	0	aosOptions.AddString("-srcwin");
123	0	aosOptions.AddString(CPLSPrintf("%d", m_window[0]));
124	0	aosOptions.AddString(CPLSPrintf("%d", m_window[1]));
125	0	aosOptions.AddString(CPLSPrintf("%d", m_window[2]));
126	0	aosOptions.AddString(CPLSPrintf("%d", m_window[3]));
127
128	0	if (!m_allowExtentOutsideSource)
129	0	{
130		// Unless we've specifically allowed the bounding box to extend beyond
131		// the source raster, raise an error.
132	0	aosOptions.AddString("-epo");
133	0	}
134
135	0	GDALTranslateOptions *psOptions =
136	0	GDALTranslateOptionsNew(aosOptions.List(), nullptr);
137
138	0	GDALDatasetH hSrcDS = GDALDataset::ToHandle(poSrcDS);
139	0	auto poRetDS = GDALDataset::FromHandle(
140	0	GDALTranslate("", hSrcDS, psOptions, nullptr));
141	0	GDALTranslateOptionsFree(psOptions);
142
143	0	const bool bOK = poRetDS != nullptr;
144	0	if (bOK)
145	0	{
146	0	m_outputDataset.Set(std::unique_ptr<GDALDataset>(poRetDS));
147	0	}
148
149	0	return bOK;
150	0	}
151
152	0	GDALGeoTransform gt;
153	0	if (poSrcDS->GetGeoTransform(gt) != CE_None)
154	0	{
155	0	ReportError(
156	0	CE_Failure, CPLE_NotSupported,
157	0	"Clipping is not supported on a raster without a geotransform");
158	0	return false;
159	0	}
160	0	if (!gt.IsAxisAligned())
161	0	{
162	0	ReportError(CE_Failure, CPLE_NotSupported,
163	0	"Clipping is not supported on a raster whose geotransform "
164	0	"has rotation terms");
165	0	return false;
166	0	}
167
168	0	auto [poClipGeom, errMsg] = GetClipGeometry();
169	0	if (!poClipGeom)
170	0	{
171	0	ReportError(CE_Failure, CPLE_AppDefined, "%s", errMsg.c_str());
172	0	return false;
173	0	}
174
175	0	auto poLikeDS = m_likeDataset.GetDatasetRef();
176	0	if (!poClipGeom->getSpatialReference() && poLikeDS &&
177	0	poLikeDS->GetLayerCount() == 0)
178	0	{
179	0	ReportError(CE_Failure, CPLE_AppDefined,
180	0	"Dataset '%s' has no CRS. Its bounds cannot be used.",
181	0	poLikeDS->GetDescription());
182	0	return false;
183	0	}
184
185	0	CPLStringList aosOptions;
186	0	aosOptions.AddString("-of");
187	0	aosOptions.AddString("VRT");
188
189	0	OGREnvelope env;
190	0	poClipGeom->getEnvelope(&env);
191
192	0	if (m_onlyBBOX)
193	0	{
194	0	auto poPoly = std::make_unique<OGRPolygon>(env);
195	0	poPoly->assignSpatialReference(poClipGeom->getSpatialReference());
196	0	poClipGeom = std::move(poPoly);
197	0	}
198
199	0	const bool bBottomUpRaster = gt.yscale > 0;
200
201	0	if (poClipGeom->IsRectangle() && !m_addAlpha && !bBottomUpRaster)
202	0	{
203	0	aosOptions.AddString("-projwin");
204	0	aosOptions.AddString(CPLSPrintf("%.17g", env.MinX));
205	0	aosOptions.AddString(CPLSPrintf("%.17g", env.MaxY));
206	0	aosOptions.AddString(CPLSPrintf("%.17g", env.MaxX));
207	0	aosOptions.AddString(CPLSPrintf("%.17g", env.MinY));
208
209	0	auto poClipGeomSRS = poClipGeom->getSpatialReference();
210	0	if (poClipGeomSRS)
211	0	{
212	0	const char *const apszOptions[] = {"FORMAT=WKT2", nullptr};
213	0	const std::string osWKT = poClipGeomSRS->exportToWkt(apszOptions);
214	0	aosOptions.AddString("-projwin_srs");
215	0	aosOptions.AddString(osWKT.c_str());
216	0	}
217
218	0	if (m_allowExtentOutsideSource)
219	0	{
220	0	aosOptions.AddString("--no-warn-about-outside-window");
221	0	}
222	0	else
223	0	{
224		// Unless we've specifically allowed the bounding box to extend beyond
225		// the source raster, raise an error.
226	0	aosOptions.AddString("-epo");
227	0	}
228
229	0	GDALTranslateOptions *psOptions =
230	0	GDALTranslateOptionsNew(aosOptions.List(), nullptr);
231
232	0	GDALDatasetH hSrcDS = GDALDataset::ToHandle(poSrcDS);
233	0	auto poRetDS = GDALDataset::FromHandle(
234	0	GDALTranslate("", hSrcDS, psOptions, nullptr));
235	0	GDALTranslateOptionsFree(psOptions);
236
237	0	const bool bOK = poRetDS != nullptr;
238	0	if (bOK)
239	0	{
240	0	m_outputDataset.Set(std::unique_ptr<GDALDataset>(poRetDS));
241	0	}
242
243	0	return bOK;
244	0	}
245	0	else
246	0	{
247	0	if (bBottomUpRaster)
248	0	{
249	0	gt.yorig += gt.yscale * poSrcDS->GetRasterYSize();
250	0	gt.yscale = -gt.yscale;
251	0	}
252
253	0	{
254	0	auto poClipGeomInSrcSRS =
255	0	std::unique_ptr<OGRGeometry>(poClipGeom->clone());
256	0	if (poClipGeom->getSpatialReference() && poSrcDS->GetSpatialRef())
257	0	poClipGeomInSrcSRS->transformTo(poSrcDS->GetSpatialRef());
258	0	poClipGeomInSrcSRS->getEnvelope(&env);
259	0	}
260
261	0	OGREnvelope rasterEnv;
262	0	poSrcDS->GetExtent(&rasterEnv, nullptr);
263	0	if (!m_allowExtentOutsideSource && !rasterEnv.Contains(env))
264	0	{
265	0	ReportError(CE_Failure, CPLE_AppDefined,
266	0	"Clipping geometry is partially or totally outside the "
267	0	"extent of the raster. You can set the "
268	0	"'allow-bbox-outside-source' argument to proceed.");
269	0	return false;
270	0	}
271
272	0	if (m_addAlpha)
273	0	{
274	0	aosOptions.AddString("-dstalpha");
275	0	}
276
277	0	aosOptions.AddString("-cutline");
278	0	aosOptions.AddString(poClipGeom->exportToWkt());
279
280	0	aosOptions.AddString("-wo");
281	0	aosOptions.AddString("CUTLINE_ALL_TOUCHED=YES");
282
283	0	auto poClipGeomSRS = poClipGeom->getSpatialReference();
284	0	if (poClipGeomSRS)
285	0	{
286	0	const char *const apszOptions[] = {"FORMAT=WKT2", nullptr};
287	0	const std::string osWKT = poClipGeomSRS->exportToWkt(apszOptions);
288	0	aosOptions.AddString("-cutline_srs");
289	0	aosOptions.AddString(osWKT.c_str());
290	0	}
291
292	0	constexpr double REL_EPS_PIXEL = 1e-3;
293	0	const double dfMinX =
294	0	gt.xorig +
295	0	floor((env.MinX - gt.xorig) / gt.xscale + REL_EPS_PIXEL) *
296	0	gt.xscale;
297	0	const double dfMinY =
298	0	gt.yorig +
299	0	ceil((env.MinY - gt.yorig) / gt.yscale - REL_EPS_PIXEL) * gt.yscale;
300	0	const double dfMaxX =
301	0	gt.xorig +
302	0	ceil((env.MaxX - gt.xorig) / gt.xscale - REL_EPS_PIXEL) * gt.xscale;
303	0	const double dfMaxY =
304	0	gt.yorig +
305	0	floor((env.MaxY - gt.yorig) / gt.yscale + REL_EPS_PIXEL) *
306	0	gt.yscale;
307
308	0	aosOptions.AddString("-te");
309	0	aosOptions.AddString(CPLSPrintf("%.17g", dfMinX));
310	0	aosOptions.AddString(
311	0	CPLSPrintf("%.17g", bBottomUpRaster ? dfMaxY : dfMinY));
312	0	aosOptions.AddString(CPLSPrintf("%.17g", dfMaxX));
313	0	aosOptions.AddString(
314	0	CPLSPrintf("%.17g", bBottomUpRaster ? dfMinY : dfMaxY));
315
316	0	aosOptions.AddString("-tr");
317	0	aosOptions.AddString(CPLSPrintf("%.17g", gt.xscale));
318	0	aosOptions.AddString(CPLSPrintf("%.17g", std::fabs(gt.yscale)));
319
320	0	GDALWarpAppOptions *psOptions =
321	0	GDALWarpAppOptionsNew(aosOptions.List(), nullptr);
322
323	0	GDALDatasetH hSrcDS = GDALDataset::ToHandle(poSrcDS);
324	0	auto poRetDS = GDALDataset::FromHandle(
325	0	GDALWarp("", nullptr, 1, &hSrcDS, psOptions, nullptr));
326	0	GDALWarpAppOptionsFree(psOptions);
327
328	0	const bool bOK = poRetDS != nullptr;
329	0	if (bOK)
330	0	{
331	0	m_outputDataset.Set(std::unique_ptr<GDALDataset>(poRetDS));
332	0	}
333
334	0	return bOK;
335	0	}
336	0	}
337
338	0	GDALRasterClipAlgorithmStandalone::~GDALRasterClipAlgorithmStandalone() =
339		default;
340
341		//! @endcond