Get true maximum value of a raster

General Usage
1

Hi,

I have raster which I would like to get the maximum value from using ArchGDAL but I get only the (I believe 98th percentile (as it is done in QGIS)


using ArchGDAL
const AG = ArchGDAL

dataset = AG.readraster("raster.tif")

GDAL Dataset (Driver: GTiff/GeoTIFF)
File(s): 
  /Volumes/Black/Git_Data/node_creation/FlowR/su30m_nepal_propa_17843_nb_propagations.tif
  /Volumes/Black/Git_Data/node_creation/FlowR/su30m_nepal_propa_17843_nb_propagations.tif.aux.xml

Dataset (width x height): 79607 x 46661 (pixels)
Number of raster bands: 1
  [GA_ReadOnly] Band 1 (Gray): 79607 x 46661 (UInt32)

band = AG.getband(dataset, 1)

[GA_ReadOnly] Band 1 (Gray): 79607 x 46661 (UInt32)
    blocksize: 79607×1, nodata: nothing, units: 1.0px + 0.0
    overviews:

AG.maximum(band)
11131.0

Now if I use GeoArray

using GeoData

raster = GDALarray("raster.tif") |> GeoArray

79607×46661×1 GeoArray{UInt32,3} with dimensions: 
  X: range(-177639.86776953447, stop=618420.1322304655, length=79607)
    Projected: Ordered Regular Intervals crs: WellKnownText,
  Y: range(3.380732533405694e6, stop=2.914132533405694e6, length=46661)
    Projected: Ordered Regular Intervals crs: WellKnownText,
  Band: 1:1 Categorical: Ordered


using Statistics
m = Int(maximum(raster))

49231

49231 is the correct maximum raster value computed outside Julia for QC

Can I get the true maximum value from ArchGDAL like

band.ComputeRasterMinMax()

in python?

Thanks!

2

Hmm interesting. So the gdalgetrastermaximum function that is called here doesn’t do any computation, but gets it from the metadata if it exists. So apparently the metadata in your file was set to the percentile and not the max.

If you want to compute the max, you can call maximum(band) (Base.maximum).

3

Indeed - thank you for the solution

I manage to get the expected result by doing:


maxBand = Int(maximum(band))
4

I just happended to bump into an error this morning:


ERROR: Failed to precompile ArchGDAL [c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3] to /Users/ad/.julia/compiled/v1.6/ArchGDAL/jl_tCKQKD.
Stacktrace:
 [1] error(s::String)
   @ Base ./error.jl:33
 [2] compilecache(pkg::Base.PkgId, path::String, internal_stderr::Base.TTY, internal_stdout::Base.TTY, ignore_loaded_modules::Bool)
   @ Base ./loading.jl:1385
 [3] compilecache(pkg::Base.PkgId, path::String)
   @ Base ./loading.jl:1329
 [4] _require(pkg::Base.PkgId)
   @ Base ./loading.jl:1043
 [5] require(uuidkey::Base.PkgId)
   @ Base ./loading.jl:936
 [6] require(into::Module, mod::Symbol)
   @ Base ./loading.jl:923

Would you know if any update happened recently?

5

Yes, that is https://github.com/JuliaGeo/GDAL.jl/issues/126. If anyone has some time to edit the registry to prevent GDAL_jll releases to pick up GEOS_jll version 3.10, that’d be very welcome.

1 Like