When working with large raster datasets in QGIS, especially Float64 GeoTIFFs with hundreds of millions of cells, you may want to locate and zoom into the pixel that holds the maximum value. This can be useful in terrain analysis, remote sensing, or any context where the peak value matters.

A common approach involves converting the raster into a point or polygon layer. However, this is resource-intensive and often impractical for large datasets. This tutorial introduces an efficient alternative using PyQGIS and NumPy to directly zoom to the maximum value without raster reclassification or vectorization.

Problem Description

You have a raster (e.g., Float64 TIFF) with more than 500 million cells. You already know the maximum value from the raster properties, but want to zoom into the corresponding location without converting the raster to points or polygons.

Although some workflows suggest raster reclassification followed by vector extraction, this is inefficient for large rasters. Instead, we can use Python to directly identify and zoom to the pixel with the highest value.

Prerequisites

  • QGIS installed with Python console enabled
  • numpy and gdal available (both come bundled with QGIS)
  • Sufficient memory to load the raster into a NumPy array

Step-by-Step Python Script

Run the following script in the QGIS Python Console. Replace "randrastfloat64" with the name of your raster layer.

from osgeo import gdal
import numpy as np

# Step 1: Load the raster layer by name
layer = QgsProject.instance().mapLayersByName("randrastfloat64")[0]  # Change this to your layer name

# Step 2: Open raster source with GDAL
dataset = gdal.Open(layer.source())

# Step 3: Read raster into a NumPy array
arr = dataset.ReadAsArray()

# Step 4: Find the maximum value
maxvalue = arr.max()
print(f"Pixel count: {np.multiply(*arr.shape):,}")
print(f"Max value is: {maxvalue}")

# Step 5: Locate indices (row, column) of the max value
i = np.where(arr == maxvalue)
irow, icol = i[0][0], i[1][0]

# Step 6: Extract georeference info
GT = dataset.GetGeoTransform()
upperleftx, pixelwidth, _, upperlefty, _, pixelheight = GT

# Step 7: Compute the center coordinate of the pixel
xcoord = upperleftx + icol * pixelwidth + pixelwidth / 2
ycoord = upperlefty + irow * pixelheight + pixelheight / 2

# Step 8: Zoom to the location using a bounding box
c = QgsCircle(QgsPoint(xcoord, ycoord), 20)  # Create a 20m radius circle
zoomrectangle = c.boundingBox()

iface.mapCanvas().setExtent(zoomrectangle)
iface.mapCanvas().refresh()

print("Done")

How It Works

  • The script reads the raster as a NumPy array and finds the maximum value.
  • It calculates the spatial coordinates of that pixel using the raster’s geotransform metadata.
  • A circle is created around the max-value cell, and the QGIS map canvas zooms to this extent.

Advantages

  • No raster reclassification or conversion to vector
  • Minimal processing overhead
  • Works well even with large raster datasets (given sufficient memory)

Tips

  • If your raster is too large to fit in memory, consider working with a downsampled version or use a tiled approach.
  • You can adapt this method to search for other statistics (e.g., minimum, mean, thresholds).

Conclusion

This tutorial offers an efficient and elegant way to zoom into the cell containing the maximum value in a large raster using Python within QGIS. It avoids the need for costly raster-to-vector conversion and provides an exact spatial location for immediate visualization.