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