From 0df47a04519b61f510fd2f9d89ffc5e5dfd894b9 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sun, 3 Nov 2024 12:49:04 +0100 Subject: [PATCH] Add GDALRasterComputeMinMaxLocation / GDALRasterBand::ComputeMinMaxLocation, and map it to SWIG and add a gdal_minmax_location.py script: ``` $ python ../swig/python/gdal-utils/osgeo_utils/samples/gdal_minmax_location.py byte.tif Minimum=74.0 at (col,line)=(9,17), (X,Y)_georef=(441290.0,3750270.0), (long,lat)_WGS84=(-117.6358076,33.8929309) Maximum=255.0 at (col,line)=(2,18), (X,Y)_georef=(440870.0,3750210.0), (long,lat)_WGS84=(-117.6403456,33.8923662) ``` --- autotest/cpp/test_gdal.cpp | 76 ++++++ autotest/gcore/basic_test.py | 19 ++ doc/source/api/python_samples.rst | 1 + gcore/gdal.h | 4 + gcore/gdal_priv.h | 3 + gcore/gdalrasterband.cpp | 244 ++++++++++++++++++ swig/include/Band.i | 17 ++ swig/include/python/gdal_python.i | 41 +++ .../samples/gdal_minmax_location.py | 128 +++++++++ 9 files changed, 533 insertions(+) create mode 100644 swig/python/gdal-utils/osgeo_utils/samples/gdal_minmax_location.py diff --git a/autotest/cpp/test_gdal.cpp b/autotest/cpp/test_gdal.cpp index 45ebfb4ab870..813d983b66d5 100644 --- a/autotest/cpp/test_gdal.cpp +++ b/autotest/cpp/test_gdal.cpp @@ -4777,4 +4777,80 @@ TEST_F(test_gdal, ReadRaster) } } +// Test GDALComputeRasterMinMaxLocation +TEST_F(test_gdal, GDALComputeRasterMinMaxLocation) +{ + GDALDatasetH hDS = GDALOpen(GCORE_DATA_DIR "byte.tif", GA_ReadOnly); + ASSERT_NE(hDS, nullptr); + GDALRasterBandH hBand = GDALGetRasterBand(hDS, 1); + double dfMin = 0; + double dfMax = 0; + int nMinX = -1; + int nMinY = -1; + int nMaxX = -1; + int nMaxY = -1; + EXPECT_EQ(GDALComputeRasterMinMaxLocation(hBand, &dfMin, &dfMax, &nMinX, + &nMinY, &nMaxX, &nMaxY), + CE_None); + EXPECT_EQ(dfMin, 74.0); + EXPECT_EQ(dfMax, 255.0); + EXPECT_EQ(nMinX, 9); + EXPECT_EQ(nMinY, 17); + EXPECT_EQ(nMaxX, 2); + EXPECT_EQ(nMaxY, 18); + GByte val = 0; + EXPECT_EQ(GDALRasterIO(hBand, GF_Read, nMinX, nMinY, 1, 1, &val, 1, 1, + GDT_Byte, 0, 0), + CE_None); + EXPECT_EQ(val, 74); + EXPECT_EQ(GDALRasterIO(hBand, GF_Read, nMaxX, nMaxY, 1, 1, &val, 1, 1, + GDT_Byte, 0, 0), + CE_None); + EXPECT_EQ(val, 255); + GDALClose(hDS); +} + +// Test GDALComputeRasterMinMaxLocation +TEST_F(test_gdal, GDALComputeRasterMinMaxLocation_with_mask) +{ + GDALDatasetUniquePtr poDS(GDALDriver::FromHandle(GDALGetDriverByName("MEM")) + ->Create("", 2, 2, 1, GDT_Byte, nullptr)); + std::array buffer = { + 2, 10, ////////////////////////////////////////////////////////// + 4, 20, ////////////////////////////////////////////////////////// + }; + GDALRasterIOExtraArg sExtraArg; + INIT_RASTERIO_EXTRA_ARG(sExtraArg); + EXPECT_EQ(poDS->GetRasterBand(1)->RasterIO( + GF_Write, 0, 0, 2, 2, buffer.data(), 2, 2, GDT_Byte, + sizeof(uint8_t), 2 * sizeof(uint8_t), &sExtraArg), + CE_None); + + poDS->GetRasterBand(1)->CreateMaskBand(0); + std::array buffer_mask = { + 0, 255, ////////////////////////////////////////////////////////// + 255, 0, ////////////////////////////////////////////////////////// + }; + EXPECT_EQ(poDS->GetRasterBand(1)->GetMaskBand()->RasterIO( + GF_Write, 0, 0, 2, 2, buffer_mask.data(), 2, 2, GDT_Byte, + sizeof(uint8_t), 2 * sizeof(uint8_t), &sExtraArg), + CE_None); + + double dfMin = 0; + double dfMax = 0; + int nMinX = -1; + int nMinY = -1; + int nMaxX = -1; + int nMaxY = -1; + EXPECT_EQ(poDS->GetRasterBand(1)->ComputeRasterMinMaxLocation( + &dfMin, &dfMax, &nMinX, &nMinY, &nMaxX, &nMaxY), + CE_None); + EXPECT_EQ(dfMin, 4); + EXPECT_EQ(dfMax, 10); + EXPECT_EQ(nMinX, 0); + EXPECT_EQ(nMinY, 1); + EXPECT_EQ(nMaxX, 1); + EXPECT_EQ(nMaxY, 0); +} + } // namespace diff --git a/autotest/gcore/basic_test.py b/autotest/gcore/basic_test.py index 3b399d820925..33f9e58a36ab 100755 --- a/autotest/gcore/basic_test.py +++ b/autotest/gcore/basic_test.py @@ -987,3 +987,22 @@ def test_colorinterp(): assert name not in d d[name] = c assert gdal.GetColorInterpretationByName(name) == c + + +def test_ComputeMinMaxLocation(): + + ds = gdal.Open("data/byte.tif") + ret = ds.GetRasterBand(1).ComputeMinMaxLocation() + assert ( + ret.min == 74 + and ret.max == 255 + and ret.minX == 9 + and ret.minY == 17 + and ret.maxX == 2 + and ret.maxY == 18 + ) + + ds = gdal.GetDriverByName("MEM").Create("", 1, 1, 1, gdal.GDT_Float64) + ds.GetRasterBand(1).Fill(float("nan")) + ret = ds.GetRasterBand(1).ComputeMinMaxLocation() + assert ret is None diff --git a/doc/source/api/python_samples.rst b/doc/source/api/python_samples.rst index a2831f63ffaa..294a2bce48c4 100644 --- a/doc/source/api/python_samples.rst +++ b/doc/source/api/python_samples.rst @@ -47,6 +47,7 @@ Python Raster Sample scripts - hsv_merge: Merge greyscale image into RGB image as intensity in HSV space. - gdal_ls: Display the list of files in a virtual directory, like /vsicurl or /vsizip - gdal_cp: Copy a virtual file + - gdal_minmax_location: returns the location where min/max values of a raster are hit. Python Vector Sample scripts ------------------------------ diff --git a/gcore/gdal.h b/gcore/gdal.h index f7a411c393cb..8638f8501b27 100644 --- a/gcore/gdal.h +++ b/gcore/gdal.h @@ -1670,6 +1670,10 @@ CPLErr CPL_DLL CPL_STDCALL GDALSetRasterScale(GDALRasterBandH hBand, CPLErr CPL_DLL CPL_STDCALL GDALComputeRasterMinMax(GDALRasterBandH hBand, int bApproxOK, double adfMinMax[2]); +CPLErr CPL_DLL GDALComputeRasterMinMaxLocation(GDALRasterBandH hBand, + double *pdfMin, double *pdfMax, + int *pnMinX, int *pnMinY, + int *pnMaxX, int *pnMaxY); CPLErr CPL_DLL CPL_STDCALL GDALFlushRasterCache(GDALRasterBandH hBand); CPLErr CPL_DLL CPL_STDCALL GDALDropRasterCache(GDALRasterBandH hBand); CPLErr CPL_DLL CPL_STDCALL GDALGetRasterHistogram( diff --git a/gcore/gdal_priv.h b/gcore/gdal_priv.h index 20dcc3734cf3..c615a3a90319 100644 --- a/gcore/gdal_priv.h +++ b/gcore/gdal_priv.h @@ -1810,6 +1810,9 @@ class CPL_DLL GDALRasterBand : public GDALMajorObject virtual CPLErr SetStatistics(double dfMin, double dfMax, double dfMean, double dfStdDev); virtual CPLErr ComputeRasterMinMax(int bApproxOK, double *adfMinMax); + virtual CPLErr ComputeRasterMinMaxLocation(double *pdfMin, double *pdfMax, + int *pnMinX, int *pnMinY, + int *pnMaxX, int *pnMaxY); // Only defined when Doxygen enabled #ifdef DOXYGEN_SKIP diff --git a/gcore/gdalrasterband.cpp b/gcore/gdalrasterband.cpp index bf2da7a2893f..36c1165cf3ca 100644 --- a/gcore/gdalrasterband.cpp +++ b/gcore/gdalrasterband.cpp @@ -38,6 +38,7 @@ #include "gdal_rat.h" #include "gdal_priv_templates.hpp" #include "gdal_interpolateatpoint.h" +#include "gdal_minmax_element.hpp" /************************************************************************/ /* GDALRasterBand() */ @@ -7415,6 +7416,249 @@ CPLErr CPL_STDCALL GDALComputeRasterMinMax(GDALRasterBandH hBand, int bApproxOK, return poBand->ComputeRasterMinMax(bApproxOK, adfMinMax); } +/************************************************************************/ +/* ComputeRasterMinMaxLocation() */ +/************************************************************************/ + +/** + * \brief Compute the min/max values for a band, and their location. + * + * If they are hit in several locations, it is not specified which one will be + * returned. + * + * @param[out] pdfMin Pointer to the minimum value. + * @param[out] pdfMax Pointer to the maximum value. + * @param[out] pnMinX Pointer to the column where the minimum value is hit. + * @param[out] pnMinY Pointer to the line where the minimum value is hit. + * @param[out] pnMaxX Pointer to the column where the maximum value is hit. + * @param[out] pnMaxY Pointer to the line where the maximum value is hit. + * + * @return CE_None in case of success, CE_Warning if there are no valid values, + * CE_Failure in case of error. + * + * @since GDAL 3.11 + */ + +CPLErr GDALRasterBand::ComputeRasterMinMaxLocation(double *pdfMin, + double *pdfMax, int *pnMinX, + int *pnMinY, int *pnMaxX, + int *pnMaxY) +{ + int nMinX = -1; + int nMinY = -1; + int nMaxX = -1; + int nMaxY = -1; + double dfMin = std::numeric_limits::infinity(); + double dfMax = -std::numeric_limits::infinity(); + if (pdfMin) + *pdfMin = dfMin; + if (pdfMax) + *pdfMax = dfMax; + if (pnMinX) + *pnMinX = nMinX; + if (pnMinY) + *pnMinY = nMinY; + if (pnMaxX) + *pnMaxX = nMaxX; + if (pnMaxY) + *pnMaxY = nMaxY; + + if (GDALDataTypeIsComplex(eDataType)) + { + CPLError(CE_Failure, CPLE_NotSupported, + "Complex data type not supported"); + return CE_Failure; + } + + int bGotNoDataValue = FALSE; + const double dfNoDataValue = GetNoDataValue(&bGotNoDataValue); + bGotNoDataValue = bGotNoDataValue && !std::isnan(dfNoDataValue); + bool bGotFloatNoDataValue = false; + float fNoDataValue = 0.0f; + ComputeFloatNoDataValue(eDataType, dfNoDataValue, bGotNoDataValue, + fNoDataValue, bGotFloatNoDataValue); + + GDALRasterBand *poMaskBand = nullptr; + if (!bGotNoDataValue) + { + const int l_nMaskFlags = GetMaskFlags(); + if (l_nMaskFlags != GMF_ALL_VALID && l_nMaskFlags != GMF_NODATA && + GetColorInterpretation() != GCI_AlphaBand) + { + poMaskBand = GetMaskBand(); + } + } + + bool bSignedByte = false; + if (eDataType == GDT_Byte) + { + EnablePixelTypeSignedByteWarning(false); + const char *pszPixelType = + GetMetadataItem("PIXELTYPE", "IMAGE_STRUCTURE"); + EnablePixelTypeSignedByteWarning(true); + bSignedByte = + pszPixelType != nullptr && EQUAL(pszPixelType, "SIGNEDBYTE"); + } + + GByte *pabyMaskData = nullptr; + if (poMaskBand) + { + pabyMaskData = + static_cast(VSI_MALLOC2_VERBOSE(nBlockXSize, nBlockYSize)); + if (!pabyMaskData) + { + return CE_Failure; + } + } + + if (!InitBlockInfo()) + return CE_Failure; + + const GIntBig nTotalBlocks = + static_cast(nBlocksPerRow) * nBlocksPerColumn; + for (GIntBig iBlock = 0; iBlock < nTotalBlocks; ++iBlock) + { + const int iYBlock = static_cast(iBlock / nBlocksPerRow); + const int iXBlock = static_cast(iBlock % nBlocksPerRow); + + GDALRasterBlock *poBlock = GetLockedBlockRef(iXBlock, iYBlock); + if (poBlock == nullptr) + { + CPLFree(pabyMaskData); + return CE_Failure; + } + + void *const pData = poBlock->GetDataRef(); + + int nXCheck = 0, nYCheck = 0; + GetActualBlockSize(iXBlock, iYBlock, &nXCheck, &nYCheck); + + if (poMaskBand && + poMaskBand->RasterIO(GF_Read, iXBlock * nBlockXSize, + iYBlock * nBlockYSize, nXCheck, nYCheck, + pabyMaskData, nXCheck, nYCheck, GDT_Byte, 0, + nBlockXSize, nullptr) != CE_None) + { + poBlock->DropLock(); + CPLFree(pabyMaskData); + return CE_Failure; + } + + if (poMaskBand || nYCheck < nBlockYSize || nXCheck < nBlockXSize) + { + for (int iY = 0; iY < nYCheck; ++iY) + { + for (int iX = 0; iX < nXCheck; ++iX) + { + const GPtrDiff_t iOffset = + iX + static_cast(iY) * nBlockXSize; + if (pabyMaskData && pabyMaskData[iOffset] == 0) + continue; + bool bValid = true; + double dfValue = GetPixelValue( + eDataType, bSignedByte, pData, iOffset, bGotNoDataValue, + dfNoDataValue, bGotFloatNoDataValue, fNoDataValue, + bValid); + if (!bValid) + continue; + if (dfValue < dfMin) + { + dfMin = dfValue; + nMinX = iXBlock * nBlockXSize + iX; + nMinY = iYBlock * nBlockYSize + iY; + } + if (dfValue > dfMax) + { + dfMax = dfValue; + nMaxX = iXBlock * nBlockXSize + iX; + nMaxY = iYBlock * nBlockYSize + iY; + } + } + } + } + else + { + const auto [pos_min, pos_max] = gdal::minmax_element( + pData, static_cast(nBlockXSize) * nBlockYSize, + bSignedByte ? GDT_Int8 : eDataType, bGotNoDataValue, + dfNoDataValue); + const int nMinXBlock = static_cast(pos_min % nBlockXSize); + const int nMinYBlock = static_cast(pos_min / nBlockXSize); + const int nMaxXBlock = static_cast(pos_max % nBlockXSize); + const int nMaxYBlock = static_cast(pos_max / nBlockXSize); + bool bValid = true; + const double dfMinValueBlock = GetPixelValue( + eDataType, bSignedByte, pData, pos_min, bGotNoDataValue, + dfNoDataValue, bGotFloatNoDataValue, fNoDataValue, bValid); + if (!bValid) + continue; + const double dfMaxValueBlock = GetPixelValue( + eDataType, bSignedByte, pData, pos_max, bGotNoDataValue, + dfNoDataValue, bGotFloatNoDataValue, fNoDataValue, bValid); + // No need to check for bValid here since if we got a valid min, + // we'll also get a valid max + if (dfMinValueBlock < dfMin) + { + dfMin = dfMinValueBlock; + nMinX = nMinXBlock; + nMinY = nMinYBlock; + } + if (dfMaxValueBlock > dfMax) + { + dfMax = dfMaxValueBlock; + nMaxX = nMaxXBlock; + nMaxY = nMaxYBlock; + } + } + + poBlock->DropLock(); + + if (eDataType == GDT_Byte && dfMin == 0 && dfMax == 255) + { + break; + } + } + + CPLFree(pabyMaskData); + + if (pdfMin) + *pdfMin = dfMin; + if (pdfMax) + *pdfMax = dfMax; + if (pnMinX) + *pnMinX = nMinX; + if (pnMinY) + *pnMinY = nMinY; + if (pnMaxX) + *pnMaxX = nMaxX; + if (pnMaxY) + *pnMaxY = nMaxY; + return nMinX < 0 ? CE_Warning : CE_None; +} + +/************************************************************************/ +/* GDALComputeRasterMinMaxLocation() */ +/************************************************************************/ + +/** + * \brief Compute the min/max values for a band, and their location. + * + * @see GDALRasterBand::ComputeRasterMinMax() + * @since GDAL 3.11 + */ + +CPLErr GDALComputeRasterMinMaxLocation(GDALRasterBandH hBand, double *pdfMin, + double *pdfMax, int *pnMinX, int *pnMinY, + int *pnMaxX, int *pnMaxY) + +{ + VALIDATE_POINTER1(hBand, "GDALComputeRasterMinMaxLocation", CE_Failure); + + GDALRasterBand *poBand = GDALRasterBand::FromHandle(hBand); + return poBand->ComputeRasterMinMaxLocation(pdfMin, pdfMax, pnMinX, pnMinY, + pnMaxX, pnMaxY); +} + /************************************************************************/ /* SetDefaultHistogram() */ /************************************************************************/ diff --git a/swig/include/Band.i b/swig/include/Band.i index ea4fd47f6e22..6c45f714e078 100644 --- a/swig/include/Band.i +++ b/swig/include/Band.i @@ -669,6 +669,23 @@ CPLErr AdviseRead( int xoff, int yoff, int xsize, int ysize, %clear (CPLErr); #endif +%apply (double *OUTPUT){double *pdfMin, double *pdfMax}; +%apply (int *OUTPUT){int *pnMinX, int *pnMinY}; +%apply (int *OUTPUT){int *pnMaxX, int *pnMaxY}; +#if !defined(SWIGPYTHON) +%apply (IF_ERROR_RETURN_NONE) { (CPLErr) }; +#endif + CPLErr ComputeMinMaxLocation( double *pdfMin, double *pdfMax, + int *pnMinX, int *pnMinY, + int *pnMaxX, int *pnMaxY ) { + return GDALComputeRasterMinMaxLocation( self, pdfMin, pdfMax, + pnMinX, pnMinY, + pnMaxX, pnMaxY ); + } +#if !defined(SWIGPYTHON) +%clear (CPLErr); +#endif + %newobject AsMDArray; GDALMDArrayHS *AsMDArray() { diff --git a/swig/include/python/gdal_python.i b/swig/include/python/gdal_python.i index fac167ee295b..3a81d2635ade 100644 --- a/swig/include/python/gdal_python.i +++ b/swig/include/python/gdal_python.i @@ -4997,6 +4997,47 @@ def InterpolateAtPoint(self, *args, **kwargs): return ret[1] %} +%feature("shadow") ComputeMinMaxLocation %{ +def ComputeMinMaxLocation(self, *args, **kwargs): + """Compute the min/max values for a band, and their location. + + If they are hit in several locations, it is not specified which one will be + returned. +. + This is a mapping of :cpp:func:`GDALRasterBand::ComputeRasterMinMaxLocation`. + + Parameters + ---------- + None + + Returns + ------- + a named tuple (min, max, minX, minY, maxX, maxY) or or ``None`` + in case of error or no valid pixel. + """ + + ret = $action(self, *args, **kwargs) + if ret[0] != CE_None: + return None + + import collections + tuple = collections.namedtuple('ComputeMinMaxLocationResult', + ['min', + 'max', + 'minX', + 'minY', + 'maxX', + 'maxY', + ]) + tuple.min = ret[1] + tuple.max = ret[2] + tuple.minX = ret[3] + tuple.minY = ret[4] + tuple.maxX = ret[5] + tuple.maxY = ret[6] + return tuple +%} + %pythoncode %{ # VSIFile: Copyright (c) 2024, Dan Baston diff --git a/swig/python/gdal-utils/osgeo_utils/samples/gdal_minmax_location.py b/swig/python/gdal-utils/osgeo_utils/samples/gdal_minmax_location.py new file mode 100644 index 000000000000..018ec00f0d29 --- /dev/null +++ b/swig/python/gdal-utils/osgeo_utils/samples/gdal_minmax_location.py @@ -0,0 +1,128 @@ +# !/usr/bin/env python3 +############################################################################### +# Project: GDAL utils +# Purpose: Get min/max location +# Author: Even Rouault +# +############################################################################### +# Copyright (c) 2024, Even Rouault +# +# SPDX-License-Identifier: MIT +############################################################################### + +import sys +import textwrap +from typing import Optional + +from osgeo import gdal, osr +from osgeo_utils.auxiliary.gdal_argparse import GDALArgumentParser, GDALScript +from osgeo_utils.auxiliary.util import PathOrDS, open_ds + + +def gdalminmaxlocation_util( + filename_or_ds: PathOrDS, + band_num: int, + open_options: Optional[dict] = None, + **kwargs, +): + ds = open_ds(filename_or_ds, open_options=open_options) + band = ds.GetRasterBand(band_num) + ret = band.ComputeMinMaxLocation() + if ret is None: + print("No valid pixels") + return 1 + gt = ds.GetGeoTransform(can_return_null=True) + if gt: + srs = ds.GetSpatialRef() + if srs: + wgs84 = osr.SpatialReference() + wgs84.SetFromUserInput("WGS84") + wgs84.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER) + ct = osr.CreateCoordinateTransformation(srs, wgs84) + georefX, georefY = gdal.ApplyGeoTransform( + gt, ret.minX + 0.5, ret.minY + 0.5 + ) + long, lat, _ = ct.TransformPoint(georefX, georefY) + print( + f"Minimum={ret.min} at (col,line)=({ret.minX},{ret.minY}), (X,Y)_georef=({georefX},{georefY}), (long,lat)_WGS84=({long:.7f},{lat:.7f})" + ) + georefX, georefY = gdal.ApplyGeoTransform( + gt, ret.maxX + 0.5, ret.maxY + 0.5 + ) + long, lat, _ = ct.TransformPoint(georefX, georefY) + print( + f"Maximum={ret.max} at (col,line)=({ret.maxX},{ret.maxY}), (X,Y)_georef=({georefX},{georefY}), (long,lat)_WGS84=({long:.7f},{lat:.7f})" + ) + else: + georefX, georefY = gdal.ApplyGeoTransform( + gt, ret.minX + 0.5, ret.minY + 0.5 + ) + print( + f"Minimum={ret.min} at (col,line)=({ret.minX},{ret.minY}), (X,Y)_georef=({georefX},{georefY})" + ) + georefX, georefY = gdal.ApplyGeoTransform( + gt, ret.maxX + 0.5, ret.maxY + 0.5 + ) + print( + f"Maximum={ret.max} at (col,line)=({ret.maxX},{ret.maxY}), (X,Y)_georef=({georefX},{georefY})" + ) + else: + print(f"Minimum={ret.min} at (col,line)=({ret.minX},{ret.minY})") + print(f"Maximum={ret.max} at (col,line)=({ret.maxX},{ret.maxY})") + + return 0 + + +class GDALMinMaxLocation(GDALScript): + def __init__(self): + super().__init__() + self.title = "Raster min/max location query tool" + self.description = textwrap.dedent( + """\ + The gdal_minmax_location utility returns the location where min/max values of a raster are hit.""" + ) + self.interactive_mode = None + + def get_parser(self, argv) -> GDALArgumentParser: + parser = self.parser + + parser.add_argument( + "-b", + dest="band_num", + metavar="band", + type=int, + default=1, + help="Selects a band to query (default: first one).", + ) + + parser.add_argument( + "-oo", + dest="open_options", + metavar="NAME=VALUE", + help="Dataset open option (format specific).", + nargs="+", + ) + + parser.add_argument( + "filename_or_ds", + metavar="filename", + type=str, + help="The source GDAL raster datasource name.", + ) + + return parser + + def augment_kwargs(self, kwargs) -> dict: + return kwargs + + def doit(self, **kwargs): + return gdalminmaxlocation_util(**kwargs) + + +def main(argv=sys.argv): + gdal.UseExceptions() + return GDALMinMaxLocation().main(argv) + + +if __name__ == "__main__": + sys.exit(main(sys.argv))