Find coordinates corresponding to max and min values using Python

Question

I'm using GDAL and numpy so far in my Python program. I am able to get the max and min values of a band within a .grib2 file with:

def getDataStats(grib2Filepath, bandNr):
    local_minimum = None
    local_maximum = None

    file = gdal.Open(grib2Filepath)
    srcband = file.GetRasterBand(bandNr)
    fileStats = srcband.GetStatistics(True, True)
    local_minimum = fileStats[0]
    local_maximum = fileStats[1]

I am trying to get the coordinates that correspond to the maximum and minimum values. I've tried:

data_array = srcband.ReadAsArray().astype(float)
(y_index, x_index) = np.nonzero(data_array == local_minimum)

But

print(f"min: {local_minimum}, max: {local_maximum}, x: {x_index}, y: {y_index}")

serves

min: -51.380010986328, max: 36.499993896484, x: [], y: []

Answer 1

from osgeo import gdal
import numpy as np

ds = gdal.Open(r"/home/bera/Desktop/gistest/dem.tif")
arr = ds.ReadAsArray()

def findit(dataset, array, value):
    """Prints the coordinates of where an input value
       is found in an array"""
    if value in array:
        i = np.where(array == value) #Find the indices of the value
        irow, icol = i[0][0], i[1][0]
        GT = dataset.GetGeoTransform()
        upperleftx, pixelwidth, upperlefty, pixelheight = GT[0], GT[1], GT[3], GT[5]
        xcoord = upperleftx+icol*pixelwidth+pixelwidth/2
        ycoord = upperlefty+irow*pixelheight+pixelheight/2
        print(xcoord, ycoord)
    else:
        print(f"{value} is not in array")

findit(ds, arr, arr.max())
#429487.00008182594 6178720.219712928
findit(ds, arr, arr.min())
#409083.7854456485 6173480.200701521
findit(ds, arr, 9999)
#9999 is not in array