pyModis/scripts/modis_mosaic.py

#!/usr/bin/env python
# script to create mosaic from several tiles of MODIS
#
#  (c) Copyright Luca Delucchi 2012
#  Authors: Luca Delucchi
#  Email: luca dot delucchi at iasma dot it
#
##################################################################
#
#  This MODIS Python script is licensed under the terms of GNU GPL 2.
#  This program is free software; you can redistribute it and/or
#  modify it under the terms of the GNU General Public License as
#  published by the Free Software Foundation; either version 2 of
#  the License, or (at your option) any later version.
#  This program is distributed in the hope that it will be useful,
#  but WITHOUT ANY WARRANTY; without even implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
#  See the GNU General Public License for more details.
#
##################################################################
"""Script to mosaic the input tiles. It is able to use MRT or GDAL as backend
"""

import os
import sys
try:
    from pymodis import optparse_gui
    WXPYTHON = True
except:
    WXPYTHON = False
from pymodis import convertmodis
from pymodis import convertmodis_gdal
from pymodis import optparse_required
from optparse import OptionGroup
try:
    import osgeo.gdal as gdal
except ImportError:
    try:
        import gdal
    except ImportError:
        raise Exception('Python GDAL library not found, please install python-gdal')

ERROR = "You have to define the name of a text file containing HDF files" \
        " (One HDF file for line)."


def main():
    """Main function"""
    # usage
    usage = "usage: %prog [options] hdflist_file"
    if 1 == len(sys.argv) and WXPYTHON:
        option_parser_class = optparse_gui.OptionParser
    else:
        option_parser_class = optparse_required.OptionParser
    parser = option_parser_class(usage=usage, description='modis_mosaic')
    groupR = OptionGroup(parser, 'General options')
    groupG = OptionGroup(parser, 'Options for GDAL')
    groupM = OptionGroup(parser, 'Options for MRT')
    # output
    groupR.add_option("-o", "--output", dest="output", required=True,
                      help="the name or prefix (for VRT) of output mosaic",
                      metavar="OUTPUT_FILE")
    # subset
    groupR.add_option("-s", "--subset", dest="subset",
                      help="a subset of product layers. The string should"
                      " be similar to: 1 0 [default: all layers]")
    # options only for GDAL
    groupG.add_option("-f", "--output-format", dest="output_format",
                      metavar="OUTPUT_FORMAT", default="GTiff",
                      help="output format supported: GTiff, HDF4Image"
                      " [default=%default]")
#    groupG.add_option("-g", "--grain", dest="grain",
#                      type="int", help="force the spatial resolution of output"
#                      " file")
    # options for set VRT
    groupG.add_option("-v", "--vrt", dest="vrt", action="store_true",
                      default=False, help="Create a GDAL VRT file. No other "
                      "GDAL options have to been set")
    # mrt path
    groupM.add_option("-m", "--mrt", dest="mrt_path", type='directory',
                      help="the path to MRT software", metavar="MRT_PATH")
    parser.add_option_group(groupR)
    parser.add_option_group(groupG)
    parser.add_option_group(groupM)
    (options, args) = parser.parse_args()
    # check the number of tiles
    if len(args) == 0 and not WXPYTHON:
        parser.print_help()
        sys.exit(1)
    if not args:
        parser.error(ERROR)
    else:
        if not isinstance(args, list):
            parser.error(ERROR)
        elif len(args) > 1:
            parser.error(ERROR)
        if not os.path.isfile(args[0]):
            parser.error(ERROR + '. ' + args[0] + ' does not exists')

    if not os.path.isfile(args[0]):
        parser.error("You have to define the name of a text file containing "
                     "HDF files. (One HDF file for line)")

    # check is a subset it is set
    if not options.subset:
        options.subset = False
    else:
        if  (options.subset.strip().startswith('(') or options.subset.strip().endswith(')')):
            parser.error('ERROR: The spectral string should be similar to: '
                         '"1 0" without "(" and ")"')
#    if not options.grain and options.vrt:
#        options.grain = False
#    elif not options.grain and options.vrt:
#        parser.error("You have to define the resolution of output file. Please"
#                     " set -g/--grain option")
    if options.mrt_path:
        modisOgg = convertmodis.createMosaic(args[0], options.output,
                                             options.mrt_path, options.subset)
        modisOgg.run()
    else:
        tiles = dict()
        dire = os.path.dirname(args[0])
        with open(args[0]) as f:
            for l in f:
                name = os.path.splitext(l.strip())[0]
                day = name.split('.')[1]
                if day not in tiles.keys():
                    tiles[day] = list()
                if '.hdf' not in name:
                    if dire not in l:
                        fname = os.path.join(dire, l.strip())
                    else:
                        fname = l.strip()
                    tiles[day].append(fname)

        for day in tiles.keys():
            modisOgg = convertmodis_gdal.createMosaicGDAL(tiles[day],
                                                          options.subset,
                                                          options.output_format)
            output = "{da}_{fi}".format(da=day,  fi=options.output)
            if options.vrt:
                modisOgg.write_vrt(output)
            else:
                modisOgg.run(output)

if __name__ == "__main__":
    gdal.AllRegister()
    argv = gdal.GeneralCmdLineProcessor(sys.argv)
    if argv is not None:
        main()