/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);
    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-02-14 06:52

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	0	AddArg("add-alpha", 0,
91	0	_("Adds an alpha mask band to the destination when the source "
92	0	"raster have none."),
93	0	&m_addAlpha);
94	0	}
95
96		/************************************************************************/
97		/* GDALRasterClipAlgorithm::RunStep() */
98		/************************************************************************/
99
100		bool GDALRasterClipAlgorithm::RunStep(GDALPipelineStepRunContext &)
101	0	{
102	0	auto poSrcDS = m_inputDataset[0].GetDatasetRef();
103	0	CPLAssert(poSrcDS);
104	0	CPLAssert(m_outputDataset.GetName().empty());
105	0	CPLAssert(!m_outputDataset.GetDatasetRef());
106
107	0	if (!m_window.empty())
108	0	{
109	0	if (m_addAlpha)
110	0	{
111	0	ReportError(CE_Failure, CPLE_NotSupported,
112	0	"'alpha' argument is not supported with 'window'");
113	0	return false;
114	0	}
115
116	0	CPLStringList aosOptions;
117	0	aosOptions.AddString("-of");
118	0	aosOptions.AddString("VRT");
119
120	0	aosOptions.AddString("-srcwin");
121	0	aosOptions.AddString(CPLSPrintf("%d", m_window[0]));
122	0	aosOptions.AddString(CPLSPrintf("%d", m_window[1]));
123	0	aosOptions.AddString(CPLSPrintf("%d", m_window[2]));
124	0	aosOptions.AddString(CPLSPrintf("%d", m_window[3]));
125
126	0	if (!m_allowExtentOutsideSource)
127	0	{
128		// Unless we've specifically allowed the bounding box to extend beyond
129		// the source raster, raise an error.
130	0	aosOptions.AddString("-epo");
131	0	}
132
133	0	GDALTranslateOptions *psOptions =
134	0	GDALTranslateOptionsNew(aosOptions.List(), nullptr);
135
136	0	GDALDatasetH hSrcDS = GDALDataset::ToHandle(poSrcDS);
137	0	auto poRetDS = GDALDataset::FromHandle(
138	0	GDALTranslate("", hSrcDS, psOptions, nullptr));
139	0	GDALTranslateOptionsFree(psOptions);
140
141	0	const bool bOK = poRetDS != nullptr;
142	0	if (bOK)
143	0	{
144	0	m_outputDataset.Set(std::unique_ptr<GDALDataset>(poRetDS));
145	0	}
146
147	0	return bOK;
148	0	}
149
150	0	GDALGeoTransform gt;
151	0	if (poSrcDS->GetGeoTransform(gt) != CE_None)
152	0	{
153	0	ReportError(
154	0	CE_Failure, CPLE_NotSupported,
155	0	"Clipping is not supported on a raster without a geotransform");
156	0	return false;
157	0	}
158	0	if (!gt.IsAxisAligned())
159	0	{
160	0	ReportError(CE_Failure, CPLE_NotSupported,
161	0	"Clipping is not supported on a raster whose geotransform "
162	0	"has rotation terms");
163	0	return false;
164	0	}
165
166	0	auto [poClipGeom, errMsg] = GetClipGeometry();
167	0	if (!poClipGeom)
168	0	{
169	0	ReportError(CE_Failure, CPLE_AppDefined, "%s", errMsg.c_str());
170	0	return false;
171	0	}
172
173	0	auto poLikeDS = m_likeDataset.GetDatasetRef();
174	0	if (!poClipGeom->getSpatialReference() && poLikeDS &&
175	0	poLikeDS->GetLayerCount() == 0)
176	0	{
177	0	ReportError(CE_Failure, CPLE_AppDefined,
178	0	"Dataset '%s' has no CRS. Its bounds cannot be used.",
179	0	poLikeDS->GetDescription());
180	0	return false;
181	0	}
182
183	0	CPLStringList aosOptions;
184	0	aosOptions.AddString("-of");
185	0	aosOptions.AddString("VRT");
186
187	0	OGREnvelope env;
188	0	poClipGeom->getEnvelope(&env);
189
190	0	if (m_onlyBBOX)
191	0	{
192	0	auto poPoly = std::make_unique<OGRPolygon>(env);
193	0	poPoly->assignSpatialReference(poClipGeom->getSpatialReference());
194	0	poClipGeom = std::move(poPoly);
195	0	}
196
197	0	const bool bBottomUpRaster = gt.yscale > 0;
198
199	0	if (poClipGeom->IsRectangle() && !m_addAlpha && !bBottomUpRaster)
200	0	{
201	0	aosOptions.AddString("-projwin");
202	0	aosOptions.AddString(CPLSPrintf("%.17g", env.MinX));
203	0	aosOptions.AddString(CPLSPrintf("%.17g", env.MaxY));
204	0	aosOptions.AddString(CPLSPrintf("%.17g", env.MaxX));
205	0	aosOptions.AddString(CPLSPrintf("%.17g", env.MinY));
206
207	0	auto poClipGeomSRS = poClipGeom->getSpatialReference();
208	0	if (poClipGeomSRS)
209	0	{
210	0	const char *const apszOptions[] = {"FORMAT=WKT2", nullptr};
211	0	const std::string osWKT = poClipGeomSRS->exportToWkt(apszOptions);
212	0	aosOptions.AddString("-projwin_srs");
213	0	aosOptions.AddString(osWKT.c_str());
214	0	}
215
216	0	if (m_allowExtentOutsideSource)
217	0	{
218	0	aosOptions.AddString("--no-warn-about-outside-window");
219	0	}
220	0	else
221	0	{
222		// Unless we've specifically allowed the bounding box to extend beyond
223		// the source raster, raise an error.
224	0	aosOptions.AddString("-epo");
225	0	}
226
227	0	GDALTranslateOptions *psOptions =
228	0	GDALTranslateOptionsNew(aosOptions.List(), nullptr);
229
230	0	GDALDatasetH hSrcDS = GDALDataset::ToHandle(poSrcDS);
231	0	auto poRetDS = GDALDataset::FromHandle(
232	0	GDALTranslate("", hSrcDS, psOptions, nullptr));
233	0	GDALTranslateOptionsFree(psOptions);
234
235	0	const bool bOK = poRetDS != nullptr;
236	0	if (bOK)
237	0	{
238	0	m_outputDataset.Set(std::unique_ptr<GDALDataset>(poRetDS));
239	0	}
240
241	0	return bOK;
242	0	}
243	0	else
244	0	{
245	0	if (bBottomUpRaster)
246	0	{
247	0	gt.yorig += gt.yscale * poSrcDS->GetRasterYSize();
248	0	gt.yscale = -gt.yscale;
249	0	}
250
251	0	{
252	0	auto poClipGeomInSrcSRS =
253	0	std::unique_ptr<OGRGeometry>(poClipGeom->clone());
254	0	if (poClipGeom->getSpatialReference() && poSrcDS->GetSpatialRef())
255	0	poClipGeomInSrcSRS->transformTo(poSrcDS->GetSpatialRef());
256	0	poClipGeomInSrcSRS->getEnvelope(&env);
257	0	}
258
259	0	OGREnvelope rasterEnv;
260	0	poSrcDS->GetExtent(&rasterEnv, nullptr);
261	0	if (!m_allowExtentOutsideSource && !rasterEnv.Contains(env))
262	0	{
263	0	ReportError(CE_Failure, CPLE_AppDefined,
264	0	"Clipping geometry is partially or totally outside the "
265	0	"extent of the raster. You can set the "
266	0	"'allow-bbox-outside-source' argument to proceed.");
267	0	return false;
268	0	}
269
270	0	if (m_addAlpha)
271	0	{
272	0	aosOptions.AddString("-dstalpha");
273	0	}
274
275	0	aosOptions.AddString("-cutline");
276	0	aosOptions.AddString(poClipGeom->exportToWkt());
277
278	0	aosOptions.AddString("-wo");
279	0	aosOptions.AddString("CUTLINE_ALL_TOUCHED=YES");
280
281	0	auto poClipGeomSRS = poClipGeom->getSpatialReference();
282	0	if (poClipGeomSRS)
283	0	{
284	0	const char *const apszOptions[] = {"FORMAT=WKT2", nullptr};
285	0	const std::string osWKT = poClipGeomSRS->exportToWkt(apszOptions);
286	0	aosOptions.AddString("-cutline_srs");
287	0	aosOptions.AddString(osWKT.c_str());
288	0	}
289
290	0	constexpr double REL_EPS_PIXEL = 1e-3;
291	0	const double dfMinX =
292	0	gt.xorig +
293	0	floor((env.MinX - gt.xorig) / gt.xscale + REL_EPS_PIXEL) *
294	0	gt.xscale;
295	0	const double dfMinY =
296	0	gt.yorig +
297	0	ceil((env.MinY - gt.yorig) / gt.yscale - REL_EPS_PIXEL) * gt.yscale;
298	0	const double dfMaxX =
299	0	gt.xorig +
300	0	ceil((env.MaxX - gt.xorig) / gt.xscale - REL_EPS_PIXEL) * gt.xscale;
301	0	const double dfMaxY =
302	0	gt.yorig +
303	0	floor((env.MaxY - gt.yorig) / gt.yscale + REL_EPS_PIXEL) *
304	0	gt.yscale;
305
306	0	aosOptions.AddString("-te");
307	0	aosOptions.AddString(CPLSPrintf("%.17g", dfMinX));
308	0	aosOptions.AddString(
309	0	CPLSPrintf("%.17g", bBottomUpRaster ? dfMaxY : dfMinY));
310	0	aosOptions.AddString(CPLSPrintf("%.17g", dfMaxX));
311	0	aosOptions.AddString(
312	0	CPLSPrintf("%.17g", bBottomUpRaster ? dfMinY : dfMaxY));
313
314	0	aosOptions.AddString("-tr");
315	0	aosOptions.AddString(CPLSPrintf("%.17g", gt.xscale));
316	0	aosOptions.AddString(CPLSPrintf("%.17g", std::fabs(gt.yscale)));
317
318	0	GDALWarpAppOptions *psOptions =
319	0	GDALWarpAppOptionsNew(aosOptions.List(), nullptr);
320
321	0	GDALDatasetH hSrcDS = GDALDataset::ToHandle(poSrcDS);
322	0	auto poRetDS = GDALDataset::FromHandle(
323	0	GDALWarp("", nullptr, 1, &hSrcDS, psOptions, nullptr));
324	0	GDALWarpAppOptionsFree(psOptions);
325
326	0	const bool bOK = poRetDS != nullptr;
327	0	if (bOK)
328	0	{
329	0	m_outputDataset.Set(std::unique_ptr<GDALDataset>(poRetDS));
330	0	}
331
332	0	return bOK;
333	0	}
334	0	}
335
336	0	GDALRasterClipAlgorithmStandalone::~GDALRasterClipAlgorithmStandalone() =
337		default;
338
339		//! @endcond