Exporting .adf files to xyz ascii format

[Frank Warmerdam]

William,

OK, first download, and unpack:

ftp://ftp.remotesensing.org/pub/gdal/openev/OpenEV_FW_142.zip

Edit the openev_fw\setfw.bat file and fix the OPENEV_FW_DIR to point to 
whereever
you installed the tree. Dump the attached gdal2xyz.py script into the 
'openev_fw' directory
(same place as the setfw.bat file). Then in a dos window:

cd <wherever you put stuff>
setfw.bat
python gdal2xyz.py demo-data\utm.tif out.txt

This should (slowly) translate the example raster file, 
demo-data\utm.tif, into
a text file in xyz format called out.txt.

You can do the same thing for your file. Use the name of the directory 
containing
the w001001.adf as the input filename.

The gdal2xyz.py script also includes some additional commandline options for
capturing subsets of your data:

python gdal2xyz.py
Usage: gdal2xyz.py [-skip factor] [-srcwin xoff yoff width height]
[-band b] srcfile [dstfile]

The gdal2xyz.py script has also been committed to CVS, and is available 
for anyone else
who wants to do something similar. I hope this helps.

Best regards,

-- 
---------------------------------------+--------------------------------------
I set the clouds in motion - turn up | Frank Warmerdam, warmerdam at p...
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush | Geospatial Programmer for Rent



-------------- next part --------------
name="gdal2xyz.py"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
filename="gdal2xyz.py"

#!/usr/bin/env python
###############################################################################
# $Id: gdal2xyz.py,v 1.2 2002/09/04 18:11:17 warmerda Exp $
#
# Project: GDAL
# Purpose: Script to translate GDAL supported raster into XYZ ASCII
# point stream.
# Author: Frank Warmerdam, warmerdam at p...
#
###############################################################################
# Copyright (c) 2002, Frank Warmerdam <warmerdam at p...>
# 
# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Library General Public
# License as published by the Free Software Foundation; either
# version 2 of the License, or (at your option) any later version.
# 
# This library is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# Library General Public License for more details.
# 
# You should have received a copy of the GNU Library General Public
# License along with this library; if not, write to the
# Free Software Foundation, Inc., 59 Temple Place - Suite 330,
# Boston, MA 02111-1307, USA.
###############################################################################
# 
# $Log: gdal2xyz.py,v $
# Revision 1.2 2002/09/04 18:11:17 warmerda
# fixed to emit center of pixel
#
# Revision 1.1 2002/09/04 17:58:07 warmerda
# New
#
#

import gdal
import sys
import Numeric

# =============================================================================
def Usage():
print 'Usage: gdal2xyz.py [-skip factor] [-srcwin xoff yoff width height]'
print ' [-band b] srcfile [dstfile]'
print
sys.exit( 1 )

# =============================================================================
#
# Program mainline.
#

if __name__ == '__main__':

srcwin = None
skip = 1 
srcfile = None
dstfile = None
band_num = 1

# Parse command line arguments.
i = 1
while i < len(sys.argv):
arg = sys.argv[i]

if arg == '-srcwin':
srcwin = (int(sys.argv[i+1]),int(sys.argv[i+2]),
int(sys.argv[i+3]),int(sys.argv[i+4]))
i = i + 4

elif arg == '-skip':
skip = int(sys.argv[i+1])
i = i + 1

elif arg == '-band':
band_num = int(sys.argv[i+1])
i = i + 1

elif arg[0] == '-':
Usage()

elif srcfile is None:
srcfile = arg

elif dstfile is None:
dstfile = arg

else:
Usage()

i = i + 1

if srcfile is None:
Usage()

# Open source file. 
srcds = gdal.Open( srcfile )
if srcds is None:
print 'Could not open %s.' % srcfile
sys.exit( 1 )

band = srcds.GetRasterBand(band_num)
if band is None:
print 'Could not get band %d' % band_num
sys.exit( 1 )

gt = srcds.GetGeoTransform()

# Collect information on all the source files.
if srcwin is None:
srcwin = (0,0,srcds.RasterXSize,srcds.RasterYSize)

# Open the output file.
if dstfile is not None:
dst_fh = open(dstfile,'wt')
else:
dst_fh = sys.stdout

# Setup an appropriate print format.
if abs(gt[0]) < 180 and abs(gt[3]) < 180 \
and abs(srcds.RasterXSize * gt[1]) < 180 \
and abs(srcds.RasterYSize * gt[5]) < 180:
format = '%.10g %.10g %g\n'
else:
format = '%.3f %.3f %g\n'

# Loop emitting data.

for y in range(srcwin[1],srcwin[1]+srcwin[3],skip):

data = band.ReadAsArray( srcwin[0], y, srcwin[2], 1 ) 
data = Numeric.reshape( data, (srcwin[2],) )

for x_i in range(0,srcwin[2],skip):

x = x_i + srcwin[0]

geo_x = gt[0] + (x+0.5) * gt[1] + (y+0.5) * gt[2]
geo_y = gt[3] + (x+0.5) * gt[4] + (y+0.5) * gt[5]

line = format % (float(geo_x),float(geo_y),data[x_i])

dst_fh.write( line )