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
33 import os.path
34 import sys
35
36 from osgeo import gdal
37 from osgeo import ogr
38
39
41 print("""
42 gdal_polygonize [-8] [-nomask] [-mask filename] raster_file [-b band|mask]
43 [-q] [-f ogr_format] out_file [layer] [fieldname]
44 """)
45 sys.exit(1)
46
47
51
52
53 -def GetExtension(filename):
54 ext = os.path.splitext(filename)[1]
55 if ext.startswith('.'):
56 ext = ext[1:]
57 return ext
58
59
76
77
79 drv_list = GetOutputDriversFor(filename)
80 ext = GetExtension(filename)
81 if not drv_list:
82 if not ext:
83 return 'ESRI Shapefile'
84 else:
85 raise Exception("Cannot guess driver for %s" % filename)
86 elif len(drv_list) > 1:
87 print("Several drivers matching %s extension. Using %s" % (ext if ext else '', drv_list[0]))
88 return drv_list[0]
89
90
92 frmt = None
93 options = []
94 quiet_flag = 0
95 src_filename = None
96 src_band_n = 1
97
98 dst_filename = None
99 dst_layername = None
100 dst_fieldname = None
101 dst_field = -1
102
103 mask = 'default'
104
105 gdal.AllRegister()
106 argv = gdal.GeneralCmdLineProcessor(argv)
107 if argv is None:
108 sys.exit(0)
109
110
111 i = 1
112 while i < len(argv):
113 arg = argv[i]
114
115 if arg == '-f' or arg == '-of':
116 i = i + 1
117 frmt = argv[i]
118
119 elif arg == '-q' or arg == '-quiet':
120 quiet_flag = 1
121
122 elif arg == '-8':
123 options.append('8CONNECTED=8')
124
125 elif arg == '-nomask':
126 mask = 'none'
127
128 elif arg == '-mask':
129 i = i + 1
130 mask = argv[i]
131
132 elif arg == '-b':
133 i = i + 1
134 if argv[i].startswith('mask'):
135 src_band_n = argv[i]
136 else:
137 src_band_n = int(argv[i])
138
139 elif src_filename is None:
140 src_filename = argv[i]
141
142 elif dst_filename is None:
143 dst_filename = argv[i]
144
145 elif dst_layername is None:
146 dst_layername = argv[i]
147
148 elif dst_fieldname is None:
149 dst_fieldname = argv[i]
150
151 else:
152 Usage()
153
154 i = i + 1
155
156 if src_filename is None or dst_filename is None:
157 Usage()
158
159 if frmt is None:
160 frmt = GetOutputDriverFor(dst_filename)
161
162 if dst_layername is None:
163 dst_layername = 'out'
164
165
166
167
168 try:
169 gdal.Polygonize
170 except AttributeError:
171 print('')
172 print('gdal.Polygonize() not available. You are likely using "old gen"')
173 print('bindings or an older version of the next gen bindings.')
174 print('')
175 sys.exit(1)
176
177
178
179
180
181 src_ds = gdal.Open(src_filename)
182
183 if src_ds is None:
184 print('Unable to open %s' % src_filename)
185 sys.exit(1)
186
187 if src_band_n == 'mask':
188 srcband = src_ds.GetRasterBand(1).GetMaskBand()
189
190 options.append('DATASET_FOR_GEOREF=' + src_filename)
191 elif isinstance(src_band_n, str) and src_band_n.startswith('mask,'):
192 srcband = src_ds.GetRasterBand(int(src_band_n[len('mask,'):])).GetMaskBand()
193
194 options.append('DATASET_FOR_GEOREF=' + src_filename)
195 else:
196 srcband = src_ds.GetRasterBand(src_band_n)
197
198 if mask == 'default':
199 maskband = srcband.GetMaskBand()
200 elif mask == 'none':
201 maskband = None
202 else:
203 mask_ds = gdal.Open(mask)
204 maskband = mask_ds.GetRasterBand(1)
205
206
207
208
209
210 try:
211 gdal.PushErrorHandler('CPLQuietErrorHandler')
212 dst_ds = ogr.Open(dst_filename, update=1)
213 gdal.PopErrorHandler()
214 except:
215 dst_ds = None
216
217
218
219
220 if dst_ds is None:
221 drv = ogr.GetDriverByName(frmt)
222 if not quiet_flag:
223 print('Creating output %s of format %s.' % (dst_filename, frmt))
224 dst_ds = drv.CreateDataSource(dst_filename)
225
226
227
228
229 try:
230 dst_layer = dst_ds.GetLayerByName(dst_layername)
231 except:
232 dst_layer = None
233
234 if dst_layer is None:
235
236 srs = src_ds.GetSpatialRef()
237 dst_layer = dst_ds.CreateLayer(dst_layername, geom_type=ogr.wkbPolygon, srs=srs)
238
239 if dst_fieldname is None:
240 dst_fieldname = 'DN'
241
242 fd = ogr.FieldDefn(dst_fieldname, ogr.OFTInteger)
243 dst_layer.CreateField(fd)
244 dst_field = 0
245 else:
246 if dst_fieldname is not None:
247 dst_field = dst_layer.GetLayerDefn().GetFieldIndex(dst_fieldname)
248 if dst_field < 0:
249 print("Warning: cannot find field '%s' in layer '%s'" % (dst_fieldname, dst_layername))
250
251
252
253
254
255 if quiet_flag:
256 prog_func = None
257 else:
258 prog_func = gdal.TermProgress_nocb
259
260 result = gdal.Polygonize(srcband, maskband, dst_layer, dst_field, options,
261 callback=prog_func)
262
263 srcband = None
264 src_ds = None
265 dst_ds = None
266 mask_ds = None
267
268 return result
269
270
271 if __name__ == '__main__':
272 sys.exit(main(sys.argv))
273