python function to read zmap plus asci grid format

read_zmapplus.py

def read_zmapplusgrid(inputfile,dtype=np.float64):
    """Read ZmapPlus grids
    Format is explained here http://lists.osgeo.org/pipermail/gdal-dev/2011-June/029173.html"""
    
    # read header, read until second '@', record header lines and content
    infile = open(inputfile,'r')
    comments=[]
    head=[]
    a=0 # count '@'s
    headers=0 # cound header+comment lines
    while a<2:
        #print "Reading header line", headers
        line=infile.readline()
        headers+=1
        if line[:1]=='@':
            a+=1
        if line[:1]=='!':
            comments.append(line)
        elif a==1 or (line[:1]=='@'):
            row=line.strip().strip('@').split(',')
            head.append(row)
        else:
            print "Header section does not seem to be defined by 2 @'s"
            return 0
    infile.close()
    
    # ok now process the header
    nodes_per_line                          = int(head[0][2])
    field_w                                 = head[1][0]
    null_value                              = head[1][1].strip()
    null_value2                             = head[1][2].strip()
    dec_places                              = head[1][3]
    strt_c                                  = head[1][4]
    no_rows                                 = int( head[2][0] )
    no_cols                                 = int( head[2][1] )
    minx                                    = np.float64(head[2][2])
    maxx                                    = np.float64(head[2][3])
    miny                                    = np.float64(head[2][4])
    maxy                                    = np.float64(head[2][5])
    
    # decide on the null value
    if not null_value: null_value=null_value2
    print "Read {} header lines, we have {} rows, {} cols, and data from {} to {} and {} to {}".format(headers,no_rows,no_cols, minx, maxx, miny, maxy )
    
    # load the values in, converting null value to NaN's
    #f = np.genfromtxt(inputfile, skip_header=headers,  dtype=np.float64, converters=dict([(x,conv) for x in range(nodes_per_line)]))
    infile = open(inputfile,'r')
    indata=infile.readlines()
    infile.close()
    l=0.
    f=np.zeros(no_cols*no_rows,dtype=dtype)
    n=0.
    for line in indata:
       l+=1
       #print l
       if l<=headers:
           continue
       else:
           sl=[]
           for x in line.split():
               #print x
               if "" in line:
                   print line
               if (x==null_value) or (dtype(x)==dtype(null_value)):
                   f[n]=np.nan
               else:
                   #print x
                   f[n]=dtype(x)
               n+=1
       if (n % 10000. == 0.) and (n>10000):
            print "Read node {} of {}".format(n,no_cols*no_rows)

    
    #f=np.array(f)
    print "Read {} nodes".format(len(f))#.shape[0]*f.shape[1])
    
    # reshape it, and swap axis. Beacuse columns major order from fortran.
    # Look it up, yuck
    f = f.reshape((no_cols,no_rows)).swapaxes(0,1) 
    #z = ma.masked_where(z==np.nan, z) # mask nan's
    
    # now generate x and y coords
    x  = np.linspace(minx, maxx, no_cols) # get x axis coords
    y  = np.linspace(maxy, miny, no_rows)
    x,y = np.meshgrid(x,y) # get x and y of every node
    return x,y,f