/****************************************************************************** * * Project: MapServer * Purpose: Commandline App to build tile index for raster files. * Author: Frank Warmerdam, warmerdam@pobox.com * ****************************************************************************** * Copyright (c) 2001, Frank Warmerdam, DM Solutions Group Inc * Copyright (c) 2007-2023, Even Rouault * * SPDX-License-Identifier: MIT ****************************************************************************/ #include "cpl_port.h" #include "cpl_conv.h" #include "cpl_minixml.h" #include "cpl_string.h" #include "gdal_utils.h" #include "gdal_priv.h" #include "gdal_utils_priv.h" #include "ogr_api.h" #include "ogrsf_frmts.h" #include "ogr_spatialref.h" #include "commonutils.h" #include "gdalargumentparser.h" #include #include #include #include #include typedef enum { FORMAT_AUTO, FORMAT_WKT, FORMAT_EPSG, FORMAT_PROJ } SrcSRSFormat; /************************************************************************/ /* GDALTileIndexRasterMetadata */ /************************************************************************/ struct GDALTileIndexRasterMetadata { OGRFieldType eType = OFTString; std::string osFieldName{}; std::string osRasterItemName{}; }; /************************************************************************/ /* GDALTileIndexOptions */ /************************************************************************/ struct GDALTileIndexOptions { bool bInvokedFromGdalRasterIndex = false; bool bOverwrite = false; bool bSkipErrors = false; std::string osFormat{}; std::string osIndexLayerName{}; std::string osLocationField = "location"; CPLStringList aosLCO{}; std::string osTargetSRS{}; bool bWriteAbsolutePath = false; bool bSkipDifferentProjection = false; std::string osSrcSRSFieldName{}; SrcSRSFormat eSrcSRSFormat = FORMAT_AUTO; double xres = std::numeric_limits::quiet_NaN(); double yres = std::numeric_limits::quiet_NaN(); double xmin = std::numeric_limits::quiet_NaN(); double ymin = std::numeric_limits::quiet_NaN(); double xmax = std::numeric_limits::quiet_NaN(); double ymax = std::numeric_limits::quiet_NaN(); std::string osBandCount{}; std::string osNodata{}; std::string osColorInterp{}; std::string osDataType{}; bool bMaskBand = false; std::vector aosMetadata{}; std::string osGTIFilename{}; bool bRecursive = false; double dfMinPixelSize = std::numeric_limits::quiet_NaN(); double dfMaxPixelSize = std::numeric_limits::quiet_NaN(); std::vector aoFetchMD{}; std::set oSetFilenameFilters{}; GDALProgressFunc pfnProgress = nullptr; void *pProgressData = nullptr; }; /************************************************************************/ /* GDALTileIndexAppOptionsGetParser() */ /************************************************************************/ static std::unique_ptr GDALTileIndexAppOptionsGetParser( GDALTileIndexOptions *psOptions, GDALTileIndexOptionsForBinary *psOptionsForBinary) { auto argParser = std::make_unique( "gdaltindex", /* bForBinary=*/psOptionsForBinary != nullptr); argParser->add_description( _("Build a tile index from a list of datasets.")); argParser->add_epilog( _("For more details, see the full documentation for gdaltindex at\n" "https://gdal.org/programs/gdaltindex.html")); // Hidden as used by gdal raster index argParser->add_argument("--invoked-from-gdal-raster-index") .store_into(psOptions->bInvokedFromGdalRasterIndex) .hidden(); // Hidden as used by gdal raster index argParser->add_argument("-skip_errors") .store_into(psOptions->bSkipErrors) .hidden(); argParser->add_argument("-overwrite") .flag() .store_into(psOptions->bOverwrite) .help(_("Overwrite the output tile index file if it already exists.")); argParser->add_argument("-recursive") .flag() .store_into(psOptions->bRecursive) .help(_("Whether directories specified in should be " "explored recursively.")); argParser->add_argument("-filename_filter") .metavar("") .append() .store_into(psOptions->oSetFilenameFilters) .help(_("Pattern that the filenames contained in directories pointed " "by should follow.")); argParser->add_argument("-min_pixel_size") .metavar("") .store_into(psOptions->dfMinPixelSize) .help(_("Minimum pixel size in term of geospatial extent per pixel " "(resolution) that a raster should have to be selected.")); argParser->add_argument("-max_pixel_size") .metavar("") .store_into(psOptions->dfMaxPixelSize) .help(_("Maximum pixel size in term of geospatial extent per pixel " "(resolution) that a raster should have to be selected.")); argParser->add_output_format_argument(psOptions->osFormat); argParser->add_argument("-tileindex") .metavar("") .store_into(psOptions->osLocationField) .help(_("Name of the layer in the tile index file.")); argParser->add_argument("-write_absolute_path") .flag() .store_into(psOptions->bWriteAbsolutePath) .help(_("Write the absolute path of the raster files in the tile index " "file.")); argParser->add_argument("-skip_different_projection") .flag() .store_into(psOptions->bSkipDifferentProjection) .help(_( "Only files with the same projection as files already inserted in " "the tile index will be inserted (unless -t_srs is specified).")); argParser->add_argument("-t_srs") .metavar("") .store_into(psOptions->osTargetSRS) .help(_("Geometries of input files will be transformed to the desired " "target coordinate reference system.")); argParser->add_argument("-src_srs_name") .metavar("") .store_into(psOptions->osSrcSRSFieldName) .help(_("Name of the field in the tile index file where the source SRS " "will be stored.")); argParser->add_argument("-src_srs_format") .metavar("{AUTO|WKT|EPSG|PROJ}") .choices("AUTO", "WKT", "EPSG", "PROJ") .action( [psOptions](const auto &f) { if (f == "WKT") psOptions->eSrcSRSFormat = FORMAT_WKT; else if (f == "EPSG") psOptions->eSrcSRSFormat = FORMAT_EPSG; else if (f == "PROJ") psOptions->eSrcSRSFormat = FORMAT_PROJ; else psOptions->eSrcSRSFormat = FORMAT_AUTO; }) .help(_("Format of the source SRS to store in the tile index file.")); argParser->add_argument("-lyr_name") .metavar("") .store_into(psOptions->osIndexLayerName) .help(_("Name of the layer in the tile index file.")); argParser->add_layer_creation_options_argument(psOptions->aosLCO); // GTI driver options argParser->add_argument("-gti_filename") .metavar("") .store_into(psOptions->osGTIFilename) .help(_("Filename of the XML Virtual Tile Index file to generate.")); // NOTE: no store_into argParser->add_argument("-tr") .metavar(" ") .nargs(2) .scan<'g', double>() .help(_("Set target resolution.")); // NOTE: no store_into argParser->add_argument("-te") .metavar(" ") .nargs(4) .scan<'g', double>() .help(_("Set target extent in SRS unit.")); argParser->add_argument("-ot") .metavar("") .store_into(psOptions->osDataType) .help(_("Output data type.")); argParser->add_argument("-bandcount") .metavar("") .store_into(psOptions->osBandCount) .help(_("Number of bands of the tiles of the tile index.")); argParser->add_argument("-nodata") .metavar("") .append() .store_into(psOptions->osNodata) .help(_("Nodata value of the tiles of the tile index.")); // Should we use choices here? argParser->add_argument("-colorinterp") .metavar("") .append() .store_into(psOptions->osColorInterp) .help(_("Color interpretation of of the tiles of the tile index: red, " "green, blue, alpha, gray, undefined.")); argParser->add_argument("-mask") .flag() .store_into(psOptions->bMaskBand) .help(_("Add a mask band to the tiles of the tile index.")); argParser->add_argument("-mo") .metavar("=") .append() .store_into(psOptions->aosMetadata) .help(_("Write an arbitrary layer metadata item, for formats that " "support layer metadata.")); // NOTE: no store_into argParser->add_argument("-fetch_md") .nargs(3) .metavar(" ") .append() .help("Fetch a metadata item from the raster tile and write it as a " "field in the tile index."); if (psOptionsForBinary) { argParser->add_quiet_argument(&psOptionsForBinary->bQuiet); argParser->add_argument("index_file") .metavar("") .store_into(psOptionsForBinary->osDest) .help(_("The name of the output file to create/append to.")); argParser->add_argument("file_or_dir") .metavar("") .nargs(argparse::nargs_pattern::at_least_one) .action([psOptionsForBinary](const std::string &s) { psOptionsForBinary->aosSrcFiles.AddString(s.c_str()); }) .help(_( "The input GDAL raster files or directory, can be multiple " "locations separated by spaces. Wildcards may also be used.")); } return argParser; } /************************************************************************/ /* GDALTileIndexAppGetParserUsage() */ /************************************************************************/ std::string GDALTileIndexAppGetParserUsage() { try { GDALTileIndexOptions sOptions; GDALTileIndexOptionsForBinary sOptionsForBinary; auto argParser = GDALTileIndexAppOptionsGetParser(&sOptions, &sOptionsForBinary); return argParser->usage(); } catch (const std::exception &err) { CPLError(CE_Failure, CPLE_AppDefined, "Unexpected exception: %s", err.what()); return std::string(); } } /************************************************************************/ /* GDALTileIndexTileIterator */ /************************************************************************/ struct GDALTileIndexTileIterator { const GDALTileIndexOptions *psOptions = nullptr; int nSrcCount = 0; const char *const *papszSrcDSNames = nullptr; std::string osCurDir{}; int iCurSrc = 0; VSIDIR *psDir = nullptr; CPL_DISALLOW_COPY_ASSIGN(GDALTileIndexTileIterator) GDALTileIndexTileIterator(const GDALTileIndexOptions *psOptionsIn, int nSrcCountIn, const char *const *papszSrcDSNamesIn) : psOptions(psOptionsIn), nSrcCount(nSrcCountIn), papszSrcDSNames(papszSrcDSNamesIn) { } void reset() { if (psDir) VSICloseDir(psDir); psDir = nullptr; iCurSrc = 0; } std::string next() { while (true) { if (!psDir) { if (iCurSrc == nSrcCount) { break; } VSIStatBufL sStatBuf; const char *pszCurName = papszSrcDSNames[iCurSrc++]; if (VSIStatL(pszCurName, &sStatBuf) == 0 && VSI_ISDIR(sStatBuf.st_mode)) { auto poSrcDS = std::unique_ptr( GDALDataset::Open(pszCurName, GDAL_OF_RASTER, nullptr, nullptr, nullptr)); if (poSrcDS) return pszCurName; osCurDir = pszCurName; psDir = VSIOpenDir( osCurDir.c_str(), /*nDepth=*/psOptions->bRecursive ? -1 : 0, nullptr); if (!psDir) { CPLError(CE_Failure, CPLE_AppDefined, "Cannot open directory %s", osCurDir.c_str()); return std::string(); } } else { return pszCurName; } } auto psEntry = VSIGetNextDirEntry(psDir); if (!psEntry) { VSICloseDir(psDir); psDir = nullptr; continue; } if (!psOptions->oSetFilenameFilters.empty()) { bool bMatchFound = false; const std::string osFilenameOnly = CPLGetFilename(psEntry->pszName); for (const auto &osFilter : psOptions->oSetFilenameFilters) { if (GDALPatternMatch(osFilenameOnly.c_str(), osFilter.c_str())) { bMatchFound = true; break; } } if (!bMatchFound) continue; } std::string osFilename = CPLFormFilenameSafe( osCurDir.c_str(), psEntry->pszName, nullptr); if (VSI_ISDIR(psEntry->nMode)) { auto poSrcDS = std::unique_ptr( GDALDataset::Open(osFilename.c_str(), GDAL_OF_RASTER, nullptr, nullptr, nullptr)); if (poSrcDS) { return osFilename; } continue; } return osFilename; } return std::string(); } }; /************************************************************************/ /* GDALTileIndex() */ /************************************************************************/ /* clang-format off */ /** * Build a tile index from a list of datasets. * * This is the equivalent of the * gdaltindex utility. * * GDALTileIndexOptions* must be allocated and freed with * GDALTileIndexOptionsNew() and GDALTileIndexOptionsFree() respectively. * * @param pszDest the destination dataset path. * @param nSrcCount the number of input datasets. * @param papszSrcDSNames the list of input dataset names * @param psOptionsIn the options struct returned by GDALTileIndexOptionsNew() or * NULL. * @param pbUsageError pointer to a integer output variable to store if any * usage error has occurred. * @return the output dataset (new dataset that must be closed using * GDALClose()) or NULL in case of error. * * @since GDAL3.9 */ /* clang-format on */ GDALDatasetH GDALTileIndex(const char *pszDest, int nSrcCount, const char *const *papszSrcDSNames, const GDALTileIndexOptions *psOptionsIn, int *pbUsageError) { return GDALTileIndexInternal(pszDest, nullptr, nullptr, nSrcCount, papszSrcDSNames, psOptionsIn, pbUsageError); } GDALDatasetH GDALTileIndexInternal(const char *pszDest, GDALDatasetH hTileIndexDS, OGRLayerH hLayer, int nSrcCount, const char *const *papszSrcDSNames, const GDALTileIndexOptions *psOptionsIn, int *pbUsageError) { if (nSrcCount == 0) { CPLError(CE_Failure, CPLE_AppDefined, "No input dataset specified."); if (pbUsageError) *pbUsageError = TRUE; return nullptr; } auto psOptions = psOptionsIn ? std::make_unique(*psOptionsIn) : std::make_unique(); GDALTileIndexTileIterator oGDALTileIndexTileIterator( psOptions.get(), nSrcCount, papszSrcDSNames); /* -------------------------------------------------------------------- */ /* Create and validate target SRS if given. */ /* -------------------------------------------------------------------- */ OGRSpatialReference oTargetSRS; if (!psOptions->osTargetSRS.empty()) { if (psOptions->bSkipDifferentProjection) { CPLError(CE_Warning, CPLE_AppDefined, "-skip_different_projections does not apply " "when -t_srs is requested."); } oTargetSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER); // coverity[tainted_data] oTargetSRS.SetFromUserInput(psOptions->osTargetSRS.c_str()); } /* -------------------------------------------------------------------- */ /* Open or create the target datasource */ /* -------------------------------------------------------------------- */ std::unique_ptr poTileIndexDSUnique; GDALDataset *poTileIndexDS = GDALDataset::FromHandle(hTileIndexDS); OGRLayer *poLayer = OGRLayer::FromHandle(hLayer); bool bExistingLayer = false; std::string osFormat; if (!hTileIndexDS) { if (psOptions->bOverwrite) { CPLPushErrorHandler(CPLQuietErrorHandler); auto hDriver = GDALIdentifyDriver(pszDest, nullptr); if (hDriver) GDALDeleteDataset(hDriver, pszDest); else VSIUnlink(pszDest); CPLPopErrorHandler(); } poTileIndexDSUnique.reset( GDALDataset::Open(pszDest, GDAL_OF_VECTOR | GDAL_OF_UPDATE, nullptr, nullptr, nullptr)); if (poTileIndexDSUnique != nullptr) { auto poDriver = poTileIndexDSUnique->GetDriver(); if (poDriver) osFormat = poDriver->GetDescription(); if (poTileIndexDSUnique->GetLayerCount() == 1) { poLayer = poTileIndexDSUnique->GetLayer(0); } else { if (psOptions->osIndexLayerName.empty()) { CPLError(CE_Failure, CPLE_AppDefined, "Multiple layers detected: -lyr_name must be " "specified."); if (pbUsageError) *pbUsageError = true; return nullptr; } CPLPushErrorHandler(CPLQuietErrorHandler); poLayer = poTileIndexDSUnique->GetLayerByName( psOptions->osIndexLayerName.c_str()); CPLPopErrorHandler(); } } else { if (psOptions->osFormat.empty()) { const auto aoDrivers = GetOutputDriversFor(pszDest, GDAL_OF_VECTOR); if (aoDrivers.empty()) { CPLError(CE_Failure, CPLE_AppDefined, "Cannot guess driver for %s", pszDest); return nullptr; } else { if (aoDrivers.size() > 1) { CPLError( CE_Warning, CPLE_AppDefined, "Several drivers matching %s extension. Using %s", CPLGetExtensionSafe(pszDest).c_str(), aoDrivers[0].c_str()); } osFormat = aoDrivers[0]; } } else { osFormat = psOptions->osFormat; } auto poDriver = GetGDALDriverManager()->GetDriverByName(osFormat.c_str()); if (poDriver == nullptr) { CPLError(CE_Warning, CPLE_AppDefined, "%s driver not available.", osFormat.c_str()); return nullptr; } poTileIndexDSUnique.reset( poDriver->Create(pszDest, 0, 0, 0, GDT_Unknown, nullptr)); if (!poTileIndexDSUnique) return nullptr; } poTileIndexDS = poTileIndexDSUnique.get(); } auto poOutDrv = poTileIndexDS->GetDriver(); if (osFormat.empty() && poOutDrv) osFormat = poOutDrv->GetDescription(); const char *pszVal = poOutDrv ? poOutDrv->GetMetadataItem(GDAL_DMD_MAX_STRING_LENGTH) : nullptr; const int nMaxFieldSize = pszVal ? atoi(pszVal) : 0; const bool bFailOnErrors = psOptions->bInvokedFromGdalRasterIndex && !psOptions->bSkipErrors; bool bSkipFirstTile = false; if (poLayer) { bExistingLayer = true; } else { std::string osLayerName; if (psOptions->osIndexLayerName.empty()) { VSIStatBuf sStat; if (EQUAL(osFormat.c_str(), "ESRI Shapefile") || VSIStat(pszDest, &sStat) == 0) { osLayerName = CPLGetBasenameSafe(pszDest); } else { CPLError(CE_Failure, CPLE_AppDefined, "-lyr_name must be specified."); if (pbUsageError) *pbUsageError = true; return nullptr; } } else { if (psOptions->bOverwrite) { for (int i = 0; i < poTileIndexDS->GetLayerCount(); ++i) { auto poExistingLayer = poTileIndexDS->GetLayer(i); if (poExistingLayer && poExistingLayer->GetName() == psOptions->osIndexLayerName) { if (poTileIndexDS->DeleteLayer(i) != OGRERR_NONE) return nullptr; break; } } } osLayerName = psOptions->osIndexLayerName; } /* get spatial reference for output file from target SRS (if set) */ /* or from first input file */ OGRSpatialReference oSRS; if (!oTargetSRS.IsEmpty()) { oSRS = oTargetSRS; } else { std::string osFilename = oGDALTileIndexTileIterator.next(); if (osFilename.empty()) { CPLError(CE_Failure, CPLE_AppDefined, "Cannot find any tile"); return nullptr; } oGDALTileIndexTileIterator.reset(); std::unique_ptr poFailureIntoWarning; if (!bFailOnErrors) poFailureIntoWarning = std::make_unique(); CPL_IGNORE_RET_VAL(poFailureIntoWarning); auto poSrcDS = std::unique_ptr(GDALDataset::Open( osFilename.c_str(), GDAL_OF_RASTER | GDAL_OF_VERBOSE_ERROR, nullptr, nullptr, nullptr)); if (!poSrcDS) { CPLError(bFailOnErrors ? CE_Failure : CE_Warning, CPLE_AppDefined, "Unable to open %s%s.", osFilename.c_str(), bFailOnErrors ? "" : ", skipping"); if (bFailOnErrors) return nullptr; bSkipFirstTile = true; } else { auto poSrcSRS = poSrcDS->GetSpatialRef(); if (poSrcSRS) oSRS = *poSrcSRS; } } poLayer = poTileIndexDS->CreateLayer( osLayerName.c_str(), oSRS.IsEmpty() ? nullptr : &oSRS, wkbPolygon, psOptions->aosLCO.List()); if (!poLayer) return nullptr; OGRFieldDefn oLocationField(psOptions->osLocationField.c_str(), OFTString); oLocationField.SetWidth(nMaxFieldSize); if (poLayer->CreateField(&oLocationField) != OGRERR_NONE) return nullptr; if (!psOptions->osSrcSRSFieldName.empty()) { OGRFieldDefn oSrcSRSField(psOptions->osSrcSRSFieldName.c_str(), OFTString); oSrcSRSField.SetWidth(nMaxFieldSize); if (poLayer->CreateField(&oSrcSRSField) != OGRERR_NONE) return nullptr; } } auto poLayerDefn = poLayer->GetLayerDefn(); for (const auto &oFetchMD : psOptions->aoFetchMD) { if (poLayerDefn->GetFieldIndex(oFetchMD.osFieldName.c_str()) < 0) { OGRFieldDefn oField(oFetchMD.osFieldName.c_str(), oFetchMD.eType); if (poLayer->CreateField(&oField) != OGRERR_NONE) return nullptr; } } if (!psOptions->osGTIFilename.empty()) { if (!psOptions->aosMetadata.empty()) { CPLError(CE_Failure, CPLE_NotSupported, "-mo is not supported when -gti_filename is used"); return nullptr; } CPLXMLNode *psRoot = CPLCreateXMLNode(nullptr, CXT_Element, "GDALTileIndexDataset"); CPLCreateXMLElementAndValue(psRoot, "IndexDataset", pszDest); CPLCreateXMLElementAndValue(psRoot, "IndexLayer", poLayer->GetName()); CPLCreateXMLElementAndValue(psRoot, "LocationField", psOptions->osLocationField.c_str()); if (!std::isnan(psOptions->xres)) { CPLCreateXMLElementAndValue(psRoot, "ResX", CPLSPrintf("%.18g", psOptions->xres)); CPLCreateXMLElementAndValue(psRoot, "ResY", CPLSPrintf("%.18g", psOptions->yres)); } if (!std::isnan(psOptions->xmin)) { CPLCreateXMLElementAndValue(psRoot, "MinX", CPLSPrintf("%.18g", psOptions->xmin)); CPLCreateXMLElementAndValue(psRoot, "MinY", CPLSPrintf("%.18g", psOptions->ymin)); CPLCreateXMLElementAndValue(psRoot, "MaxX", CPLSPrintf("%.18g", psOptions->xmax)); CPLCreateXMLElementAndValue(psRoot, "MaxY", CPLSPrintf("%.18g", psOptions->ymax)); } int nBandCount = 0; if (!psOptions->osBandCount.empty()) { nBandCount = atoi(psOptions->osBandCount.c_str()); } else { if (!psOptions->osDataType.empty()) { nBandCount = std::max( nBandCount, CPLStringList(CSLTokenizeString2( psOptions->osDataType.c_str(), ", ", 0)) .size()); } if (!psOptions->osNodata.empty()) { nBandCount = std::max( nBandCount, CPLStringList(CSLTokenizeString2( psOptions->osNodata.c_str(), ", ", 0)) .size()); } if (!psOptions->osColorInterp.empty()) { nBandCount = std::max(nBandCount, CPLStringList( CSLTokenizeString2( psOptions->osColorInterp.c_str(), ", ", 0)) .size()); } } for (int i = 0; i < nBandCount; ++i) { auto psBand = CPLCreateXMLNode(psRoot, CXT_Element, "Band"); CPLAddXMLAttributeAndValue(psBand, "band", CPLSPrintf("%d", i + 1)); if (!psOptions->osDataType.empty()) { const CPLStringList aosTokens( CSLTokenizeString2(psOptions->osDataType.c_str(), ", ", 0)); if (aosTokens.size() == 1) CPLAddXMLAttributeAndValue(psBand, "dataType", aosTokens[0]); else if (i < aosTokens.size()) CPLAddXMLAttributeAndValue(psBand, "dataType", aosTokens[i]); } if (!psOptions->osNodata.empty()) { const CPLStringList aosTokens( CSLTokenizeString2(psOptions->osNodata.c_str(), ", ", 0)); if (aosTokens.size() == 1) CPLCreateXMLElementAndValue(psBand, "NoDataValue", aosTokens[0]); else if (i < aosTokens.size()) CPLCreateXMLElementAndValue(psBand, "NoDataValue", aosTokens[i]); } if (!psOptions->osColorInterp.empty()) { const CPLStringList aosTokens(CSLTokenizeString2( psOptions->osColorInterp.c_str(), ", ", 0)); if (aosTokens.size() == 1) CPLCreateXMLElementAndValue(psBand, "ColorInterp", aosTokens[0]); else if (i < aosTokens.size()) CPLCreateXMLElementAndValue(psBand, "ColorInterp", aosTokens[i]); } } if (psOptions->bMaskBand) { CPLCreateXMLElementAndValue(psRoot, "MaskBand", "true"); } int res = CPLSerializeXMLTreeToFile(psRoot, psOptions->osGTIFilename.c_str()); CPLDestroyXMLNode(psRoot); if (!res) return nullptr; } else { poLayer->SetMetadataItem("LOCATION_FIELD", psOptions->osLocationField.c_str()); if (!std::isnan(psOptions->xres)) { poLayer->SetMetadataItem("RESX", CPLSPrintf("%.18g", psOptions->xres)); poLayer->SetMetadataItem("RESY", CPLSPrintf("%.18g", psOptions->yres)); } if (!std::isnan(psOptions->xmin)) { poLayer->SetMetadataItem("MINX", CPLSPrintf("%.18g", psOptions->xmin)); poLayer->SetMetadataItem("MINY", CPLSPrintf("%.18g", psOptions->ymin)); poLayer->SetMetadataItem("MAXX", CPLSPrintf("%.18g", psOptions->xmax)); poLayer->SetMetadataItem("MAXY", CPLSPrintf("%.18g", psOptions->ymax)); } if (!psOptions->osBandCount.empty()) { poLayer->SetMetadataItem("BAND_COUNT", psOptions->osBandCount.c_str()); } if (!psOptions->osDataType.empty()) { poLayer->SetMetadataItem("DATA_TYPE", psOptions->osDataType.c_str()); } if (!psOptions->osNodata.empty()) { poLayer->SetMetadataItem("NODATA", psOptions->osNodata.c_str()); } if (!psOptions->osColorInterp.empty()) { poLayer->SetMetadataItem("COLOR_INTERPRETATION", psOptions->osColorInterp.c_str()); } if (psOptions->bMaskBand) { poLayer->SetMetadataItem("MASK_BAND", "YES"); } const CPLStringList aosMetadata(psOptions->aosMetadata); for (const auto &[pszKey, pszValue] : cpl::IterateNameValue(aosMetadata)) { poLayer->SetMetadataItem(pszKey, pszValue); } } const int ti_field = poLayerDefn->GetFieldIndex(psOptions->osLocationField.c_str()); if (ti_field < 0) { CPLError(CE_Failure, CPLE_AppDefined, "Unable to find field `%s' in file `%s'.", psOptions->osLocationField.c_str(), pszDest); return nullptr; } int i_SrcSRSName = -1; if (!psOptions->osSrcSRSFieldName.empty()) { i_SrcSRSName = poLayerDefn->GetFieldIndex(psOptions->osSrcSRSFieldName.c_str()); if (i_SrcSRSName < 0) { CPLError(CE_Failure, CPLE_AppDefined, "Unable to find field `%s' in file `%s'.", psOptions->osSrcSRSFieldName.c_str(), pszDest); return nullptr; } } // Load in memory existing file names in tile index. std::set oSetExistingFiles; OGRSpatialReference oAlreadyExistingSRS; if (bExistingLayer) { for (auto &&poFeature : poLayer) { if (poFeature->IsFieldSetAndNotNull(ti_field)) { if (oSetExistingFiles.empty()) { auto poSrcDS = std::unique_ptr(GDALDataset::Open( poFeature->GetFieldAsString(ti_field), GDAL_OF_RASTER, nullptr, nullptr, nullptr)); if (poSrcDS) { auto poSrcSRS = poSrcDS->GetSpatialRef(); if (poSrcSRS) oAlreadyExistingSRS = *poSrcSRS; } } oSetExistingFiles.insert(poFeature->GetFieldAsString(ti_field)); } } } std::string osCurrentPath; if (psOptions->bWriteAbsolutePath) { char *pszCurrentPath = CPLGetCurrentDir(); if (pszCurrentPath) { osCurrentPath = pszCurrentPath; } else { CPLError(CE_Warning, CPLE_AppDefined, "This system does not support the CPLGetCurrentDir call. " "The option -bWriteAbsolutePath will have no effect."); } CPLFree(pszCurrentPath); } const bool bIsGTIContext = !std::isnan(psOptions->xres) || !std::isnan(psOptions->xmin) || !psOptions->osBandCount.empty() || !psOptions->osNodata.empty() || !psOptions->osColorInterp.empty() || !psOptions->osDataType.empty() || psOptions->bMaskBand || !psOptions->aosMetadata.empty() || !psOptions->osGTIFilename.empty(); /* -------------------------------------------------------------------- */ /* loop over GDAL files, processing. */ /* -------------------------------------------------------------------- */ int iCur = 0; int nTotal = nSrcCount + 1; while (true) { const std::string osSrcFilename = oGDALTileIndexTileIterator.next(); if (osSrcFilename.empty()) break; if (bSkipFirstTile) { bSkipFirstTile = false; continue; } std::string osFileNameToWrite; VSIStatBuf sStatBuf; // Make sure it is a file before building absolute path name. if (!osCurrentPath.empty() && CPLIsFilenameRelative(osSrcFilename.c_str()) && VSIStat(osSrcFilename.c_str(), &sStatBuf) == 0) { osFileNameToWrite = CPLProjectRelativeFilenameSafe( osCurrentPath.c_str(), osSrcFilename.c_str()); } else { osFileNameToWrite = osSrcFilename.c_str(); } // Checks that file is not already in tileindex. if (oSetExistingFiles.find(osFileNameToWrite) != oSetExistingFiles.end()) { CPLError(CE_Warning, CPLE_AppDefined, "File %s is already in tileindex. Skipping it.", osFileNameToWrite.c_str()); continue; } std::unique_ptr poSrcDS; { std::unique_ptr poFailureIntoWarning; if (!bFailOnErrors) poFailureIntoWarning = std::make_unique(); CPL_IGNORE_RET_VAL(poFailureIntoWarning); poSrcDS.reset(GDALDataset::Open( osSrcFilename.c_str(), GDAL_OF_RASTER | GDAL_OF_VERBOSE_ERROR, nullptr, nullptr, nullptr)); if (poSrcDS == nullptr) { CPLError(bFailOnErrors ? CE_Failure : CE_Warning, CPLE_AppDefined, "Unable to open %s%s.", osSrcFilename.c_str(), bFailOnErrors ? "" : ", skipping"); if (bFailOnErrors) return nullptr; continue; } } GDALGeoTransform gt; if (poSrcDS->GetGeoTransform(gt) != CE_None) { CPLError(bFailOnErrors ? CE_Failure : CE_Warning, CPLE_AppDefined, "It appears no georeferencing is available for\n" "`%s'%s.", osSrcFilename.c_str(), bFailOnErrors ? "" : ", skipping"); if (bFailOnErrors) return nullptr; continue; } auto poSrcSRS = poSrcDS->GetSpatialRef(); // If not set target srs, test that the current file uses same // projection as others. if (oTargetSRS.IsEmpty()) { if (!oAlreadyExistingSRS.IsEmpty()) { if (poSrcSRS == nullptr || !poSrcSRS->IsSame(&oAlreadyExistingSRS)) { CPLError( CE_Warning, CPLE_AppDefined, "%s is not using the same projection system " "as other files in the tileindex.\n" "This may cause problems when using it in MapServer " "for example.\n" "Use -t_srs option to set target projection system. %s", osSrcFilename.c_str(), psOptions->bSkipDifferentProjection ? "Skipping this file." : ""); if (psOptions->bSkipDifferentProjection) { continue; } } } else { if (poSrcSRS) oAlreadyExistingSRS = *poSrcSRS; } } const int nXSize = poSrcDS->GetRasterXSize(); const int nYSize = poSrcDS->GetRasterYSize(); if (nXSize == 0 || nYSize == 0) { CPLError(bFailOnErrors ? CE_Failure : CE_Warning, CPLE_AppDefined, "%s has 0 width or height%s", osSrcFilename.c_str(), bFailOnErrors ? "" : ", skipping"); if (bFailOnErrors) return nullptr; continue; } double adfX[5] = {0.0, 0.0, 0.0, 0.0, 0.0}; double adfY[5] = {0.0, 0.0, 0.0, 0.0, 0.0}; adfX[0] = gt[0] + 0 * gt[1] + 0 * gt[2]; adfY[0] = gt[3] + 0 * gt[4] + 0 * gt[5]; adfX[1] = gt[0] + nXSize * gt[1] + 0 * gt[2]; adfY[1] = gt[3] + nXSize * gt[4] + 0 * gt[5]; adfX[2] = gt[0] + nXSize * gt[1] + nYSize * gt[2]; adfY[2] = gt[3] + nXSize * gt[4] + nYSize * gt[5]; adfX[3] = gt[0] + 0 * gt[1] + nYSize * gt[2]; adfY[3] = gt[3] + 0 * gt[4] + nYSize * gt[5]; adfX[4] = gt[0] + 0 * gt[1] + 0 * gt[2]; adfY[4] = gt[3] + 0 * gt[4] + 0 * gt[5]; // If set target srs, compute the reprojected extent. // Use GDALWarp() with VRT output so the stored bounding box matches // what the GTI reader computes in its Pass 1 (GetSourceDesc). // A plain 4-corner transform or standalone GDALSuggestedWarpOutput2() // can underestimate the extent because GDALWarp() may adjust the // resolution after calling GDALSuggestedWarpOutput2(), producing a // different pixel count and thus a different extent. if (!oTargetSRS.IsEmpty() && poSrcSRS) { if (!poSrcSRS->IsSame(&oTargetSRS)) { char *pszDstWKT = nullptr; oTargetSRS.exportToWkt(&pszDstWKT); CPLStringList aosWarpArgs; aosWarpArgs.AddString("-of"); aosWarpArgs.AddString("VRT"); aosWarpArgs.AddString("-t_srs"); aosWarpArgs.AddString(pszDstWKT); CPLFree(pszDstWKT); GDALWarpAppOptions *psWarpOptions = GDALWarpAppOptionsNew(aosWarpArgs.List(), nullptr); GDALDatasetH hSrcDS = GDALDataset::ToHandle(poSrcDS.get()); GDALDatasetH ahSrcDS[] = {hSrcDS}; int bUsageError = false; auto poWarpDS = std::unique_ptr( GDALDataset::FromHandle(GDALWarp( "", nullptr, 1, ahSrcDS, psWarpOptions, &bUsageError))); GDALWarpAppOptionsFree(psWarpOptions); if (!poWarpDS) { CPLError(bFailOnErrors ? CE_Failure : CE_Warning, CPLE_AppDefined, "unable to compute reprojected extent from " "source SRS `%s' to target SRS `%s' for " "file `%s'%s", poSrcDS->GetProjectionRef(), psOptions->osTargetSRS.c_str(), osFileNameToWrite.c_str(), bFailOnErrors ? "" : ", skipping"); if (bFailOnErrors) return nullptr; continue; } GDALGeoTransform warpGT; poWarpDS->GetGeoTransform(warpGT); const int nDstPixels = poWarpDS->GetRasterXSize(); const int nDstLines = poWarpDS->GetRasterYSize(); const double dfDstMinX = warpGT.xorig; const double dfDstMaxY = warpGT.yorig; const double dfDstMaxX = warpGT.xorig + nDstPixels * warpGT.xscale; const double dfDstMinY = warpGT.yorig + nDstLines * warpGT.yscale; adfX[0] = dfDstMinX; adfY[0] = dfDstMaxY; adfX[1] = dfDstMaxX; adfY[1] = dfDstMaxY; adfX[2] = dfDstMaxX; adfY[2] = dfDstMinY; adfX[3] = dfDstMinX; adfY[3] = dfDstMinY; adfX[4] = dfDstMinX; adfY[4] = dfDstMaxY; } } else if (bIsGTIContext && !oAlreadyExistingSRS.IsEmpty() && (poSrcSRS == nullptr || !poSrcSRS->IsSame(&oAlreadyExistingSRS))) { CPLError( CE_Failure, CPLE_AppDefined, "%s is not using the same projection system " "as other files in the tileindex. This is not compatible of " "GTI use. Use -t_srs option to reproject tile extents " "to a common SRS.", osSrcFilename.c_str()); return nullptr; } const double dfMinX = std::min(std::min(adfX[0], adfX[1]), std::min(adfX[2], adfX[3])); const double dfMinY = std::min(std::min(adfY[0], adfY[1]), std::min(adfY[2], adfY[3])); const double dfMaxX = std::max(std::max(adfX[0], adfX[1]), std::max(adfX[2], adfX[3])); const double dfMaxY = std::max(std::max(adfY[0], adfY[1]), std::max(adfY[2], adfY[3])); const double dfRes = sqrt((dfMaxX - dfMinX) * (dfMaxY - dfMinY) / nXSize / nYSize); if (!std::isnan(psOptions->dfMinPixelSize) && dfRes < psOptions->dfMinPixelSize) { CPLError(CE_Warning, CPLE_AppDefined, "%s has %f as pixel size (< %f). Skipping", osSrcFilename.c_str(), dfRes, psOptions->dfMinPixelSize); continue; } if (!std::isnan(psOptions->dfMaxPixelSize) && dfRes > psOptions->dfMaxPixelSize) { CPLError(CE_Warning, CPLE_AppDefined, "%s has %f as pixel size (> %f). Skipping", osSrcFilename.c_str(), dfRes, psOptions->dfMaxPixelSize); continue; } auto poFeature = std::make_unique(poLayerDefn); poFeature->SetField(ti_field, osFileNameToWrite.c_str()); if (i_SrcSRSName >= 0 && poSrcSRS) { const char *pszAuthorityCode = poSrcSRS->GetAuthorityCode(nullptr); const char *pszAuthorityName = poSrcSRS->GetAuthorityName(nullptr); if (psOptions->eSrcSRSFormat == FORMAT_AUTO) { if (pszAuthorityName != nullptr && pszAuthorityCode != nullptr) { poFeature->SetField(i_SrcSRSName, CPLSPrintf("%s:%s", pszAuthorityName, pszAuthorityCode)); } else if (nMaxFieldSize == 0 || strlen(poSrcDS->GetProjectionRef()) <= static_cast(nMaxFieldSize)) { poFeature->SetField(i_SrcSRSName, poSrcDS->GetProjectionRef()); } else { char *pszProj4 = nullptr; if (poSrcSRS->exportToProj4(&pszProj4) == OGRERR_NONE) { poFeature->SetField(i_SrcSRSName, pszProj4); } else { poFeature->SetField(i_SrcSRSName, poSrcDS->GetProjectionRef()); } CPLFree(pszProj4); } } else if (psOptions->eSrcSRSFormat == FORMAT_WKT) { if (nMaxFieldSize == 0 || strlen(poSrcDS->GetProjectionRef()) <= static_cast(nMaxFieldSize)) { poFeature->SetField(i_SrcSRSName, poSrcDS->GetProjectionRef()); } else { CPLError(CE_Warning, CPLE_AppDefined, "Cannot write WKT for file %s as it is too long!", osFileNameToWrite.c_str()); } } else if (psOptions->eSrcSRSFormat == FORMAT_PROJ) { char *pszProj4 = nullptr; if (poSrcSRS->exportToProj4(&pszProj4) == OGRERR_NONE) { poFeature->SetField(i_SrcSRSName, pszProj4); } CPLFree(pszProj4); } else if (psOptions->eSrcSRSFormat == FORMAT_EPSG) { if (pszAuthorityName != nullptr && pszAuthorityCode != nullptr) poFeature->SetField(i_SrcSRSName, CPLSPrintf("%s:%s", pszAuthorityName, pszAuthorityCode)); } } for (const auto &oFetchMD : psOptions->aoFetchMD) { if (EQUAL(oFetchMD.osRasterItemName.c_str(), "{PIXEL_SIZE}")) { poFeature->SetField(oFetchMD.osFieldName.c_str(), dfRes); continue; } const char *pszMD = poSrcDS->GetMetadataItem(oFetchMD.osRasterItemName.c_str()); if (pszMD) { if (EQUAL(oFetchMD.osRasterItemName.c_str(), "TIFFTAG_DATETIME")) { int nYear, nMonth, nDay, nHour, nMin, nSec; if (sscanf(pszMD, "%04d:%02d:%02d %02d:%02d:%02d", &nYear, &nMonth, &nDay, &nHour, &nMin, &nSec) == 6) { poFeature->SetField( oFetchMD.osFieldName.c_str(), CPLSPrintf("%04d/%02d/%02d %02d:%02d:%02d", nYear, nMonth, nDay, nHour, nMin, nSec)); continue; } } poFeature->SetField(oFetchMD.osFieldName.c_str(), pszMD); } } auto poPoly = std::make_unique(); auto poRing = std::make_unique(); for (int k = 0; k < 5; k++) poRing->addPoint(adfX[k], adfY[k]); poPoly->addRing(std::move(poRing)); poFeature->SetGeometryDirectly(poPoly.release()); if (poLayer->CreateFeature(poFeature.get()) != OGRERR_NONE) { CPLError(CE_Failure, CPLE_AppDefined, "Failed to create feature in tile index."); return nullptr; } ++iCur; if (psOptions->pfnProgress && !psOptions->pfnProgress(static_cast(iCur) / nTotal, "", psOptions->pProgressData)) { return nullptr; } if (iCur >= nSrcCount) ++nTotal; } if (psOptions->pfnProgress) psOptions->pfnProgress(1.0, "", psOptions->pProgressData); if (poTileIndexDSUnique) return GDALDataset::ToHandle(poTileIndexDSUnique.release()); else return GDALDataset::ToHandle(poTileIndexDS); } /************************************************************************/ /* SanitizeSRS */ /************************************************************************/ static char *SanitizeSRS(const char *pszUserInput) { OGRSpatialReferenceH hSRS; char *pszResult = nullptr; CPLErrorReset(); hSRS = OSRNewSpatialReference(nullptr); if (OSRSetFromUserInput(hSRS, pszUserInput) == OGRERR_NONE) OSRExportToWkt(hSRS, &pszResult); else { CPLError(CE_Failure, CPLE_AppDefined, "Translating SRS failed:\n%s", pszUserInput); } OSRDestroySpatialReference(hSRS); return pszResult; } /************************************************************************/ /* GDALTileIndexOptionsNew() */ /************************************************************************/ /** * Allocates a GDALTileIndexOptions struct. * * @param papszArgv NULL terminated list of options (potentially including * filename and open options too), or NULL. The accepted options are the ones of * the gdaltindex utility. * @param psOptionsForBinary (output) may be NULL (and should generally be * NULL), otherwise (gdaltindex_bin.cpp use case) must be allocated with * GDALTileIndexOptionsForBinaryNew() prior to this function. Will be filled * with potentially present filename, open options,... * @return pointer to the allocated GDALTileIndexOptions struct. Must be freed * with GDALTileIndexOptionsFree(). * * @since GDAL 3.9 */ GDALTileIndexOptions * GDALTileIndexOptionsNew(char **papszArgv, GDALTileIndexOptionsForBinary *psOptionsForBinary) { auto psOptions = std::make_unique(); /* -------------------------------------------------------------------- */ /* Parse arguments. */ /* -------------------------------------------------------------------- */ CPLStringList aosArgv; if (papszArgv) { const int nArgc = CSLCount(papszArgv); for (int i = 0; i < nArgc; i++) { aosArgv.AddString(papszArgv[i]); } } try { auto argParser = GDALTileIndexAppOptionsGetParser(psOptions.get(), psOptionsForBinary); argParser->parse_args_without_binary_name(aosArgv.List()); // Check all no store_into args if (auto oTr = argParser->present>("-tr")) { psOptions->xres = (*oTr)[0]; psOptions->yres = (*oTr)[1]; } if (auto oTargetExtent = argParser->present>("-te")) { psOptions->xmin = (*oTargetExtent)[0]; psOptions->ymin = (*oTargetExtent)[1]; psOptions->xmax = (*oTargetExtent)[2]; psOptions->ymax = (*oTargetExtent)[3]; } if (auto fetchMd = argParser->present>("-fetch_md")) { CPLAssert(fetchMd->size() % 3 == 0); // Loop for (size_t i = 0; i < fetchMd->size(); i += 3) { OGRFieldType type; const auto &typeName{fetchMd->at(i + 2)}; if (typeName == "String") { type = OFTString; } else if (typeName == "Integer") { type = OFTInteger; } else if (typeName == "Integer64") { type = OFTInteger64; } else if (typeName == "Real") { type = OFTReal; } else if (typeName == "Date") { type = OFTDate; } else if (typeName == "DateTime") { type = OFTDateTime; } else { CPLError(CE_Failure, CPLE_AppDefined, "-fetch_md requires a valid type name as third " "argument: %s was given.", fetchMd->at(i).c_str()); return nullptr; } const GDALTileIndexRasterMetadata oMD{type, fetchMd->at(i + 1), fetchMd->at(i)}; psOptions->aoFetchMD.push_back(std::move(oMD)); } } // Check -t_srs if (!psOptions->osTargetSRS.empty()) { auto sanitized{SanitizeSRS(psOptions->osTargetSRS.c_str())}; if (sanitized) { psOptions->osTargetSRS = sanitized; CPLFree(sanitized); } else { // Error was already reported by SanitizeSRS, just return nullptr psOptions->osTargetSRS.clear(); return nullptr; } } } catch (const std::exception &error) { CPLError(CE_Failure, CPLE_AppDefined, "%s", error.what()); return nullptr; } return psOptions.release(); } /************************************************************************/ /* GDALTileIndexOptionsFree() */ /************************************************************************/ /** * Frees the GDALTileIndexOptions struct. * * @param psOptions the options struct for GDALTileIndex(). * * @since GDAL 3.9 */ void GDALTileIndexOptionsFree(GDALTileIndexOptions *psOptions) { delete psOptions; } /************************************************************************/ /* GDALTileIndexOptionsSetProgress() */ /************************************************************************/ /** * Set a progress function. * * @param psOptions the options struct for GDALTileIndex(). * @param pfnProgress the progress callback. * @param pProgressData the user data for the progress callback. * * @since GDAL 3.11 */ void GDALTileIndexOptionsSetProgress(GDALTileIndexOptions *psOptions, GDALProgressFunc pfnProgress, void *pProgressData) { psOptions->pfnProgress = pfnProgress; psOptions->pProgressData = pProgressData; } #undef CHECK_HAS_ENOUGH_ADDITIONAL_ARGS