Package osgeo :: Package utils :: Module gdal2xyz |
[hide private]
[frames] | no frames]

# Source Code for Module osgeo.utils.gdal2xyz

      1  #!/usr/bin/env python3
      2  ###############################################################################
      3  # $Id$
      4  #
      5  # Project:  GDAL
      6  # Purpose:  Script to translate GDAL supported raster into XYZ ASCII
      7  #           point stream.
      8  # Author:   Frank Warmerdam, warmerdam@pobox.com
      9  #
     10  ###############################################################################
     11  # Copyright (c) 2002, Frank Warmerdam <warmerdam@pobox.com>
     12  #
     13  # Permission is hereby granted, free of charge, to any person obtaining a
     14  # copy of this software and associated documentation files (the "Software"),
     15  # to deal in the Software without restriction, including without limitation
     16  # the rights to use, copy, modify, merge, publish, distribute, sublicense,
     17  # and/or sell copies of the Software, and to permit persons to whom the
     18  # Software is furnished to do so, subject to the following conditions:
     19  #
     20  # The above copyright notice and this permission notice shall be included
     21  # in all copies or substantial portions of the Software.
     22  #
     23  # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
     24  # OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
     25  # FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
     26  # THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
     27  # LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
     28  # FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
     29  # DEALINGS IN THE SOFTWARE.
     30  ###############################################################################
     31
     32  import sys
     33
     34  from osgeo import gdal
     35
     36  try:
     37      import numpy as Numeric
     38  except ImportError:
     39      import Numeric
     40
     41  # =============================================================================
     42
     43
     44 -def Usage():
     45      print('Usage: gdal2xyz.py [-skip factor] [-srcwin xoff yoff width height]')
     46      print('                   [-band b] [-csv] srcfile [dstfile]')
     47      print('')
     48      sys.exit(1)
     49
     50
     51 -def main(argv):
     52      srcwin = None
     53      skip = 1
     54      srcfile = None
     55      dstfile = None
     56      band_nums = []
     57      delim = ' '
     58
     59      gdal.AllRegister()
     60      argv = gdal.GeneralCmdLineProcessor(argv)
     61      if argv is None:
     62          sys.exit(0)
     63
     64      # Parse command line arguments.
     65      i = 1
     66      while i < len(argv):
     67          arg = argv[i]
     68
     69          if arg == '-srcwin':
     70              srcwin = (int(argv[i + 1]), int(argv[i + 2]),
     71                        int(argv[i + 3]), int(argv[i + 4]))
     72              i = i + 4
     73
     74          elif arg == '-skip':
     75              skip = int(argv[i + 1])
     76              i = i + 1
     77
     78          elif arg == '-band':
     79              band_nums.append(int(argv[i + 1]))
     80              i = i + 1
     81
     82          elif arg == '-csv':
     83              delim = ','
     84
     85          elif arg[0] == '-':
     86              Usage()
     87
     88          elif srcfile is None:
     89              srcfile = arg
     90
     91          elif dstfile is None:
     92              dstfile = arg
     93
     94          else:
     95              Usage()
     96
     97          i = i + 1
     98
     99      if srcfile is None:
    100          Usage()
    101
    102      if band_nums == []:
    103          band_nums = [1]
    104      # Open source file.
    105      srcds = gdal.Open(srcfile)
    106      if srcds is None:
    107          print('Could not open %s.' % srcfile)
    108          sys.exit(1)
    109
    110      bands = []
    111      for band_num in band_nums:
    112          band = srcds.GetRasterBand(band_num)
    113          if band is None:
    114              print('Could not get band %d' % band_num)
    115              sys.exit(1)
    116          bands.append(band)
    117
    118      gt = srcds.GetGeoTransform()
    119
    120      # Collect information on all the source files.
    121      if srcwin is None:
    122          srcwin = (0, 0, srcds.RasterXSize, srcds.RasterYSize)
    123
    124      # Open the output file.
    125      if dstfile is not None:
    126          dst_fh = open(dstfile, 'wt')
    127      else:
    128          dst_fh = sys.stdout
    129
    130      dt = srcds.GetRasterBand(1).DataType
    131      if dt == gdal.GDT_Int32 or dt == gdal.GDT_UInt32:
    132          band_format = (("%d" + delim) * len(bands)).rstrip(delim) + '\n'
    133      else:
    134          band_format = (("%g" + delim) * len(bands)).rstrip(delim) + '\n'
    135
    136      # Setup an appropriate print format.
    137      if abs(gt[0]) < 180 and abs(gt[3]) < 180 \
    138         and abs(srcds.RasterXSize * gt[1]) < 180 \
    139         and abs(srcds.RasterYSize * gt[5]) < 180:
    140          frmt = '%.10g' + delim + '%.10g' + delim + '%s'
    141      else:
    142          frmt = '%.3f' + delim + '%.3f' + delim + '%s'
    143
    144      # Loop emitting data.
    145
    146      for y in range(srcwin[1], srcwin[1] + srcwin[3], skip):
    147
    148          data = []
    149          for band in bands:
    150
    151              band_data = band.ReadAsArray(srcwin[0], y, srcwin[2], 1)
    152              band_data = Numeric.reshape(band_data, (srcwin[2],))
    153              data.append(band_data)
    154
    155          for x_i in range(0, srcwin[2], skip):
    156
    157              x = x_i + srcwin[0]
    158
    159              geo_x = gt[0] + (x + 0.5) * gt[1] + (y + 0.5) * gt[2]
    160              geo_y = gt[3] + (x + 0.5) * gt[4] + (y + 0.5) * gt[5]
    161
    162              x_i_data = []
    163              for i in range(len(bands)):
    164                  x_i_data.append(data[i][x_i])
    165
    166              band_str = band_format % tuple(x_i_data)
    167
    168              line = frmt % (float(geo_x), float(geo_y), band_str)
    169
    170              dst_fh.write(line)
    171
    172
    173  if __name__ == '__main__':
    174      sys.exit(main(sys.argv))
    175