Skip to content

Instantly share code, notes, and snippets.

@wassname
Created June 8, 2014 21:45
Show Gist options
  • Save wassname/526d5fde3f3cbeb67da8 to your computer and use it in GitHub Desktop.
Save wassname/526d5fde3f3cbeb67da8 to your computer and use it in GitHub Desktop.
python function to read zmap plus asci grid format
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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment