4

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.

5

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.

  • I prefer your solution because it can be done in a few lines and my code is already very long. I found only one problem, that is, if I want the minimum value your code will give me -9999 which theoretically is right but in practice it represents a NODATA. How to solve this problem with a few lines of code? – Domenico Fuoco Aug 1 '17 at 17:51
  • @mapr Maybe you can help me because you've had a similar problem – Domenico Fuoco Aug 1 '17 at 18:59
  • 1
    @DomenicoFuoco Good spot. I've edited my answer to include a method that worked for me in regards to masking NODATA. – Ali Aug 1 '17 at 19:01
  • I'm sorry but whit this metod the results of 'min = x.min()' is "masked" – Domenico Fuoco Aug 1 '17 at 19:27
  • 1
    @DomenicoFuoco Sorry, that was a typo, fixed now – Ali Aug 1 '17 at 19:32
4

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] )
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)
3

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))
  • Your solution is more elegant and functional, thank you! – Domenico Fuoco Aug 2 '17 at 9:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.