[gdal-dev] Problem with gdal_pansharpen tool

Lorenzo Bernardi lorenzo.bernardi at nais-solutions.it
Fri Apr 22 08:37:19 PDT 2016


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 undestrand that the problem is the bands on the VRTDataset but for now i
didn't find any solutions.

I thank you in advance for any suggestions.

Thank you,

Lorenzo Bernardi


-- 

Lorenzo Bernardi



-- 

Lorenzo Bernardi
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20160422/813c45b9/attachment-0001.html>


More information about the gdal-dev mailing list