diff --git a/apps/CMakeLists.txt b/apps/CMakeLists.txt index 47492b836e30..40af26a6f8d6 100644 --- a/apps/CMakeLists.txt +++ b/apps/CMakeLists.txt @@ -13,6 +13,7 @@ add_library( gdalalg_raster.cpp gdalalg_raster_buildvrt.cpp gdalalg_raster_info.cpp + gdalalg_raster_clip.cpp gdalalg_raster_convert.cpp gdalalg_raster_edit.cpp gdalalg_raster_pipeline.cpp diff --git a/apps/gdalalg_raster.cpp b/apps/gdalalg_raster.cpp index b2975197a303..bb0ad3fe8e7a 100644 --- a/apps/gdalalg_raster.cpp +++ b/apps/gdalalg_raster.cpp @@ -14,6 +14,7 @@ #include "gdalalg_raster_buildvrt.h" #include "gdalalg_raster_info.h" +#include "gdalalg_raster_clip.h" #include "gdalalg_raster_convert.h" #include "gdalalg_raster_edit.h" #include "gdalalg_raster_pipeline.h" @@ -39,6 +40,7 @@ class GDALRasterAlgorithm final : public GDALAlgorithm { RegisterSubAlgorithm(); RegisterSubAlgorithm(); + RegisterSubAlgorithm(); RegisterSubAlgorithm(); RegisterSubAlgorithm(); RegisterSubAlgorithm(); diff --git a/apps/gdalalg_raster_clip.cpp b/apps/gdalalg_raster_clip.cpp new file mode 100644 index 000000000000..cfa806749324 --- /dev/null +++ b/apps/gdalalg_raster_clip.cpp @@ -0,0 +1,145 @@ +/****************************************************************************** + * + * Project: GDAL + * Purpose: "clip" step of "raster pipeline", or "gdal raster clip" standalone + * Author: Even Rouault + * + ****************************************************************************** + * Copyright (c) 2025, Even Rouault + * + * SPDX-License-Identifier: MIT + ****************************************************************************/ + +#include "gdalalg_raster_clip.h" + +#include "gdal_priv.h" +#include "gdal_utils.h" + +#include + +//! @cond Doxygen_Suppress + +#ifndef _ +#define _(x) (x) +#endif + +/************************************************************************/ +/* GDALRasterClipAlgorithm::GDALRasterClipAlgorithm() */ +/************************************************************************/ + +GDALRasterClipAlgorithm::GDALRasterClipAlgorithm(bool standaloneStep) + : GDALRasterPipelineStepAlgorithm(NAME, DESCRIPTION, HELP_URL, + standaloneStep) +{ + AddBBOXArg(&m_bbox, _("Clipping bounding box as xmin,ymin,xmax,ymax")) + .SetMutualExclusionGroup("exclusion-group"); + AddArg("bbox-crs", 0, _("CRS of clipping bounding box"), &m_bboxCrs) + .SetIsCRSArg() + .AddHiddenAlias("bbox_srs"); + AddArg("like", 0, _("Raster dataset to use as a template for bounds"), + &m_likeDataset, GDAL_OF_RASTER) + .SetMetaVar("DATASET") + .SetMutualExclusionGroup("exclusion-group"); +} + +/************************************************************************/ +/* GDALRasterClipAlgorithm::RunStep() */ +/************************************************************************/ + +bool GDALRasterClipAlgorithm::RunStep(GDALProgressFunc, void *) +{ + CPLAssert(m_inputDataset.GetDatasetRef()); + CPLAssert(m_outputDataset.GetName().empty()); + CPLAssert(!m_outputDataset.GetDatasetRef()); + + CPLStringList aosOptions; + aosOptions.AddString("-of"); + aosOptions.AddString("VRT"); + if (!m_bbox.empty()) + { + aosOptions.AddString("-projwin"); + aosOptions.AddString(CPLSPrintf("%.17g", m_bbox[0])); // minx + aosOptions.AddString(CPLSPrintf("%.17g", m_bbox[3])); // maxy + aosOptions.AddString(CPLSPrintf("%.17g", m_bbox[2])); // maxx + aosOptions.AddString(CPLSPrintf("%.17g", m_bbox[1])); // miny + + if (!m_bboxCrs.empty()) + { + aosOptions.AddString("-projwin_srs"); + aosOptions.AddString(m_bboxCrs.c_str()); + } + } + else if (auto poLikeDS = m_likeDataset.GetDatasetRef()) + { + 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 (!poLikeSRS) + { + ReportError(CE_Failure, CPLE_AppDefined, + "Dataset '%s' has no SRS. Its bounds cannot be used.", + poLikeDS->GetDescription()); + return false; + } + 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]; + const double dfMinX = + std::min(std::min(dfTLX, dfTRX), std::min(dfBLX, dfBRX)); + const double dfMinY = + std::min(std::min(dfTLY, dfTRY), std::min(dfBLY, dfBRY)); + const double dfMaxX = + std::max(std::max(dfTLX, dfTRX), std::max(dfBLX, dfBRX)); + const double dfMaxY = + std::max(std::max(dfTLY, dfTRY), std::max(dfBLY, dfBRY)); + + aosOptions.AddString("-projwin"); + aosOptions.AddString(CPLSPrintf("%.17g", dfMinX)); + aosOptions.AddString(CPLSPrintf("%.17g", dfMaxY)); + aosOptions.AddString(CPLSPrintf("%.17g", dfMaxX)); + aosOptions.AddString(CPLSPrintf("%.17g", dfMinY)); + + const char *const apszOptions[] = {"FORMAT=WKT2", nullptr}; + const std::string osWKT = poLikeSRS->exportToWkt(apszOptions); + aosOptions.AddString("-projwin_srs"); + aosOptions.AddString(osWKT.c_str()); + } + else + { + ReportError(CE_Failure, CPLE_AppDefined, + "Either --bbox or --like must be specified"); + return false; + } + + GDALTranslateOptions *psOptions = + GDALTranslateOptionsNew(aosOptions.List(), nullptr); + + GDALDatasetH hSrcDS = GDALDataset::ToHandle(m_inputDataset.GetDatasetRef()); + auto poRetDS = + GDALDataset::FromHandle(GDALTranslate("", hSrcDS, psOptions, nullptr)); + GDALTranslateOptionsFree(psOptions); + + const bool bOK = poRetDS != nullptr; + if (bOK) + { + m_outputDataset.Set(std::unique_ptr(poRetDS)); + } + + return bOK; +} + +//! @endcond diff --git a/apps/gdalalg_raster_clip.h b/apps/gdalalg_raster_clip.h new file mode 100644 index 000000000000..7ac795bd9890 --- /dev/null +++ b/apps/gdalalg_raster_clip.h @@ -0,0 +1,62 @@ +/****************************************************************************** + * + * Project: GDAL + * Purpose: "clip" step of "raster pipeline", or "gdal raster clip" standalone + * Author: Even Rouault + * + ****************************************************************************** + * Copyright (c) 2025, Even Rouault + * + * SPDX-License-Identifier: MIT + ****************************************************************************/ + +#ifndef GDALALG_RASTER_CLIP_INCLUDED +#define GDALALG_RASTER_CLIP_INCLUDED + +#include "gdalalg_raster_pipeline.h" + +//! @cond Doxygen_Suppress + +/************************************************************************/ +/* GDALRasterClipAlgorithm */ +/************************************************************************/ + +class GDALRasterClipAlgorithm /* non final */ + : public GDALRasterPipelineStepAlgorithm +{ + public: + static constexpr const char *NAME = "clip"; + static constexpr const char *DESCRIPTION = "Clip a raster dataset."; + static constexpr const char *HELP_URL = "/programs/gdal_raster_clip.html"; + + static std::vector GetAliases() + { + return {}; + } + + explicit GDALRasterClipAlgorithm(bool standaloneStep = false); + + private: + bool RunStep(GDALProgressFunc pfnProgress, void *pProgressData) override; + + std::vector m_bbox{}; + std::string m_bboxCrs{}; + GDALArgDatasetValue m_likeDataset{}; +}; + +/************************************************************************/ +/* GDALRasterClipAlgorithmStandalone */ +/************************************************************************/ + +class GDALRasterClipAlgorithmStandalone final : public GDALRasterClipAlgorithm +{ + public: + GDALRasterClipAlgorithmStandalone() + : GDALRasterClipAlgorithm(/* standaloneStep = */ true) + { + } +}; + +//! @endcond + +#endif /* GDALALG_RASTER_CLIP_INCLUDED */ diff --git a/apps/gdalalg_raster_pipeline.cpp b/apps/gdalalg_raster_pipeline.cpp index 0ab163fa531c..45b4214f7f45 100644 --- a/apps/gdalalg_raster_pipeline.cpp +++ b/apps/gdalalg_raster_pipeline.cpp @@ -12,6 +12,7 @@ #include "gdalalg_raster_pipeline.h" #include "gdalalg_raster_read.h" +#include "gdalalg_raster_clip.h" #include "gdalalg_raster_edit.h" #include "gdalalg_raster_reproject.h" #include "gdalalg_raster_write.h" @@ -161,6 +162,7 @@ GDALRasterPipelineAlgorithm::GDALRasterPipelineAlgorithm( m_stepRegistry.Register(); m_stepRegistry.Register(); + m_stepRegistry.Register(); m_stepRegistry.Register(); m_stepRegistry.Register(); } diff --git a/autotest/utilities/test_gdalalg_raster_pipeline.py b/autotest/utilities/test_gdalalg_raster_pipeline.py index e75519dc9dfe..55ede9b0fa23 100755 --- a/autotest/utilities/test_gdalalg_raster_pipeline.py +++ b/autotest/utilities/test_gdalalg_raster_pipeline.py @@ -539,3 +539,156 @@ def test_gdalalg_raster_pipeline_reproject_almost_all_args(tmp_vsimem): (-117.641, 0.0005, 0.0, 33.9008, 0.0, -0.0004), rel=1e-8 ) assert ds.GetRasterBand(1).Checksum() == 8515 + + +def test_gdalalg_raster_pipeline_clip_missing_bbox_or_like(tmp_vsimem): + + out_filename = str(tmp_vsimem / "out.tif") + + pipeline = get_pipeline_alg() + with pytest.raises( + Exception, match="clip: Either --bbox or --like must be specified" + ): + pipeline.ParseRunAndFinalize( + [ + "read", + "../gcore/data/byte.tif", + "!", + "clip", + "!", + "write", + "--overwrite", + out_filename, + ] + ) + + +def test_gdalalg_raster_pipeline_clip(tmp_vsimem): + + out_filename = str(tmp_vsimem / "out.tif") + + pipeline = get_pipeline_alg() + assert pipeline.ParseRunAndFinalize( + [ + "read", + "../gcore/data/byte.tif", + "!", + "clip", + "--bbox=440780,3750200,441860,3751260", + "!", + "write", + "--overwrite", + out_filename, + ] + ) + + with gdal.Open(out_filename) as ds: + assert ds.RasterXSize == 18 + assert ds.RasterYSize == 18 + assert ds.GetSpatialRef().GetAuthorityCode(None) == "26711" + assert ds.GetGeoTransform() == pytest.approx( + (440780, 60, 0, 3751260, 0, -60), rel=1e-8 + ) + assert ds.GetRasterBand(1).Checksum() == 3695 + + out2_filename = str(tmp_vsimem / "out2.tif") + + pipeline = get_pipeline_alg() + assert pipeline.ParseRunAndFinalize( + [ + "read", + "../gcore/data/byte.tif", + "!", + "clip", + "--like", + out_filename, + "!", + "write", + "--overwrite", + out2_filename, + ] + ) + + with gdal.Open(out_filename) as ds: + assert ds.RasterXSize == 18 + assert ds.RasterYSize == 18 + assert ds.GetSpatialRef().GetAuthorityCode(None) == "26711" + assert ds.GetGeoTransform() == pytest.approx( + (440780, 60, 0, 3751260, 0, -60), rel=1e-8 + ) + assert ds.GetRasterBand(1).Checksum() == 3695 + + +def test_gdalalg_raster_pipeline_clip_like_error(tmp_vsimem): + + ref_filename = str(tmp_vsimem / "out.tif") + gdal.GetDriverByName("GTiff").Create(ref_filename, 1, 1) + + out_filename = str(tmp_vsimem / "out.tif") + + pipeline = get_pipeline_alg() + with pytest.raises(Exception, match="has no geotransform matrix"): + pipeline.ParseRunAndFinalize( + [ + "read", + "../gcore/data/byte.tif", + "!", + "clip", + "--like", + ref_filename, + "!", + "write", + "--overwrite", + out_filename, + ] + ) + + with gdal.GetDriverByName("GTiff").Create(ref_filename, 1, 1) as ds: + ds.SetGeoTransform([2, 1, 0, 49, 0, -1]) + + pipeline = get_pipeline_alg() + with pytest.raises(Exception, match="has no SRS"): + pipeline.ParseRunAndFinalize( + [ + "read", + "../gcore/data/byte.tif", + "!", + "clip", + "--like", + ref_filename, + "!", + "write", + "--overwrite", + out_filename, + ] + ) + + +def test_gdalalg_raster_pipeline_clip_bbox_crs(tmp_vsimem): + + out_filename = str(tmp_vsimem / "out.tif") + + pipeline = get_pipeline_alg() + assert pipeline.ParseRunAndFinalize( + [ + "read", + "../gcore/data/byte.tif", + "!", + "clip", + "--bbox=-117.631,33.89,-117.628,33.9005", + "--bbox-crs=NAD27", + "!", + "write", + "--overwrite", + out_filename, + ] + ) + + with gdal.Open(out_filename) as ds: + assert ds.RasterXSize == 5 + assert ds.RasterYSize == 19 + assert ds.GetSpatialRef().GetAuthorityCode(None) == "26711" + print(ds.GetGeoTransform()) + assert ds.GetGeoTransform() == pytest.approx( + (441620.0, 60.0, 0.0, 3751140.0, 0.0, -60.0), rel=1e-8 + ) diff --git a/doc/source/conf.py b/doc/source/conf.py index 625663f20718..f2fc81f3a83d 100644 --- a/doc/source/conf.py +++ b/doc/source/conf.py @@ -229,6 +229,13 @@ [author_evenr], 1, ), + ( + "programs/gdal_raster_clip", + "gdal-raster-clip", + "Clip a raster dataset", + [author_evenr], + 1, + ), ( "programs/gdal_raster_convert", "gdal-raster-convert", diff --git a/doc/source/programs/gdal_raster.rst b/doc/source/programs/gdal_raster.rst index 67cd7bd12fe1..25d4dfcae415 100644 --- a/doc/source/programs/gdal_raster.rst +++ b/doc/source/programs/gdal_raster.rst @@ -20,6 +20,7 @@ Synopsis Usage: gdal raster where is one of: - buildvrt: Build a virtual dataset (VRT). + - clip: Clip a raster dataset. - convert: Convert a raster dataset. - info: Return information on a raster dataset. - pipeline: Process a raster dataset. @@ -30,6 +31,7 @@ Available sub-commands - :ref:`gdal_raster_buildvrt_subcommand` - :ref:`gdal_raster_info_subcommand` +- :ref:`gdal_raster_clip_subcommand` - :ref:`gdal_raster_convert_subcommand` - :ref:`gdal_raster_pipeline_subcommand` - :ref:`gdal_raster_reproject_subcommand` diff --git a/doc/source/programs/gdal_raster_clip.rst b/doc/source/programs/gdal_raster_clip.rst new file mode 100644 index 000000000000..eda664b5bb74 --- /dev/null +++ b/doc/source/programs/gdal_raster_clip.rst @@ -0,0 +1,117 @@ +.. _gdal_raster_clip_subcommand: + +================================================================================ +"gdal raster clip" sub-command +================================================================================ + +.. versionadded:: 3.11 + +.. only:: html + + Clip a raster dataset. + +.. Index:: gdal raster clip + +Synopsis +-------- + +.. code-block:: + + Usage: gdal raster clip [OPTIONS] + + Clip a raster dataset. + + Positional arguments: + -i, --input Input raster dataset [required] + -o, --output Output raster 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: + -f, --of, --format, --output-format Output format + --co, --creation-option = Creation option [may be repeated] + --overwrite Whether overwriting existing output is allowed + --bbox Clipping bounding box as xmin,ymin,xmax,ymax + Mutually exclusive with --like + --bbox-crs CRS of clipping bounding box + --like Raster dataset to use as a template for bounds + Mutually exclusive with --bbox + + Advanced Options: + --if, --input-format Input formats [may be repeated] + --oo, --open-option Open options [may be repeated] + + +Description +----------- + +:program:`gdal raster clip` can be used to clip a raster dataset using +georeferenced coordinates. + +Either :option:`--bbox` or :option:`--like` must be specified. + +The output dataset is in the same SRS as the input one, and the original +resolution is preserved. Bounds are rounded to match whole pixel locations +(i.e. there is no resampling involved) + +``clip`` can also be used as a step of :ref:`gdal_raster_pipeline_subcommand`. + +Standard options +++++++++++++++++ + +.. include:: gdal_options/of_raster_create_copy.rst + +.. include:: gdal_options/co.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. + +.. 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 raster reprojection. + Instead, the bounds are reprojected from the bbox-crs to the CRS of the + input dataset. + +.. option:: --like + + Raster dataset to use as a template for bounds. This option is mutually + exclusive with :option:`--bbox` and :option:`--bbox-crs`. + +Advanced options +++++++++++++++++ + +.. include:: gdal_options/oo.rst + +.. include:: gdal_options/if.rst + +Examples +-------- + +.. example:: + :title: Clip a GeoTIFF file to the bounding box from longitude 2, latitude 49, to longitude 3, lattidue 50 in WGS 84 + + .. code-block:: bash + + $ gdal raster clip --bbox=2,49,3,50 --bbox-crs=EPSG:4326 in.tif out.tif --overwrite + +.. example:: + :title: Clip a GeoTIFF file using the bounds of :file:`reference.tif` + + .. code-block:: bash + + $ gdal raster clip --like=reference.tif in.tif out.tif --overwrite diff --git a/doc/source/programs/gdal_raster_pipeline.rst b/doc/source/programs/gdal_raster_pipeline.rst index 62dff6e37ce4..d788909d0e81 100644 --- a/doc/source/programs/gdal_raster_pipeline.rst +++ b/doc/source/programs/gdal_raster_pipeline.rst @@ -49,6 +49,21 @@ 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 raster dataset. + + Options: + --bbox Clipping bounding box as xmin,ymin,xmax,ymax + Mutually exclusive with --like + --bbox-crs CRS of clipping bounding box + --like Raster dataset to use as a template for bounds + Mutually exclusive with --bbox + +Details for options can be found in :ref:`gdal_raster_clip_subcommand`. + * edit [OPTIONS] .. code-block:: diff --git a/doc/source/programs/index.rst b/doc/source/programs/index.rst index c2160d4b7ff5..5434f9b9751b 100644 --- a/doc/source/programs/index.rst +++ b/doc/source/programs/index.rst @@ -31,6 +31,7 @@ single :program:`gdal` program that accepts commands and subcommands. gdal_raster gdal_raster_buildvrt gdal_raster_info + gdal_raster_clip gdal_raster_convert gdal_raster_edit gdal_raster_pipeline @@ -48,6 +49,7 @@ single :program:`gdal` program that accepts commands and subcommands. - :ref:`gdal_raster_command`: Entry point for raster commands - :ref:`gdal_raster_buildvrt_subcommand`: buildvrt: Build a virtual dataset (VRT) - :ref:`gdal_raster_info_subcommand`: Get information on a raster dataset + - :ref:`gdal_raster_clip_subcommand`: Clip a raster dataset - :ref:`gdal_raster_convert_subcommand`: Convert a raster dataset - :ref:`gdal_raster_edit_subcommand`: Edit in place a raster dataset - :ref:`gdal_raster_pipeline_subcommand`: Process a raster dataset