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)


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


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

1 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)
        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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.