Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

I have a raster dataset in ASCII format. Using Python, I need to determine the min and max values inside the dataset. I've been told that the header information is key, which contains things like number of rows/columns, cell size, and etc. My question... can't you simply skip the header information and read the entire dataset to determine the min and max values?

This is what I'm trying to do. I'm skipping the first couple of lines which contain the header information, and trying to determine the values from there on. The following is sort of what I have, but need some guidance as I am new to Python. Please excuse my mistakes.

raster_file = open('data.asc', 'r') # Open the file
data = raster_file.readlines()[4:] # Read the lines in the file, and skip the first six lines

for lines in data:
    print max(data) # Find the max value in data
    print min(data) # Find hte min value in data

Any suggestions?

share|improve this question
2  
Are you using open source or ESRI stack? –  underdark Oct 23 '11 at 9:31

5 Answers 5

You can use numpy. See the example below. A numpy masked array can be generated accounting for the no data values. See the numpy help topic for mafromtxt and genfromtxt

Below is a small ascii file with a nodata value of -999

ncols          3
nrows          3
xllcorner      0
yllcorner      0
cellsize       1
NODATA_value   -999
0 1 2
-999 4 5 
6 7 8

>>> import numpy as np
>>> ascii_file = "c:/temp/Ascii_3x3_1nodata.asc"
>>> an_array = np.mafromtxt(ascii_file, 'float', '#', None, 6, None, '-999')

>>> print an_array

[[0.0 1.0 2.0]
 [-- 4.0 5.0]  
 [6.0 7.0 8.0]]

>>>

from there it is simply a matter to determine the statistics you want

>>> print an_array.min()
0.0
>>> print an_array.max()
8.0
>>> print an_array.mean()
4.125
>>> 
share|improve this answer
    
Thank you Dan. I will give this a try. Is there an alternative way... maybe without the numpy module? –  kolor Oct 24 '11 at 6:04

You want raster data statistics.
See what you are doing in the gui first (for homework.)

Then you can use a python window or a script.

import arcpy
arcpy.CalculateStatistics_management("c:/data/image.tif", "4", "6", "0;255;21")
share|improve this answer
    
Once you calculate statistics, you can always access the statistics through the raster object property too. e.g. r = arcpy.Raster("c:/data/image.tif"), r.mean, r.minimum, r.maximum –  blord-castillo Oct 23 '11 at 14:41
    
@blord-castillo Cool! Didn't know that. Thanks for the tip :) –  kolor Oct 24 '11 at 6:06

edit - now tested and working -
edit - improved code

import sys

class Ascii_file(object):
    def __init__(self,file):
        self.raster_file = open(file, 'r') # Open the file
        self.max=sys.float_info.min
        self.min=sys.float_info.max
    def __minmax(self,value):
        if value>self.max:self.max=value
        if value<self.min:self.min=value
    def getMinMax(self):
        data = self.raster_file.readlines()
        data_values=data[6:]
        nodata=float(data[5].split()[1])
        for line in data_values:
            values=line.split(" ")
            for value in values:
                value=float(value)
                if value==nodata:continue
                else: self.__minmax(value)
        return self.min, self.max

if __name__=="__main__":
    myfile = Ascii_file('data.asc')
    print myfile.getMinMax()
share|improve this answer
    
This is sort of what I was trying earlier, but I keep getting errors when I use the split method: AttributeError: 'list' object has no attribute 'split' –  kolor Oct 25 '11 at 19:47
    
I feel like the line data = raster_file.readlines()[4:] doesn't actually work when it comes to specifying the range. I fixed the error which I was having in the previous comment. This was done by adding num = data[7] in the 3rd line. It was then split using values = num.split() and was able to find the max/min, but for that particular line only. How can I find the max/min from the whole document? –  kolor Oct 26 '11 at 0:42
    
oh, my mistake, "data" is a list, "lines" is the string. I've edited the code... I tested it with an asc file. Just copy and paste, pay atention to the indentation. –  Pablo Oct 26 '11 at 12:11
1  
You can drop the if check==True block by initializing your min/max values. You'll want to initialize min to sys.float_info.max and max to sys.float_info.min. –  Sasa Ivetic Oct 26 '11 at 14:07
2  
You have to initialize max to sys.float_info.min, and min to sys.float_info.max. That you your initial min will be the largest possible value, and any value you compare to it will be smaller, and thus become the new min. Same goes with your max value, it will be the smallest possible value, and any value you compare to it will be larger, and thus the new max. –  Sasa Ivetic Oct 26 '11 at 14:41

If you don't want to use numpy (and you really should, it's perfect for this sort of thing), then you will need to:

  • initialise your maximum variable to a very large negative number and your minimum variable to a very large positive number
  • split each line to get a list of strings and use list comprehension to convert it to a list of floats
  • finally use something like maximum = max(maximum, max(myfloatlist)) and an equivalent for the minimum value.
share|improve this answer

I just did this the other day. I used arcpy.RasterToNumPyArray, converted the numpy array to a list, then iterated through my list via a list comprehension to find the min and max values. But I have to admit, I like some of the other solutions here better, as mine relies on being on the ESRI stack.

import arcpy
import numpy
myArray = arcpy.RasterToNumPyArray(r"D:\NED_93512417\NED_93512417_3DEM_RPRJ.TIF")
p = myArray.tolist()

max_elev = max([item for sublist in p for item in sublist])
min_elev = min([item for sublist in p for item in sublist])
share|improve this answer
    
isn't myArray.min()/myArray.max() simpler/faster? –  Mike T Oct 24 '11 at 17:13
    
@Chad, if you already have the numpy array, then there is no need to convert to a list, simply use the min(), max() etc functions in my thread above. As you also note, no implied access to Arcpy was indicated. –  Dan Patterson Oct 24 '11 at 17:15

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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