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 os
33 import os.path
34 import sys
35 from osgeo import gdal
36
37
41
42
43 -def GetExtension(filename):
44 ext = os.path.splitext(filename)[1]
45 if ext.startswith('.'):
46 ext = ext[1:]
47 return ext
48
49
72
73
75 drv_list = GetOutputDriversFor(filename)
76 ext = GetExtension(filename)
77 if not drv_list:
78 if not ext:
79 return 'GTiff'
80 else:
81 raise Exception("Cannot guess driver for %s" % filename)
82 elif len(drv_list) > 1:
83 print("Several drivers matching %s extension. Using %s" % (ext if ext else '', drv_list[0]))
84 return drv_list[0]
85
86
88 print('Usage: gdal_pansharpen [--help-general] pan_dataset {spectral_dataset[,band=num]}+ out_dataset')
89 print(' [-of format] [-b band]* [-w weight]*')
90 print(' [-r {nearest,bilinear,cubic,cubicspline,lanczos,average}]')
91 print(' [-threads {ALL_CPUS|number}] [-bitdepth val] [-nodata val]')
92 print(' [-spat_adjust {union,intersection,none,nonewithoutwarning}]')
93 print(' [-verbose_vrt] [-co NAME=VALUE]* [-q]')
94 print('')
95 print('Create a dataset resulting from a pansharpening operation.')
96 return -1
97
98
100
101 argv = gdal.GeneralCmdLineProcessor(argv)
102 if argv is None:
103 return -1
104
105 pan_name = None
106 last_name = None
107 spectral_ds = []
108 spectral_bands = []
109 out_name = None
110 bands = []
111 weights = []
112 frmt = None
113 creation_options = []
114 callback = gdal.TermProgress_nocb
115 resampling = None
116 spat_adjust = None
117 verbose_vrt = False
118 num_threads = None
119 bitdepth = None
120 nodata = None
121
122 i = 1
123 argc = len(argv)
124 while i < argc:
125 if (argv[i] == '-of' or argv[i] == '-f') and i < len(argv) - 1:
126 frmt = argv[i + 1]
127 i = i + 1
128 elif argv[i] == '-r' and i < len(argv) - 1:
129 resampling = argv[i + 1]
130 i = i + 1
131 elif argv[i] == '-spat_adjust' and i < len(argv) - 1:
132 spat_adjust = argv[i + 1]
133 i = i + 1
134 elif argv[i] == '-b' and i < len(argv) - 1:
135 bands.append(int(argv[i + 1]))
136 i = i + 1
137 elif argv[i] == '-w' and i < len(argv) - 1:
138 weights.append(float(argv[i + 1]))
139 i = i + 1
140 elif argv[i] == '-co' and i < len(argv) - 1:
141 creation_options.append(argv[i + 1])
142 i = i + 1
143 elif argv[i] == '-threads' and i < len(argv) - 1:
144 num_threads = argv[i + 1]
145 i = i + 1
146 elif argv[i] == '-bitdepth' and i < len(argv) - 1:
147 bitdepth = argv[i + 1]
148 i = i + 1
149 elif argv[i] == '-nodata' and i < len(argv) - 1:
150 nodata = argv[i + 1]
151 i = i + 1
152 elif argv[i] == '-q':
153 callback = None
154 elif argv[i] == '-verbose_vrt':
155 verbose_vrt = True
156 elif argv[i][0] == '-':
157 sys.stderr.write('Unrecognized option : %s\n' % argv[i])
158 return Usage()
159 elif pan_name is None:
160 pan_name = argv[i]
161 pan_ds = gdal.Open(pan_name)
162 if pan_ds is None:
163 return 1
164 else:
165 if last_name is not None:
166 pos = last_name.find(',band=')
167 if pos > 0:
168 spectral_name = last_name[0:pos]
169 ds = gdal.Open(spectral_name)
170 if ds is None:
171 return 1
172 band_num = int(last_name[pos + len(',band='):])
173 band = ds.GetRasterBand(band_num)
174 spectral_ds.append(ds)
175 spectral_bands.append(band)
176 else:
177 spectral_name = last_name
178 ds = gdal.Open(spectral_name)
179 if ds is None:
180 return 1
181 for j in range(ds.RasterCount):
182 spectral_ds.append(ds)
183 spectral_bands.append(ds.GetRasterBand(j + 1))
184
185 last_name = argv[i]
186
187 i = i + 1
188
189 if pan_name is None or not spectral_bands:
190 return Usage()
191 out_name = last_name
192
193 if frmt is None:
194 frmt = GetOutputDriverFor(out_name)
195
196 if not bands:
197 bands = [j + 1 for j in range(len(spectral_bands))]
198 else:
199 for band in bands:
200 if band < 0 or band > len(spectral_bands):
201 print('Invalid band number in -b: %d' % band)
202 return 1
203
204 if weights and len(weights) != len(spectral_bands):
205 print('There must be as many -w values specified as input spectral bands')
206 return 1
207
208 vrt_xml = """<VRTDataset subClass="VRTPansharpenedDataset">\n"""
209 if bands != [j + 1 for j in range(len(spectral_bands))]:
210 for i, band in enumerate(bands):
211 sband = spectral_bands[band - 1]
212 datatype = gdal.GetDataTypeName(sband.DataType)
213 colorname = gdal.GetColorInterpretationName(sband.GetColorInterpretation())
214 vrt_xml += """ <VRTRasterBand dataType="%s" band="%d" subClass="VRTPansharpenedRasterBand">
215 <ColorInterp>%s</ColorInterp>
216 </VRTRasterBand>\n""" % (datatype, i + 1, colorname)
217
218 vrt_xml += """ <PansharpeningOptions>\n"""
219
220 if weights:
221 vrt_xml += """ <AlgorithmOptions>\n"""
222 vrt_xml += """ <Weights>"""
223 for i, weight in enumerate(weights):
224 if i > 0:
225 vrt_xml += ","
226 vrt_xml += "%.16g" % weight
227 vrt_xml += "</Weights>\n"
228 vrt_xml += """ </AlgorithmOptions>\n"""
229
230 if resampling is not None:
231 vrt_xml += ' <Resampling>%s</Resampling>\n' % resampling
232
233 if num_threads is not None:
234 vrt_xml += ' <NumThreads>%s</NumThreads>\n' % num_threads
235
236 if bitdepth is not None:
237 vrt_xml += ' <BitDepth>%s</BitDepth>\n' % bitdepth
238
239 if nodata is not None:
240 vrt_xml += ' <NoData>%s</NoData>\n' % nodata
241
242 if spat_adjust is not None:
243 vrt_xml += ' <SpatialExtentAdjustment>%s</SpatialExtentAdjustment>\n' % spat_adjust
244
245 pan_relative = '0'
246 if frmt.upper() == 'VRT':
247 if not os.path.isabs(pan_name):
248 pan_relative = '1'
249 pan_name = os.path.relpath(pan_name, os.path.dirname(out_name))
250
251 vrt_xml += """ <PanchroBand>
252 <SourceFilename relativeToVRT="%s">%s</SourceFilename>
253 <SourceBand>1</SourceBand>
254 </PanchroBand>\n""" % (pan_relative, pan_name)
255
256 for i, sband in enumerate(spectral_bands):
257 dstband = ''
258 for j, band in enumerate(bands):
259 if i + 1 == band:
260 dstband = ' dstBand="%d"' % (j + 1)
261 break
262
263 ms_relative = '0'
264 ms_name = spectral_ds[i].GetDescription()
265 if frmt.upper() == 'VRT':
266 if not os.path.isabs(ms_name):
267 ms_relative = '1'
268 ms_name = os.path.relpath(ms_name, os.path.dirname(out_name))
269
270 vrt_xml += """ <SpectralBand%s>
271 <SourceFilename relativeToVRT="%s">%s</SourceFilename>
272 <SourceBand>%d</SourceBand>
273 </SpectralBand>\n""" % (dstband, ms_relative, ms_name, sband.GetBand())
274
275 vrt_xml += """ </PansharpeningOptions>\n"""
276 vrt_xml += """</VRTDataset>\n"""
277
278 if frmt.upper() == 'VRT':
279 f = gdal.VSIFOpenL(out_name, 'wb')
280 if f is None:
281 print('Cannot create %s' % out_name)
282 return 1
283 gdal.VSIFWriteL(vrt_xml, 1, len(vrt_xml), f)
284 gdal.VSIFCloseL(f)
285 if verbose_vrt:
286 vrt_ds = gdal.Open(out_name, gdal.GA_Update)
287 vrt_ds.SetMetadata(vrt_ds.GetMetadata())
288 else:
289 vrt_ds = gdal.Open(out_name)
290 if vrt_ds is None:
291 return 1
292
293 return 0
294
295 vrt_ds = gdal.Open(vrt_xml)
296 out_ds = gdal.GetDriverByName(frmt).CreateCopy(out_name, vrt_ds, 0, creation_options, callback=callback)
297 if out_ds is None:
298 return 1
299 return 0
300
301
303 return gdal_pansharpen(argv)
304
305
306 if __name__ == '__main__':
307 sys.exit(main(sys.argv))
308