#!/usr/bin/python2
# -*- coding: utf-8 -*-
# Copyright (c) 2000, Atlantis Scientific Inc. (www.atlsci.com)
# Copyright (C) 2005  Gabriel Ebner <ge@gabrielebner.at>
# Copyright (c) 2009, Even Rouault <even dot rouault at mines-paris dot org>
#
# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Library General Public
# License as published by the Free Software Foundation; either
# version 2 of the License, or (at your option) any later version.
#
# This library is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
# Library General Public License for more details.
#
# You should have received a copy of the GNU Library General Public
# License along with this library; if not, write to the
# Free Software Foundation, Inc., 59 Temple Place - Suite 330,
# Boston, MA 02111-1307, USA.
#
# modified output of geotransform to retain same numerical
# precision as the gdal 'C' utilities
# Norman Vine  nhv@cape.com  03-Oct-2005 6:23:45 am

import sys

from osgeo import gdal

# =============================================================================


def names_to_fileinfos(names):
    """
    Translate a list of GDAL filenames, into file_info objects.

    names -- list of valid GDAL dataset names.

    Returns a list of file_info objects.  There may be less file_info objects
    than names if some of the names could not be opened as GDAL files.
    """

    file_infos = []
    for name in names:
        fi = file_info()
        if fi.init_from_name(name):
            file_infos.append(fi)
        else:
            print('Can not open dataset "%s", skipped' % name)

    return file_infos

# *****************************************************************************


