Max and Min value from an ascii grid

Question

I'm working in Qgis with python console. I have a file .asc whose structure is:

Ncols 49

Nrows 52

Xllcorner 2611187.75

Yllcorner 4336914.199

Cellsize 1000

NODATA_value -9999

0.5081 0.679 -9999 ...

I have to find the Max and the Min value in the file and I tried:

GridFile = open('filepath \ file.asc', 'r')
MaxValue = max (GridFile)
MinValue = min (GridFile)

But the procedure returns the following error

ValueError: max () is an empty sequence
ValueError: min () is an empty sequence

What I understand is that I have to skip the first 6 lines and find the maximum and the minimum of successive values, but I do not know how to do it.

Answer 1

One quick approach is to use Numpy's genfromtext function. It allows you skip lines and specify a nodata value. Based on the asc file structure (pre-edit as @Joseph has done):

import numpy as np

asc_file = r".../yourfile.asc"

x = np.genfromtxt(asc_file,skip_header=12,missing_values='-9999')

print x.max()

This returns the value 0.7315

---

EDIT

In regards to the nodata value not actually being ignored I found that using the usemask parameter worked better. It returns a masked array ensuring that subsequent operations will ignore the nodata values, or missing_values. Alternatively, instead of genfromtxt you can use mafromtext which accepts the same arguments but sets usemask=True by default. The following worked for me:

# usemask=True so masked array is returned
x = np.genfromtxt(asc_file,skip_header=12,missing_values='-9999',usemask=True)

# Equivalent using mafromtext
x = np.mafromtxt(asc_file,skip_header=12,missing_values='-9999')

# Now any calls to min/max will not include the NODATA values
print x.min()
print x.max()

The Numpy docs have a page which explains the use of missing_values and other parameters in more detail.

---

Of course if you require more functionality in terms of raster processing and interrogation then GDAL is the way to go, as @iant and @Luke have mentioned.

Answer 2

The obvious (to me) way to do this is to use GDAL. So looking in the GDAL python cookbook gives:

from osgeo import gdal
gtif = gdal.Open( "INPUT.asc" )
print gtif.GetMetadata()

gets you started. Further down the page you will find:

from osgeo import gdal
import sys
gdal.UseExceptions()

def Usage():
    print("""
    $ getrasterband.py [ band number ] input-raster
    """)
    sys.exit(1)

def main( band_num, input_file ):
    src_ds = gdal.Open( input_file )
    if src_ds is None:
        print 'Unable to open %s' % input_file
        sys.exit(1)

    try:
        srcband = src_ds.GetRasterBand(band_num)
    except RuntimeError, e:
        print 'No band %i found' % band_num
        print e
        sys.exit(1)


    print "[ NO DATA VALUE ] = ", srcband.GetNoDataValue()
    print "[ MIN ] = ", srcband.GetMinimum()
    print "[ MAX ] = ", srcband.GetMaximum()
    print "[ SCALE ] = ", srcband.GetScale()
    print "[ UNIT TYPE ] = ", srcband.GetUnitType()
    ctable = srcband.GetColorTable()

    if ctable is None:
        print 'No ColorTable found'
        sys.exit(1)

    print "[ COLOR TABLE COUNT ] = ", ctable.GetCount()
    for i in range( 0, ctable.GetCount() ):
        entry = ctable.GetColorEntry( i )
        if not entry:
            continue
        print "[ COLOR ENTRY RGB ] = ", ctable.GetColorEntryAsRGB( i, entry )

if __name__ == '__main__':

    if len( sys.argv ) < 3:
        print """
        [ ERROR ] you must supply at least two arguments:
        1) the band number to retrieve and 2) input raster
        """
        Usage()

    main( int(sys.argv[1]), sys.argv[2] )

Answer 3

Not the most elegant of solutions but you could try using the following code which:

• Finds the line beginning with a value
• Appends the line into a list
• Create a final list containing all values
• Returns the maximum value

Using the structure you posted (before it was edited), the maximum value was found to be 0.7315.

---

Here is the code:

# Create list to contain all values
list_of_values = []
# Open your file
with open("path/to/file.asc") as f:
    # Return each line as strings in a list
    mylist = f.read().splitlines()
    # Remove newlines marked by ''
    new_list = [x for x in mylist if x != '']
    # Iterate through the list  
    for x in new_list:
        # Find lines beginning with a value
        if x[0].isdigit():
            # Append this line of values to list
            list_of_values.append(x.split())

# Convert any nested lists to a single list
final_list = [values for sublist in list_of_values for values in sublist]
# Print maximum value
print max(final_list)

Answer 4

As noted by @iant

The obvious (to me) way to do this is to use GDAL

Some simple GDAL code would be something like:

from osgeo import gdal
asc = "input.asc"

ds = gdal.Open(asc)
rb = ds.GetRasterBand(1)
min, max, mean, stdv = rb.ComputeStatistics(False)
print((min,max))