Geografic-Information-System/scripts/i.image.mosaic/i.image.mosaic.py

#!/usr/bin/env python3

# written by Markus Neteler 18. August 1998 / 20. Jan. 1999
#            neteler geog.uni-hannover.de
# mosaic code from Felix Gershunov (Felix spsl.nsc.ru)
# updated for GRASS 5.7 by Michael Barton 2004/04/05
# converted to Python by Glynn Clements
#
# COPYRIGHT:    (C) 1999,2007,2008 by the GRASS Development Team
#
#               This program is free software under the GNU General Public
#               License (>=v2). Read the file COPYING that comes with GRASS
#               for details.
#
# TODO: - implement g.findfile for 3 and 4 maps (currently only current mapset supported)
#            [done for 2 maps]
#       - fix isnull() in r.mapcalc for 3 and 4 maps composites
#            [done for 2 maps]
#       - fix color table length (currently only 256 cols supported, make
#         flexible)
#            [done for 2 maps]
# --------------------------------------------------


# %module
# % description: Mosaics several images and extends colormap.
# % keyword: imagery
# % keyword: geometry
# % keyword: mosaicking
# %end
# %option G_OPT_R_INPUTS
# %end
# %option G_OPT_R_OUTPUT
# %end
from __future__ import print_function

import grass.script as gscript


def copy_colors(fh, map, offset):
    p = gscript.pipe_command("r.colors.out", map=map)
    for line in p.stdout:
        f = gscript.decode(line).rstrip("\r\n").split(" ")
        if offset:
            if f[0] in ["nv", "default"]:
                continue
            f[0] = str(float(f[0]) + offset)
        fh.write(gscript.encode(" ".join(f) + "\n"))
    p.wait()


def get_limit(map):
    return gscript.raster_info(map)["max"]


def make_expression(i, count):
    if i > count:
        return "null()"
    else:
        e = make_expression(i + 1, count)
        return "if(isnull($image%d),%s,$image%d+$offset%d)" % (i, e, i, i)


def main():
    images = options["input"].split(",")
    output = options["output"]

    count = len(images)
    msg = _("Do not forget to set region properly to cover all images.")
    gscript.warning(msg)

    offset = 0
    offsets = []
    parms = {}
    for n, img in enumerate(images):
        offsets.append(offset)
        parms["image%d" % (n + 1)] = img
        parms["offset%d" % (n + 1)] = offset
        offset += get_limit(img) + 1

    gscript.message(_("Mosaicing %d images...") % count)

    gscript.mapcalc("$output = " + make_expression(1, count), output=output, **parms)

    # modify the color table:
    p = gscript.feed_command("r.colors", map=output, rules="-")
    for img, offset in zip(images, offsets):
        print(img, offset)
        copy_colors(p.stdin, img, offset)
    p.stdin.close()
    p.wait()

    gscript.message(_("Done. Raster map <%s> created.") % output)

    # write cmd history:
    gscript.raster_history(output)


if __name__ == "__main__":
    options, flags = gscript.parser()
    main()