class file_info(object):
    """A class holding information about a GDAL file."""

    def __init__(self):
        self.band_types = [None]
        self.block_sizes = [None]
        self.nodata = [None]
        self.cts = [None]
        self.color_interps = [None]

        self.filename = None
        self.bands = None
        self.xsize = None
        self.ysize = None
        self.projection = None
        self.geotransform = None
        self.ulx = None
        self.uly = None
        self.lrx = None
        self.lry = None

    def init_from_name(self, filename):
        """
        Initialize file_info from filename

        filename -- Name of file to read.

        Returns 1 on success or 0 if the file can't be opened.
        """
        fh = gdal.Open(filename)
        if fh is None:
            return False

        self.filename = filename
        self.bands = fh.RasterCount
        self.xsize = fh.RasterXSize
        self.ysize = fh.RasterYSize
        self.projection = fh.GetProjection()
        self.geotransform = fh.GetGeoTransform()
        self.ulx = self.geotransform[0]
        self.uly = self.geotransform[3]
        self.lrx = self.ulx + self.geotransform[1] * self.xsize
        self.lry = self.uly + self.geotransform[5] * self.ysize

        self.band_types = [None]
        self.block_sizes = [None]
        self.nodata = [None]
        self.cts = [None]
        self.color_interps = [None]
        for i in range(1, fh.RasterCount + 1):
            band = fh.GetRasterBand(i)
            self.band_types.append(band.DataType)
            self.block_sizes.append(band.GetBlockSize())
            if band.GetNoDataValue() is not None:
                self.nodata.append(band.GetNoDataValue())
            self.color_interps.append(band.GetRasterColorInterpretation())
            ct = band.GetRasterColorTable()
            if ct is not None:
                self.cts.append(ct.Clone())
            else:
                self.cts.append(None)

        return True

    def write_source(self, t_fh, t_geotransform, xsize, ysize, s_band):
        t_ulx = t_geotransform[0]
        t_uly = t_geotransform[3]
        t_lrx = t_geotransform[0] + xsize * t_geotransform[1]
        t_lry = t_geotransform[3] + ysize * t_geotransform[5]

        # figure out intersection region
        tgw_ulx = max(t_ulx, self.ulx)
        tgw_lrx = min(t_lrx, self.lrx)
        if t_geotransform[5] < 0:
            tgw_uly = min(t_uly, self.uly)
            tgw_lry = max(t_lry, self.lry)
        else:
            tgw_uly = max(t_uly, self.uly)
            tgw_lry = min(t_lry, self.lry)

        # do they even intersect?
        if tgw_ulx >= tgw_lrx:
            return 1
        if t_geotransform[5] < 0 and tgw_uly <= tgw_lry:
            return 1
        if t_geotransform[5] > 0 and tgw_uly >= tgw_lry:
            return 1

        # compute target window in pixel coordinates.
        tw_xoff = int((tgw_ulx - t_geotransform[0]) / t_geotransform[1] + 0.1)
        tw_yoff = int((tgw_uly - t_geotransform[3]) / t_geotransform[5] + 0.1)
        tw_xsize = int((tgw_lrx - t_geotransform[0]) / t_geotransform[1] + 0.5) \
            - tw_xoff
        tw_ysize = int((tgw_lry - t_geotransform[3]) / t_geotransform[5] + 0.5) \
            - tw_yoff

        if tw_xsize < 1 or tw_ysize < 1:
            return 1

        # Compute source window in pixel coordinates.
        sw_xoff = int((tgw_ulx - self.geotransform[0]) / self.geotransform[1])
        sw_yoff = int((tgw_uly - self.geotransform[3]) / self.geotransform[5])
        sw_xsize = int((tgw_lrx - self.geotransform[0]) /
                       self.geotransform[1] + 0.5) - sw_xoff
        sw_ysize = int((tgw_lry - self.geotransform[3]) /
                       self.geotransform[5] + 0.5) - sw_yoff

        if sw_xsize < 1 or sw_ysize < 1:
            return 1

        t_fh.write('\t\t<SimpleSource>\n')
        t_fh.write(('\t\t\t<SourceFilename relativeToVRT="1">%s' +
                    '</SourceFilename>\n') % self.filename)
        t_fh.write('\t\t\t<SourceBand>%i</SourceBand>\n' % s_band)
        data_type_name = gdal.GetDataTypeName(self.band_types[s_band])
        block_size_x, block_size_y = self.block_sizes[s_band]
        t_fh.write('\t\t\t<SourceProperties RasterXSize="%i" RasterYSize="%i"'
                   ' DataType="%s" BlockXSize="%i" BlockYSize="%i"/>\n'
                   % (self.xsize, self.ysize, data_type_name, block_size_x, block_size_y))
        t_fh.write('\t\t\t<SrcRect xOff="%i" yOff="%i" xSize="%i" ySize="%i"/>\n'
                   % (sw_xoff, sw_yoff, sw_xsize, sw_ysize))
        t_fh.write('\t\t\t<DstRect xOff="%i" yOff="%i" xSize="%i" ySize="%i"/>\n'
                   % (tw_xoff, tw_yoff, tw_xsize, tw_ysize))
        t_fh.write('\t\t</SimpleSource>\n')

# =============================================================================


def Usage():
    print('Usage: gdal_vrtmerge.py [-o out_filename] [-separate] [-pct]')
    print('           [-ul_lr ulx uly lrx lry] [-ot datatype] [-i input_file_list')
    print('           | input_files]')

# =============================================================================
#
# Program mainline.
#


