1

I try to use gdal.BuildVRT in Python Script to calculate the maximum cell value for several overlapping raster. I know that I could use gdal.BuildVR from the QGIS library like this:

processing.run("gdal:buildvirtualraster", {
        'INPUT':list_c,
        'RESOLUTION':0,
        'SEPARATE':False,
        'PROJ_DIFFERENCE':False,
        'ADD_ALPHA':False,
        'ASSIGN_CRS':EPSG,
        'RESAMPLING':None,
        'SRC_NODATA':'255',
        'EXTRA':'',
        'OUTPUT':my_output_path})

but this doesn’t allow me to calculate the maximum value. I found this solution

but I cant make it work:

vrt_options = gdal.BuildVRTOptions(resampleAlg='max', addAlpha=False)
my_output_pat = path_to_my_output + '/' + 'output.vrt'
KRK_agg = gdal.BuildVRT(my_output_pat, list_c, options=vrt_options)
KRK_agg.FlushCache()

Any solutions?

PolyGeo
  • 65,136
  • 29
  • 109
  • 338
Miron
  • 1,231
  • 7
  • 20

1 Answers1

0

This will return the maximum value of a VRT (or a geoTiff if you use that as your input raster).

from osgeo import gdal
#The input raster
inras =r"K:\lostCreek\DEM\demgroup1\demgroup1.vrt"
#Create the gdal object...
ds = gdal.Open(inras)

def get_max_value(ds):
    '''returns the maximum values of a single band raster'''
    band = ds.GetRasterBand(1)
    stat = band.GetStatistics(True, True)
    return (band.GetNoDataValue(), stat[0], stat[1])[2]
raster_max_value = get_max_value(ds)
print(raster_max_value)
GBG
  • 9,737
  • 1
  • 15
  • 44
  • Thanks, but I am looking for something else. I edited my question. I am trying to calculate the maximum cell value of several overlapping raster. – Miron Mar 24 '22 at 08:34
  • Add them together with the raster calculator then use the code above on the final output. – GBG Mar 24 '22 at 14:38