Line data    Source code

       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  GDAL
       4             :  * Purpose:  gdal "raster create" subcommand
       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_create.h"
      14             : 
      15             : #include "cpl_conv.h"
      16             : #include "gdal_priv.h"
      17             : #include "gdal_utils.h"
      18             : #include "ogr_spatialref.h"
      19             : 
      20             : //! @cond Doxygen_Suppress
      21             : 
      22             : #ifndef _
      23             : #define _(x) (x)
      24             : #endif
      25             : 
      26             : /************************************************************************/
      27             : /*        GDALRasterCreateAlgorithm::GDALRasterCreateAlgorithm()        */
      28             : /************************************************************************/
      29             : 
      30         147 : GDALRasterCreateAlgorithm::GDALRasterCreateAlgorithm(
      31         147 :     bool standaloneStep) noexcept
      32             :     : GDALRasterPipelineStepAlgorithm(
      33             :           NAME, DESCRIPTION, HELP_URL,
      34         441 :           ConstructorOptions()
      35         147 :               .SetStandaloneStep(standaloneStep)
      36         147 :               .SetAddDefaultArguments(false)
      37         147 :               .SetAutoOpenInputDatasets(true)
      38         294 :               .SetInputDatasetHelpMsg("Template raster dataset")
      39         294 :               .SetInputDatasetAlias("like")
      40         294 :               .SetInputDatasetMetaVar("TEMPLATE-DATASET")
      41         147 :               .SetInputDatasetRequired(false)
      42         147 :               .SetInputDatasetPositional(false)
      43         441 :               .SetInputDatasetMaxCount(1))
      44             : {
      45         147 :     AddRasterInputArgs(false, false);
      46         147 :     if (standaloneStep)
      47             :     {
      48         107 :         AddProgressArg();
      49         107 :         AddRasterOutputArgs(false);
      50             :     }
      51             : 
      52         294 :     AddArg("size", 0, _("Output size in pixels"), &m_size)
      53         147 :         .SetMinCount(2)
      54         147 :         .SetMaxCount(2)
      55         147 :         .SetMinValueIncluded(0)
      56         147 :         .SetRepeatedArgAllowed(false)
      57         147 :         .SetDisplayHintAboutRepetition(false)
      58         147 :         .SetMetaVar("<width>,<height>");
      59         294 :     AddArg("band-count", 0, _("Number of bands"), &m_bandCount)
      60         147 :         .SetDefault(m_bandCount)
      61         147 :         .SetMinValueIncluded(0);
      62         147 :     AddOutputDataTypeArg(&m_type).SetDefault(m_type);
      63             : 
      64         147 :     AddNodataArg(&m_nodata, /* noneAllowed = */ true);
      65             : 
      66         147 :     AddArg("burn", 0, _("Burn value"), &m_burnValues);
      67         294 :     AddArg("crs", 0, _("Set CRS"), &m_crs)
      68         294 :         .AddHiddenAlias("a_srs")
      69         147 :         .SetIsCRSArg(/*noneAllowed=*/true);
      70         147 :     AddBBOXArg(&m_bbox);
      71             : 
      72             :     {
      73         294 :         auto &arg = AddArg("metadata", 0, _("Add metadata item"), &m_metadata)
      74         294 :                         .SetMetaVar("<KEY>=<VALUE>")
      75         147 :                         .SetPackedValuesAllowed(false);
      76          24 :         arg.AddValidationAction([this, &arg]()
      77         171 :                                 { return ParseAndValidateKeyValue(arg); });
      78         147 :         arg.AddHiddenAlias("mo");
      79             :     }
      80             : 
      81         147 :     const auto inputArg = GetArg(GDAL_ARG_NAME_INPUT);
      82         147 :     CPLAssertNotNull(inputArg);
      83             : 
      84             :     AddArg("copy-metadata", 0, _("Copy metadata from input dataset"),
      85         294 :            &m_copyMetadata)
      86         147 :         .AddDirectDependency(*inputArg);
      87             :     AddArg("copy-overviews", 0,
      88         294 :            _("Create same overview levels as input dataset"), &m_copyOverviews)
      89         147 :         .AddDirectDependency(*inputArg);
      90         147 : }
      91             : 
      92             : /************************************************************************/
      93             : /*                 GDALRasterCreateAlgorithm::RunImpl()                 */
      94             : /************************************************************************/
      95             : 
      96          55 : bool GDALRasterCreateAlgorithm::RunImpl(GDALProgressFunc pfnProgress,
      97             :                                         void *pProgressData)
      98             : {
      99          55 :     GDALPipelineStepRunContext stepCtxt;
     100          55 :     stepCtxt.m_pfnProgress = pfnProgress;
     101          55 :     stepCtxt.m_pProgressData = pProgressData;
     102          55 :     return RunPreStepPipelineValidations() && RunStep(stepCtxt);
     103             : }
     104             : 
     105             : /************************************************************************/
     106             : /*                 GDALRasterCreateAlgorithm::RunStep()                 */
     107             : /************************************************************************/
     108             : 
     109          58 : bool GDALRasterCreateAlgorithm::RunStep(GDALPipelineStepRunContext &)
     110             : {
     111          58 :     CPLAssert(!m_outputDataset.GetDatasetRef());
     112             : 
     113          58 :     if (m_standaloneStep)
     114             :     {
     115          55 :         if (m_format.empty())
     116             :         {
     117             :             const auto aosFormats =
     118             :                 CPLStringList(GDALGetOutputDriversForDatasetName(
     119          15 :                     m_outputDataset.GetName().c_str(), GDAL_OF_RASTER,
     120             :                     /* bSingleMatch = */ true,
     121          15 :                     /* bWarn = */ true));
     122          15 :             if (aosFormats.size() != 1)
     123             :             {
     124           1 :                 ReportError(CE_Failure, CPLE_AppDefined,
     125             :                             "Cannot guess driver for %s",
     126           1 :                             m_outputDataset.GetName().c_str());
     127           1 :                 return false;
     128             :             }
     129          14 :             m_format = aosFormats[0];
     130             :         }
     131             :     }
     132             :     else
     133             :     {
     134           3 :         m_format = "MEM";
     135             :     }
     136             : 
     137         114 :     OGRSpatialReference oSRS;
     138             : 
     139          57 :     GDALGeoTransform gt;
     140          57 :     bool bGTValid = false;
     141             : 
     142         114 :     CPLStringList aosCreationOptions(m_creationOptions);
     143             : 
     144          57 :     GDALDataset *poSrcDS = m_inputDataset.empty()
     145          57 :                                ? nullptr
     146          20 :                                : m_inputDataset.front().GetDatasetRef();
     147          57 :     if (poSrcDS)
     148             :     {
     149          20 :         if (m_size.empty())
     150             :         {
     151          54 :             m_size = std::vector<int>{poSrcDS->GetRasterXSize(),
     152          36 :                                       poSrcDS->GetRasterYSize()};
     153             :         }
     154             : 
     155          20 :         if (!GetArg("band-count")->IsExplicitlySet())
     156             :         {
     157          19 :             m_bandCount = poSrcDS->GetRasterCount();
     158             :         }
     159             : 
     160          20 :         if (!GetArg("datatype")->IsExplicitlySet())
     161             :         {
     162          19 :             if (m_bandCount > 0)
     163             :             {
     164             :                 m_type = GDALGetDataTypeName(
     165          19 :                     poSrcDS->GetRasterBand(1)->GetRasterDataType());
     166             :             }
     167             :         }
     168             : 
     169          20 :         if (m_crs.empty())
     170             :         {
     171          18 :             if (const auto poSRS = poSrcDS->GetSpatialRef())
     172          18 :                 oSRS = *poSRS;
     173             :         }
     174             : 
     175          20 :         if (m_bbox.empty())
     176             :         {
     177          19 :             bGTValid = poSrcDS->GetGeoTransform(gt) == CE_None;
     178             :         }
     179             : 
     180          20 :         if (m_nodata.empty() && m_bandCount > 0)
     181             :         {
     182          19 :             int bNoData = false;
     183             :             const double dfNoData =
     184          19 :                 poSrcDS->GetRasterBand(1)->GetNoDataValue(&bNoData);
     185          19 :             if (bNoData)
     186          10 :                 m_nodata = CPLSPrintf("%.17g", dfNoData);
     187             :         }
     188             : 
     189             :         // Replicate tiling of input datasets for a few popular output formats,
     190             :         // when compatible, and when the user hasn't specified creation options
     191             :         // affecting tiling.
     192          20 :         int nBlockXSize = 0, nBlockYSize = 0;
     193          20 :         if (m_bandCount > 0)
     194          20 :             poSrcDS->GetRasterBand(1)->GetBlockSize(&nBlockXSize, &nBlockYSize);
     195             : 
     196          20 :         if (EQUAL(m_format.c_str(), "GTIFF") &&
     197           3 :             aosCreationOptions.FetchNameValue("TILED") == nullptr &&
     198           2 :             aosCreationOptions.FetchNameValue("BLOCKXSIZE") == nullptr &&
     199          25 :             aosCreationOptions.FetchNameValue("BLOCKYSIZE") == nullptr &&
     200           2 :             m_bandCount > 0)
     201             :         {
     202           2 :             if (nBlockXSize != poSrcDS->GetRasterXSize() &&
     203           2 :                 (nBlockXSize % 16) == 0 && (nBlockYSize % 16) == 0)
     204             :             {
     205           1 :                 aosCreationOptions.SetNameValue("TILED", "YES");
     206             :                 aosCreationOptions.SetNameValue("BLOCKXSIZE",
     207           1 :                                                 CPLSPrintf("%d", nBlockXSize));
     208             :                 aosCreationOptions.SetNameValue("BLOCKYSIZE",
     209           1 :                                                 CPLSPrintf("%d", nBlockYSize));
     210             :             }
     211             :         }
     212          18 :         else if (EQUAL(m_format.c_str(), "COG") &&
     213          19 :                  aosCreationOptions.FetchNameValue("BLOCKSIZE") == nullptr &&
     214           1 :                  m_bandCount > 0)
     215             :         {
     216           1 :             if (nBlockXSize != poSrcDS->GetRasterXSize() &&
     217           2 :                 nBlockXSize == nBlockYSize && nBlockXSize >= 128 &&
     218           1 :                 (nBlockXSize % 16) == 0)
     219             :             {
     220             :                 aosCreationOptions.SetNameValue("BLOCKSIZE",
     221           1 :                                                 CPLSPrintf("%d", nBlockXSize));
     222             :             }
     223             :         }
     224          17 :         else if (EQUAL(m_format.c_str(), "GPKG") &&
     225           3 :                  aosCreationOptions.FetchNameValue("BLOCKSIZE") == nullptr &&
     226           2 :                  aosCreationOptions.FetchNameValue("BLOCKXSIZE") == nullptr &&
     227          22 :                  aosCreationOptions.FetchNameValue("BLOCKYSIZE") == nullptr &&
     228           2 :                  m_bandCount > 0)
     229             :         {
     230           2 :             if (nBlockXSize != poSrcDS->GetRasterXSize() &&
     231           2 :                 nBlockXSize >= 256 && nBlockXSize <= 4096 &&
     232           4 :                 nBlockYSize >= 256 && nBlockYSize <= 4096)
     233             :             {
     234             :                 aosCreationOptions.SetNameValue("BLOCKXSIZE",
     235           1 :                                                 CPLSPrintf("%d", nBlockXSize));
     236             :                 aosCreationOptions.SetNameValue("BLOCKYSIZE",
     237           1 :                                                 CPLSPrintf("%d", nBlockYSize));
     238             :             }
     239             :         }
     240             :     }
     241             : 
     242          57 :     if (m_size.empty())
     243             :     {
     244           1 :         ReportError(CE_Failure, CPLE_IllegalArg,
     245             :                     "Argument 'size' should be specified, or 'like' dataset "
     246             :                     "should be specified");
     247           1 :         return false;
     248             :     }
     249             : 
     250             :     // Guess the size from bbox if only one of the two is specified
     251         112 :     if (m_size.size() == 2 && (m_size[0] == 0 || m_size[1] == 0) &&
     252           3 :         !(m_size[0] == 0 && m_size[1] == 0) && m_bbox.size() == 4 &&
     253         112 :         (m_bbox[3] - m_bbox[1] != 0) && (m_bbox[2] - m_bbox[0] != 0))
     254             :     {
     255           2 :         constexpr double EPSILON = 1e-5;
     256           2 :         const double ratio = (m_bbox[2] - m_bbox[0]) / (m_bbox[3] - m_bbox[1]);
     257           2 :         if (m_size[0] == 0)
     258             :         {
     259           1 :             m_size[0] =
     260           1 :                 static_cast<int>(std::ceil(m_size[1] * ratio - EPSILON));
     261             :         }
     262           1 :         else if (m_size[1] == 0)
     263             :         {
     264           1 :             m_size[1] =
     265           1 :                 static_cast<int>(std::ceil(m_size[0] / ratio - EPSILON));
     266             :         }
     267             :     }
     268             : 
     269          70 :     if (!m_burnValues.empty() && m_burnValues.size() != 1 &&
     270          14 :         static_cast<int>(m_burnValues.size()) != m_bandCount)
     271             :     {
     272           2 :         if (m_bandCount == 1)
     273             :         {
     274           1 :             ReportError(CE_Failure, CPLE_IllegalArg,
     275             :                         "One value should be provided for argument "
     276             :                         "'burn', given there is one band");
     277             :         }
     278             :         else
     279             :         {
     280           1 :             ReportError(CE_Failure, CPLE_IllegalArg,
     281             :                         "One or %d values should be provided for argument "
     282             :                         "'burn', given there are %d bands",
     283             :                         m_bandCount, m_bandCount);
     284             :         }
     285           2 :         return false;
     286             :     }
     287             : 
     288          54 :     auto poDriver = GetGDALDriverManager()->GetDriverByName(m_format.c_str());
     289          54 :     if (!poDriver)
     290             :     {
     291             :         // shouldn't happen given checks done in GDALAlgorithm
     292           0 :         ReportError(CE_Failure, CPLE_AppDefined, "Cannot find driver %s",
     293             :                     m_format.c_str());
     294           0 :         return false;
     295             :     }
     296             : 
     297          54 :     if (m_appendRaster)
     298             :     {
     299           2 :         if (poDriver->GetMetadataItem(GDAL_DCAP_CREATE_SUBDATASETS) == nullptr)
     300             :         {
     301           1 :             ReportError(CE_Failure, CPLE_NotSupported,
     302             :                         "-append option not supported for driver %s",
     303           1 :                         poDriver->GetDescription());
     304           1 :             return false;
     305             :         }
     306           1 :         aosCreationOptions.SetNameValue("APPEND_SUBDATASET", "YES");
     307             :     }
     308             : 
     309             :     auto poRetDS = std::unique_ptr<GDALDataset>(poDriver->Create(
     310         159 :         m_outputDataset.GetName().c_str(), m_size[0], m_size[1], m_bandCount,
     311         106 :         GDALGetDataTypeByName(m_type.c_str()), aosCreationOptions.List()));
     312          53 :     if (!poRetDS)
     313             :     {
     314           1 :         return false;
     315             :     }
     316             : 
     317          52 :     if (!m_crs.empty() && m_crs != "none" && m_crs != "null")
     318             :     {
     319          17 :         oSRS.SetFromUserInput(m_crs.c_str());
     320          17 :         oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
     321             :     }
     322             : 
     323          52 :     if (!oSRS.IsEmpty())
     324             :     {
     325          35 :         if (poRetDS->SetSpatialRef(&oSRS) != CE_None)
     326             :         {
     327           1 :             ReportError(CE_Failure, CPLE_AppDefined, "Setting CRS failed");
     328           1 :             return false;
     329             :         }
     330             :     }
     331             : 
     332          51 :     if (!m_bbox.empty())
     333             :     {
     334          23 :         if (poRetDS->GetRasterXSize() == 0 || poRetDS->GetRasterYSize() == 0)
     335             :         {
     336           1 :             ReportError(CE_Failure, CPLE_AppDefined,
     337             :                         "Cannot set extent because one of dataset height or "
     338             :                         "width is null");
     339           1 :             return false;
     340             :         }
     341          22 :         bGTValid = true;
     342          22 :         gt.xorig = m_bbox[0];
     343          22 :         gt.xscale = (m_bbox[2] - m_bbox[0]) / poRetDS->GetRasterXSize();
     344          22 :         gt.xrot = 0;
     345          22 :         gt.yorig = m_bbox[3];
     346          22 :         gt.yrot = 0;
     347          22 :         gt.yscale = -(m_bbox[3] - m_bbox[1]) / poRetDS->GetRasterYSize();
     348             :     }
     349          50 :     if (bGTValid)
     350             :     {
     351          41 :         if (poRetDS->SetGeoTransform(gt) != CE_None)
     352             :         {
     353           1 :             ReportError(CE_Failure, CPLE_AppDefined, "Setting extent failed");
     354           1 :             return false;
     355             :         }
     356             :     }
     357             : 
     358          49 :     if (!m_nodata.empty() && !EQUAL(m_nodata.c_str(), "none"))
     359             :     {
     360          68 :         for (int i = 0; i < poRetDS->GetRasterCount(); ++i)
     361             :         {
     362          45 :             bool bCannotBeExactlyRepresented = false;
     363          45 :             if (poRetDS->GetRasterBand(i + 1)->SetNoDataValueAsString(
     364          45 :                     m_nodata.c_str(), &bCannotBeExactlyRepresented) != CE_None)
     365             :             {
     366           1 :                 if (bCannotBeExactlyRepresented)
     367             :                 {
     368           1 :                     ReportError(CE_Failure, CPLE_AppDefined,
     369             :                                 "Setting nodata value failed as it cannot be "
     370             :                                 "represented on its data type");
     371             :                 }
     372             :                 else
     373             :                 {
     374           0 :                     ReportError(CE_Failure, CPLE_AppDefined,
     375             :                                 "Setting nodata value failed");
     376             :                 }
     377           1 :                 return false;
     378             :             }
     379             :         }
     380             :     }
     381             : 
     382          48 :     if (m_copyMetadata)
     383             :     {
     384             : 
     385             :         // This should never happen because of the dependency set
     386           1 :         CPLAssertNotNull(poSrcDS);
     387             : 
     388             :         {
     389           1 :             const CPLStringList aosDomains(poSrcDS->GetMetadataDomainList());
     390           4 :             for (const char *domain : aosDomains)
     391             :             {
     392           3 :                 if (!EQUAL(domain, GDAL_MDD_IMAGE_STRUCTURE))
     393             :                 {
     394           4 :                     if (poRetDS->SetMetadata(poSrcDS->GetMetadata(domain),
     395           4 :                                              domain) != CE_None)
     396             :                     {
     397           0 :                         ReportError(CE_Failure, CPLE_AppDefined,
     398             :                                     "Cannot copy '%s' metadata domain", domain);
     399           0 :                         return false;
     400             :                     }
     401             :                 }
     402             :             }
     403             :         }
     404           3 :         for (int i = 0; i < m_bandCount; ++i)
     405             :         {
     406             :             const CPLStringList aosDomains(
     407           2 :                 poSrcDS->GetRasterBand(i + 1)->GetMetadataDomainList());
     408           3 :             for (const char *domain : aosDomains)
     409             :             {
     410           1 :                 if (!EQUAL(domain, GDAL_MDD_IMAGE_STRUCTURE))
     411             :                 {
     412           2 :                     if (poRetDS->GetRasterBand(i + 1)->SetMetadata(
     413           1 :                             poSrcDS->GetRasterBand(i + 1)->GetMetadata(domain),
     414           2 :                             domain) != CE_None)
     415             :                     {
     416           0 :                         ReportError(
     417             :                             CE_Failure, CPLE_AppDefined,
     418             :                             "Cannot copy '%s' metadata domain for band %d",
     419             :                             domain, i + 1);
     420           0 :                         return false;
     421             :                     }
     422             :                 }
     423             :             }
     424             :         }
     425             :     }
     426             : 
     427          96 :     const CPLStringList aosMD(m_metadata);
     428          60 :     for (const auto &[key, value] : cpl::IterateNameValue(aosMD))
     429             :     {
     430          12 :         if (poRetDS->SetMetadataItem(key, value) != CE_None)
     431             :         {
     432           0 :             ReportError(CE_Failure, CPLE_AppDefined,
     433             :                         "SetMetadataItem('%s', '%s') failed", key, value);
     434           0 :             return false;
     435             :         }
     436             :     }
     437             : 
     438          48 :     if (m_copyOverviews && m_bandCount > 0)
     439             :     {
     440             :         // This should never happen because of the dependency set
     441           2 :         CPLAssertNotNull(poSrcDS);
     442             : 
     443           3 :         if (poSrcDS->GetRasterXSize() != poRetDS->GetRasterXSize() ||
     444           1 :             poSrcDS->GetRasterYSize() != poRetDS->GetRasterYSize())
     445             :         {
     446           1 :             ReportError(CE_Failure, CPLE_AppDefined,
     447             :                         "Argument 'copy-overviews' can only be set when the "
     448             :                         "input and output datasets have the same dimension");
     449           1 :             return false;
     450             :         }
     451             :         const int nOverviewCount =
     452           1 :             poSrcDS->GetRasterBand(1)->GetOverviewCount();
     453           1 :         std::vector<int> anLevels;
     454           2 :         for (int i = 0; i < nOverviewCount; ++i)
     455             :         {
     456           1 :             const auto poOvrBand = poSrcDS->GetRasterBand(1)->GetOverview(i);
     457           1 :             const int nOvrFactor = GDALComputeOvFactor(
     458             :                 poOvrBand->GetXSize(), poSrcDS->GetRasterXSize(),
     459           1 :                 poOvrBand->GetYSize(), poSrcDS->GetRasterYSize());
     460           1 :             anLevels.push_back(nOvrFactor);
     461             :         }
     462           2 :         if (poRetDS->BuildOverviews(
     463           1 :                 "NONE", nOverviewCount, anLevels.data(),
     464             :                 /* nListBands = */ 0, /* panBandList = */ nullptr,
     465             :                 /* pfnProgress = */ nullptr, /* pProgressData = */ nullptr,
     466           1 :                 /* options = */ nullptr) != CE_None)
     467             :         {
     468           0 :             ReportError(CE_Failure, CPLE_AppDefined,
     469             :                         "Creating overview(s) failed");
     470           0 :             return false;
     471             :         }
     472             :     }
     473             : 
     474          47 :     if (!m_burnValues.empty())
     475             :     {
     476          42 :         for (int i = 0; i < m_bandCount; ++i)
     477             :         {
     478          27 :             const int burnValueIdx = m_burnValues.size() == 1 ? 0 : i;
     479          27 :             const auto poDstBand = poRetDS->GetRasterBand(i + 1);
     480             :             // cppcheck-suppress negativeContainerIndex
     481          27 :             if (poDstBand->Fill(m_burnValues[burnValueIdx]) != CE_None)
     482             :             {
     483           0 :                 ReportError(CE_Failure, CPLE_AppDefined,
     484             :                             "Setting burn value failed");
     485           0 :                 return false;
     486             :             }
     487             :         }
     488          15 :         if (poRetDS->FlushCache(false) != CE_None)
     489             :         {
     490           0 :             ReportError(CE_Failure, CPLE_AppDefined,
     491             :                         "Setting burn value failed");
     492           0 :             return false;
     493             :         }
     494             :     }
     495             : 
     496          47 :     m_outputDataset.Set(std::move(poRetDS));
     497             : 
     498          47 :     return true;
     499             : }
     500             : 
     501             : GDALRasterCreateAlgorithmStandalone::~GDALRasterCreateAlgorithmStandalone() =
     502             :     default;
     503             : 
     504             : //! @endcond