if __name__ == '__main__':

    names = []
    out_file = 'out.vrt'

    ulx = None
    psize_x = None
    separate = False
    pre_init = None

    argv = gdal.GeneralCmdLineProcessor(sys.argv)
    if argv is None:
        sys.exit(0)

    # Parse command line arguments.
    i = 1
    while i < len(argv):
        arg = argv[i]

        if arg == '-o':
            i = i + 1
            out_file = argv[i]

        elif arg == '-i':
            i = i + 1
            in_file_list = open(argv[i])
            names.extend(in_file_list.read().split())

        elif arg == '-separate':
            separate = True

        elif arg == '-ul_lr':
            ulx = float(argv[i + 1])
            uly = float(argv[i + 2])
            lrx = float(argv[i + 3])
            lry = float(argv[i + 4])
            i = i + 4

        elif arg[:1] == '-':
            print('Unrecognized command option: ', arg)
            Usage()
            sys.exit(1)

        else:
            names.append(arg)

        i = i + 1

    if not names:
        print('No input files selected.')
        Usage()
        sys.exit(1)

    # Collect information on all the source files.
    file_infos = names_to_fileinfos(names)
    if not file_infos:
        print('Nothing to process, exiting.')
        sys.exit(1)

    if ulx is None:
        ulx = file_infos[0].ulx
        uly = file_infos[0].uly
        lrx = file_infos[0].lrx
        lry = file_infos[0].lry

        for fi in file_infos:
            ulx = min(ulx, fi.ulx)
            uly = max(uly, fi.uly)
            lrx = max(lrx, fi.lrx)
            lry = min(lry, fi.lry)

    if psize_x is None:
        psize_x = file_infos[0].geotransform[1]
        psize_y = file_infos[0].geotransform[5]

    projection = file_infos[0].projection

    for fi in file_infos:
        if fi.geotransform[1] != psize_x or fi.geotransform[5] != psize_y:
            print("All files must have the same scale; %s does not"
                  % fi.filename)
            sys.exit(1)

        if fi.geotransform[2] != 0 or fi.geotransform[4] != 0:
            print("No file must be rotated; %s is" % fi.filename)
            sys.exit(1)

        if fi.projection != projection:
            print("All files must be in the same projection; %s is not"
                  % fi.filename)
            sys.exit(1)

    geotransform = (ulx, psize_x, 0.0, uly, 0.0, psize_y)

    xsize = int(((lrx - ulx) / geotransform[1]) + 0.5)
    ysize = int(((lry - uly) / geotransform[5]) + 0.5)

    if separate:
        bands = len(file_infos)
    else:
        bands = file_infos[0].bands

    t_fh = open(out_file, 'w')
    t_fh.write('<VRTDataset rasterXSize="%i" rasterYSize="%i">\n'
               % (xsize, ysize))
    t_fh.write('\t<GeoTransform>%24.16f, %24.16f, %24.16f, %24.16f, %24.16f, %24.16f</GeoTransform>\n'
               % geotransform)

    if projection:
        t_fh.write('\t<SRS>%s</SRS>\n' % projection)

    if separate:
        band_n = 0
        for fi in file_infos:
            band_n = band_n + 1
            if len(fi.band_types) != 2:
                print('File %s has %d bands. Only first band will be taken '
                      'into account' % (fi.filename, len(fi.band_types) - 1))
            dataType = gdal.GetDataTypeName(fi.band_types[1])

            t_fh.write('\t<VRTRasterBand dataType="%s" band="%i">\n'
                       % (dataType, band_n))
            t_fh.write('\t\t<ColorInterp>%s</ColorInterp>\n' %
                       gdal.GetColorInterpretationName(fi.color_interps[1]))
            fi.write_source(t_fh, geotransform, xsize, ysize, 1)
            t_fh.write('\t</VRTRasterBand>\n')
    else:
        for band in range(1, bands + 1):
            dataType = gdal.GetDataTypeName(file_infos[0].band_types[band])

            t_fh.write('\t<VRTRasterBand dataType="%s" band="%i">\n'
                       % (dataType, band))
            if file_infos[0].nodata != [None]:
                t_fh.write('\t\t<NoDataValue>%f</NoDataValue>\n' %
                           file_infos[0].nodata[band])
            t_fh.write('\t\t<ColorInterp>%s</ColorInterp>\n' %
                       gdal.GetColorInterpretationName(
                           file_infos[0].color_interps[band]))

            ct = file_infos[0].cts[band]
            if ct is not None:
                t_fh.write('\t\t<ColorTable>\n')
                for i in range(ct.GetCount()):
                    t_fh.write(
                        '\t\t\t<Entry c1="%i" c2="%i" c3="%i" c4="%i"/>\n'
                        % ct.GetColorEntry(i))
                t_fh.write('\t\t</ColorTable>\n')

            for fi in file_infos:
                fi.write_source(t_fh, geotransform, xsize, ysize, band)

            t_fh.write('\t</VRTRasterBand>\n')

    t_fh.write('</VRTDataset>\n')