[gdal-dev] How to locate where the raster min and max values are?

[Even Rouault]

if you are just interested in any point where the minimum is reached:

import numpy as np
from osgeo import gdal
ds = gdal.Open('P3412A.tif', gdal.GA_ReadOnly)
rb = ds.GetRasterBand(1)
img_array = rb.ReadAsArray()
vmin = img_array.min()
(line, col) = np.where(img_array==vmin)
line = line[0]
col = col[0]
print(f"col={col}, line={line}")
gt = ds.GetGeoTransform()
X = gt[0] + gt[1] * col + gt[2] * line
Y = gt[3] + gt[4] * col + gt[5] * line
print(f"X={X}, Y={Y}")

Le 30/10/2024 à 01:31, Rahkonen Jukka via gdal-dev a écrit :
>
> Hi,
>
> I would like to know the georeferenced coordinates of the min and max 
> values of a DEM file. Even better if I could forward them into a 
> vector file. If the minimum or maximum happens to be on a flat area 
> like seabed I would be happy with the first pixel at the moment.
>
> By copy-pasting from How do I open geotiff images with GDAL in Python? 
> - Stack Overflow 
> <https://stackoverflow.com/questions/41996079/how-do-i-open-geotiff-images-with-gdal-in-python> 
> and How to find the indexes of the minimum or maximum value(s) in a 
> matrix using python ? 
> <https://en.moonbooks.org/Articles/How-to-find-the-indexes-of-the-minimum-or-maximum-values-in-a-matrix-using-python-/> 
> I think I managed to get the correct points as numpy indexes
>
> >>> import numpy as np
>
> >>> from osgeo import gdal
>
> >>> ds = gdal.Open('P3412A.tif', gdal.GA_ReadOnly)
>
> >>> rb = ds.GetRasterBand(1)
>
> >>> img_array = rb.ReadAsArray()
>
> >>> vmin = img_array.min()
>
> >>> vmax = img_array.max()
>
> >>> vmin
>
> -0.929
>
> >>> vmax
>
> 17.246
>
> >>>
>
> >>> np.where(img_array==vmin)
>
> (array([1504], dtype=intg64), array([1189], dtype=int64))
>
> >>> np.where(img_array==vmax)
>
> (array([1545], dtype=int64), array([2423], dtype=int64))
>
> >>>
>
> But now I have no idea about how to get the georeferenced coordinates.
>
> The task feels rather simple and I was sure that someone has already 
> made an utility or a QGIS plugin, but all I have found yet is for R. I 
> was thinking that perhaps some of the gdaldem modes could be misused 
> for this purpose, but I believe they cannot. For QGIS I found advice 
> to use an obvious but  clumsy method of polygonising the raster and 
> finding the extremes from the vector data. And one OpenJUMP developer 
> took the challenge and wrote a prototype with Java but it is not 
> complete yet.
>
> -Jukka Rahkonen-
>
>
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/gdal-dev

-- 
http://www.spatialys.com
My software is free, but my time generally not.
Butcher of all kinds of standards, open or closed formats. At the end, this is just about bytes.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20241030/ebc66dd0/attachment.htm>