1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
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
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
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
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
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
121 if srcwin is None:
122 srcwin = (0, 0, srcds.RasterXSize, srcds.RasterYSize)
123
124
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
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
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