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

[Ujaval Gandhi]

Another approach using rasterio


import rasterio
import numpy as np


dataset = rasterio.open('P3412A.tif')
band1 = dataset.read(1)
rows, cols = np.where(band1 == np.max(band1))
lon, lat = rasterio.transform.xy(dataset.transform, row[0], col[0])
print(lat, lon)
---
Ujaval Gandhi
Spatial Thoughts
www.spatialthoughts.com
[https://mailtrack.io/l/79bb80f67e4a9e655cddd40fc6a12779870dac63?url=http%3A%2F%2Fwww.spatialthoughts.com&u=8747767&signature=198a984b1b993e73]




[data:image/gif;base64,R0lGODlhAQABAIAAAP///wAAACwAAAAAAQABAAACAkQBADs=3D]
On Wed, Oct 30, 2024 at 8:06 AM Scott via gdal-dev <gdal-dev at lists.osgeo.org>
wrote:


> It ain't pretty or efficient, but it's cheap. Here's min value. Remove
> 'r' from sort for max value:
> 
> gdal2xyz.py -csv -skipnodata source.tif /dev/stdout | grep -v done |
> sort -rnk 3,3 -t "," | tail -1
> 
> result:
> -116.9916667,36.54166667,46
> 
> On 10/29/24 17:31, Rahkonen Jukka via gdal-dev wrote:
> > 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
> [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-/
> [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 [gdal-dev at lists.osgeo.org]
> > https://lists.osgeo.org/mailman/listinfo/gdal-dev
> [https://lists.osgeo.org/mailman/listinfo/gdal-dev]
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org [gdal-dev at lists.osgeo.org]
> https://lists.osgeo.org/mailman/listinfo/gdal-dev
> [https://lists.osgeo.org/mailman/listinfo/gdal-dev]
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20241030/82a9cc80/attachment.htm>