[OSGeo-Discuss] (no subject)
Lorenzo Bernardi
lorenzo.bernardi at nais-solutions.it
Fri Apr 22 08:09:05 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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/discuss/attachments/20160422/ebbf4ac9/attachment.html>
More information about the Discuss
mailing list