[gdal-dev] Problem with gdal_pansharpen tool
Even Rouault
even.rouault at spatialys.com
Fri Apr 22 08:54:04 PDT 2016
Le vendredi 22 avril 2016 17:37:19, Lorenzo Bernardi a écrit :
> Hi,
> i'm trying to perfom a pansharpening process using the gdal_pansharpen
> tool. I used the following code:
>
> #!/usr/bin/env python
> # -*- coding: utf-8 -*-
>
>
> import sys, os
> from osgeo import gdal
>
> def Usage():
> print('Usage: gdal_pansharpen [--help-general] pan_dataset
> {spectral_dataset[,band=num]}+ out_dataset')
> print(' [-of format] [-b band]* [-w weight]*')
> print(' [-r
> {nearest,bilinear,cubic,cubicspline,lanczos,average}]')
> print(' [-threads {ALL_CPUS|number}]
> [-bitdepth val] [-nodata val]')
> print(' [-spat_adjust
> {union,intersection,none,nonewithoutwarning}]')
> print(' [-verbose_vrt] [-co NAME=VALUE]* [-q]')
> print('')
> print('Create a dataset resulting from a pansharpening operation.')
> return -1
>
>
> def gdal_pansharpen(argv):
>
> argv = gdal.GeneralCmdLineProcessor( argv )
> if argv is None:
> return -1
>
> pan_name = None
> last_name = None
> spectral_ds = []
> spectral_bands = []
> out_name = None
> bands = []
> weights = []
> format = 'VRT'
> creation_options = []
> callback = gdal.TermProgress
> resampling = None
> spat_adjust = None
> verbose_vrt = False
> num_threads = None
> bitdepth = None
> nodata = None
>
> i = 1
> argc = len(argv)
> while i < argc:
> if argv[i] == '-of' and i < len(argv)-1:
> format = argv[i+1]
> print("format " + format)
> i = i + 1
> elif argv[i] == '-r' and i < len(argv)-1:
> resampling = argv[i+1]
> print("resampling " + resampling)
> i = i + 1
> elif argv[i] == '-spat_adjust' and i < len(argv)-1:
> spat_adjust = argv[i+1]
> i = i + 1
> elif argv[i] == '-b' and i < len(argv)-1:
> bands.append(int(argv[i+1]))
> i = i + 1
> elif argv[i] == '-w' and i < len(argv)-1:
> weights.append(float(argv[i+1]))
> i = i + 1
> elif argv[i] == '-co' and i < len(argv)-1:
> creation_options.append(argv[i+1])
> i = i + 1
> elif argv[i] == '-threads' and i < len(argv)-1:
> num_threads = argv[i+1]
> i = i + 1
> elif argv[i] == '-bitdepth' and i < len(argv)-1:
> bitdepth = argv[i+1]
> i = i + 1
> elif argv[i] == '-nodata' and i < len(argv)-1:
> nodata = argv[i+1]
> i = i + 1
> elif argv[i] == '-q':
> callback = None
> elif argv[i] == '-verbose_vrt':
> verbose_vrt = True
> elif argv[i][0] == '-':
> sys.stderr.write('Unrecognized option : %s\n' % argv[i])
> return Usage()
> elif pan_name is None:
> pan_name = argv[i]
> pan_ds = gdal.Open(pan_name)
> if pan_ds is None:
> return 1
> else:
> if last_name is not None:
> pos = last_name.find(',band=')
> if pos > 0:
> spectral_name = last_name[0:pos]
> ds = gdal.Open(spectral_name)
> if ds is None:
> return 1
> band_num = int(last_name[pos+len(',band='):])
> band = ds.GetRasterBand(band_num)
> spectral_ds.append(ds)
> spectral_bands.append(band)
> else:
> spectral_name = last_name
> ds = gdal.Open(spectral_name)
> if ds is None:
> return 1
> for j in range(ds.RasterCount):
> spectral_ds.append(ds)
> spectral_bands.append(ds.GetRasterBand(j+1))
>
> last_name = argv[i]
>
> i = i + 1
> if pan_name is None or len(spectral_bands) == 0:
> return Usage()
> out_name = last_name
> if len(bands) == 0:
> bands = [ j+1 for j in range(len(spectral_bands)) ]
> else:
> for i in range(len(bands)):
> if bands[i] < 0 or bands[i] > len(spectral_bands):
> print('Invalid band number in -b: %d' % bands[i])
> return 1
>
> if len(weights) != 0 and len(weights) != len(spectral_bands):
> print('There must be as many -w values specified as input
> spectral bands')
> return 1
>
>
>
>
> vrt_xml = """<VRTDataset subClass="VRTPansharpenedDataset">\n"""
> if bands != [ j+1 for j in range(len(spectral_bands)) ]:
> for i in range(len(bands)):
> band = spectral_bands[bands[i]-1]
> datatype = gdal.GetDataTypeName(band.DataType)
> colorname =
> gdal.GetColorInterpretationName(band.GetColorInterpretation())
> vrt_xml += """ <VRTRasterBand dataType="%s" band="%d"
> subClass="VRTPansharpenedRasterBand">
> <ColorInterp>%s</ColorInterp>
> </VRTRasterBand>\n""" % (datatype, i+1, colorname)
>
> vrt_xml += """ <PansharpeningOptions>\n"""
> if len(weights) != 0:
> vrt_xml += """ <AlgorithmOptions>\n"""
> vrt_xml += """ <Weights>"""
> for i in range(len(weights)):
> if i > 0: vrt_xml += ","
> vrt_xml += "%.16g" % weights[i]
> vrt_xml += "</Weights>\n"
> vrt_xml += """ </AlgorithmOptions>\n"""
>
> if resampling is not None:
> vrt_xml += ' <Resampling>%s</Resampling>\n' % resampling
>
> if num_threads is not None:
> vrt_xml += ' <NumThreads>%s</NumThreads>\n' % num_threads
>
> if bitdepth is not None:
> vrt_xml += ' <BitDepth>%s</BitDepth>\n' % bitdepth
>
> if nodata is not None:
> vrt_xml += ' <NoData>%s</NoData>\n' % nodata
>
> if spat_adjust is not None:
> vrt_xml += '
> <SpatialExtentAdjustment>%s</SpatialExtentAdjustment>\n' % spat_adjust
>
> pan_relative='0'
> print("panname: "+pan_name)
> print(str(os.path.isabs(pan_name)))
> print(str(os.path.relpath(pan_name, os.path.dirname(out_name))))
> if format.upper() == 'VRT':
> if not os.path.isabs(pan_name):
> pan_relative='1'
> pan_name = os.path.relpath(pan_name, os.path.dirname(out_name))
>
> vrt_xml += """ <PanchroBand>
> <SourceFilename relativeToVRT="%s">%s</SourceFilename>
> <SourceBand>1</SourceBand>
> </PanchroBand>\n""" % (pan_relative, pan_name)
> print(vrt_xml)
> for i in range(len(spectral_bands)):
> dstband = ''
> for j in range(len(bands)):
> if i + 1 == bands[j]:
> dstband = ' dstBand="%d"' % (j+1)
> break
>
> ms_relative='0'
> ms_name = spectral_ds[i].GetDescription()
> if format.upper() == 'VRT':
> if not os.path.isabs(ms_name):
> ms_relative='1'
> ms_name = os.path.relpath(ms_name,
> os.path.dirname(out_name))
>
> vrt_xml += """ <SpectralBand%s>
> <SourceFilename relativeToVRT="%s">%s</SourceFilename>
> <SourceBand>%d</SourceBand>
> </SpectralBand>\n""" % (dstband, ms_relative, ms_name,
> spectral_bands[i].GetBand())
>
> vrt_xml += """ </PansharpeningOptions>\n"""
> vrt_xml += """</VRTDataset>\n"""
> if format.upper() == 'VRT':
> f = gdal.VSIFOpenL(out_name, 'wb')
> if f is None:
> print('Cannot create %s' % out_name)
> return 1
> gdal.VSIFWriteL(vrt_xml, 1, len(vrt_xml), f)
> gdal.VSIFCloseL(f)
> if verbose_vrt:
> vrt_ds = gdal.Open(out_name, gdal.GA_Update)
> vrt_ds.SetMetadata(vrt_ds.GetMetadata())
> else:
> vrt_ds = gdal.Open(out_name)
> if vrt_ds is None:
> return 1
> return 0
> vrt_ds = gdal.Open(vrt_xml)
> out_ds = gdal.GetDriverByName(format).CreateCopy(out_name, vrt_ds,
> 0, creation_options, callback = callback)
>
> if out_ds is None:
> return 1
> return 0
> def main():
> return gdal_pansharpen(sys.argv)
>
> if __name__ == '__main__':
> sys.exit(gdal_pansharpen(sys.argv))
>
>
> My input file is a 7 bands geotiff file and the band 8 of an Landsat8
> images.
>
> If i use the script as follow:
>
> gdal_pansharpen panchro.tif multispectral.tif pansharpend_output.tif
>
> it returns with this error:
>
> ERROR 1: Missing one of rasterXSize, rasterYSize or bands on VRTDataset.
I guess it is because your GDAL python bindings are linking against GDAL <
2.1.0 .
--
Spatialys - Geospatial professional services
http://www.spatialys.com
More information about the gdal-dev
mailing list