1#!/usr/bin/env python3
2###############################################################################
3# $Id: gdal_merge.py 05a98f65afa7e5048cae16a58c2268510cd89103 2021-01-15 16:16:00 +0200 Idan Miara $
4#
5# Project:  InSAR Peppers
6# Purpose:  Module to extract data from many rasters into one output.
7# Author:   Frank Warmerdam, warmerdam@pobox.com
8#
9###############################################################################
10# Copyright (c) 2000, Atlantis Scientific Inc. (www.atlsci.com)
11# Copyright (c) 2009-2011, Even Rouault <even dot rouault at spatialys.com>
12#
13# This library is free software; you can redistribute it and/or
14# modify it under the terms of the GNU Library General Public
15# License as published by the Free Software Foundation; either
16# version 2 of the License, or (at your option) any later version.
17#
18# This library is distributed in the hope that it will be useful,
19# but WITHOUT ANY WARRANTY; without even the implied warranty of
20# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
21# Library General Public License for more details.
22#
23# You should have received a copy of the GNU Library General Public
24# License along with this library; if not, write to the
25# Free Software Foundation, Inc., 59 Temple Place - Suite 330,
26# Boston, MA 02111-1307, USA.
27###############################################################################
28# changes 29Apr2011
29# If the input image is a multi-band one, use all the channels in
30# building the stack.
31# anssi.pekkarinen@fao.org
32
33import math
34import sys
35import time
36
37from osgeo import gdal
38from osgeo_utils.auxiliary.util import GetOutputDriverFor
39
40progress = gdal.TermProgress_nocb
41
42__version__ = '$id$'[5:-1]
43
44
45# =============================================================================
46def raster_copy(s_fh, s_xoff, s_yoff, s_xsize, s_ysize, s_band_n,
47                t_fh, t_xoff, t_yoff, t_xsize, t_ysize, t_band_n,
48                nodata=None, verbose=0):
49
50    if verbose != 0:
51        print('Copy %d,%d,%d,%d to %d,%d,%d,%d.'
52              % (s_xoff, s_yoff, s_xsize, s_ysize,
53                 t_xoff, t_yoff, t_xsize, t_ysize))
54
55    if nodata is not None:
56        return raster_copy_with_nodata(
57            s_fh, s_xoff, s_yoff, s_xsize, s_ysize, s_band_n,
58            t_fh, t_xoff, t_yoff, t_xsize, t_ysize, t_band_n,
59            nodata)
60
61    s_band = s_fh.GetRasterBand(s_band_n)
62    m_band = None
63    # Works only in binary mode and doesn't take into account
64    # intermediate transparency values for compositing.
65    if s_band.GetMaskFlags() != gdal.GMF_ALL_VALID:
66        m_band = s_band.GetMaskBand()
67    elif s_band.GetColorInterpretation() == gdal.GCI_AlphaBand:
68        m_band = s_band
69    if m_band is not None:
70        return raster_copy_with_mask(
71            s_fh, s_xoff, s_yoff, s_xsize, s_ysize, s_band_n,
72            t_fh, t_xoff, t_yoff, t_xsize, t_ysize, t_band_n,
73            m_band)
74
75    s_band = s_fh.GetRasterBand(s_band_n)
76    t_band = t_fh.GetRasterBand(t_band_n)
77
78    data = s_band.ReadRaster(s_xoff, s_yoff, s_xsize, s_ysize,
79                             t_xsize, t_ysize, t_band.DataType)
80    t_band.WriteRaster(t_xoff, t_yoff, t_xsize, t_ysize,
81                       data, t_xsize, t_ysize, t_band.DataType)
82
83    return 0
84
85# =============================================================================
86
87
88def raster_copy_with_nodata(s_fh, s_xoff, s_yoff, s_xsize, s_ysize, s_band_n,
89                            t_fh, t_xoff, t_yoff, t_xsize, t_ysize, t_band_n,
90                            nodata):
91    import numpy as np
92
93    s_band = s_fh.GetRasterBand(s_band_n)
94    t_band = t_fh.GetRasterBand(t_band_n)
95
96    data_src = s_band.ReadAsArray(s_xoff, s_yoff, s_xsize, s_ysize,
97                                  t_xsize, t_ysize)
98    data_dst = t_band.ReadAsArray(t_xoff, t_yoff, t_xsize, t_ysize)
99
100    if not np.isnan(nodata):
101        nodata_test = np.equal(data_src, nodata)
102    else:
103        nodata_test = np.isnan(data_src)
104
105    to_write = np.choose(nodata_test, (data_src, data_dst))
106
107    t_band.WriteArray(to_write, t_xoff, t_yoff)
108
109    return 0
110
111# =============================================================================
112
113
114def raster_copy_with_mask(s_fh, s_xoff, s_yoff, s_xsize, s_ysize, s_band_n,
115                          t_fh, t_xoff, t_yoff, t_xsize, t_ysize, t_band_n,
116                          m_band):
117    import numpy as np
118
119    s_band = s_fh.GetRasterBand(s_band_n)
120    t_band = t_fh.GetRasterBand(t_band_n)
121
122    data_src = s_band.ReadAsArray(s_xoff, s_yoff, s_xsize, s_ysize,
123                                  t_xsize, t_ysize)
124    data_mask = m_band.ReadAsArray(s_xoff, s_yoff, s_xsize, s_ysize,
125                                   t_xsize, t_ysize)
126    data_dst = t_band.ReadAsArray(t_xoff, t_yoff, t_xsize, t_ysize)
127
128    mask_test = np.equal(data_mask, 0)
129    to_write = np.choose(mask_test, (data_src, data_dst))
130
131    t_band.WriteArray(to_write, t_xoff, t_yoff)
132
133    return 0
134
135# =============================================================================
136
137
138def names_to_fileinfos(names):
139    """
140    Translate a list of GDAL filenames, into file_info objects.
141
142    names -- list of valid GDAL dataset names.
143
144    Returns a list of file_info objects.  There may be less file_info objects
145    than names if some of the names could not be opened as GDAL files.
146    """
147
148    file_infos = []
149    for name in names:
150        fi = file_info()
151        if fi.init_from_name(name) == 1:
152            file_infos.append(fi)
153
154    return file_infos
155
156# *****************************************************************************
157
158
159class file_info(object):
160    """A class holding information about a GDAL file."""
161
162    def __init__(self):
163        self.band_type = None
164        self.bands = None
165        self.ct = None
166        self.filename = None
167        self.geotransform = None
168        self.lrx = None
169        self.lry = None
170        self.projection = None
171        self.ulx = None
172        self.uly = None
173        self.xsize = None
174        self.ysize = None
175
176    def init_from_name(self, filename):
177        """
178        Initialize file_info from filename
179
180        filename -- Name of file to read.
181
182        Returns 1 on success or 0 if the file can't be opened.
183        """
184        fh = gdal.Open(filename)
185        if fh is None:
186            return 0
187
188        self.filename = filename
189        self.bands = fh.RasterCount
190        self.xsize = fh.RasterXSize
191        self.ysize = fh.RasterYSize
192        self.band_type = fh.GetRasterBand(1).DataType
193        self.projection = fh.GetProjection()
194        self.geotransform = fh.GetGeoTransform()
195        self.ulx = self.geotransform[0]
196        self.uly = self.geotransform[3]
197        self.lrx = self.ulx + self.geotransform[1] * self.xsize
198        self.lry = self.uly + self.geotransform[5] * self.ysize
199
200        ct = fh.GetRasterBand(1).GetRasterColorTable()
201        if ct is not None:
202            self.ct = ct.Clone()
203        else:
204            self.ct = None
205
206        return 1
207
208    def report(self):
209        print('Filename: ' + self.filename)
210        print('File Size: %dx%dx%d'
211              % (self.xsize, self.ysize, self.bands))
212        print('Pixel Size: %f x %f'
213              % (self.geotransform[1], self.geotransform[5]))
214        print('UL:(%f,%f)   LR:(%f,%f)'
215              % (self.ulx, self.uly, self.lrx, self.lry))
216
217    def copy_into(self, t_fh, s_band=1, t_band=1, nodata_arg=None, verbose=0):
218        """
219        Copy this files image into target file.
220
221        This method will compute the overlap area of the file_info objects
222        file, and the target gdal.Dataset object, and copy the image data
223        for the common window area.  It is assumed that the files are in
224        a compatible projection ... no checking or warping is done.  However,
225        if the destination file is a different resolution, or different
226        image pixel type, the appropriate resampling and conversions will
227        be done (using normal GDAL promotion/demotion rules).
228
229        t_fh -- gdal.Dataset object for the file into which some or all
230        of this file may be copied.
231
232        Returns 1 on success (or if nothing needs to be copied), and zero one
233        failure.
234        """
235        t_geotransform = t_fh.GetGeoTransform()
236        t_ulx = t_geotransform[0]
237        t_uly = t_geotransform[3]
238        t_lrx = t_geotransform[0] + t_fh.RasterXSize * t_geotransform[1]
239        t_lry = t_geotransform[3] + t_fh.RasterYSize * t_geotransform[5]
240
241        # figure out intersection region
242        tgw_ulx = max(t_ulx, self.ulx)
243        tgw_lrx = min(t_lrx, self.lrx)
244        if t_geotransform[5] < 0:
245            tgw_uly = min(t_uly, self.uly)
246            tgw_lry = max(t_lry, self.lry)
247        else:
248            tgw_uly = max(t_uly, self.uly)
249            tgw_lry = min(t_lry, self.lry)
250
251        # do they even intersect?
252        if tgw_ulx >= tgw_lrx:
253            return 1
254        if t_geotransform[5] < 0 and tgw_uly <= tgw_lry:
255            return 1
256        if t_geotransform[5] > 0 and tgw_uly >= tgw_lry:
257            return 1
258
259        # compute target window in pixel coordinates.
260        tw_xoff = int((tgw_ulx - t_geotransform[0]) / t_geotransform[1] + 0.1)
261        tw_yoff = int((tgw_uly - t_geotransform[3]) / t_geotransform[5] + 0.1)
262        tw_xsize = int((tgw_lrx - t_geotransform[0]) / t_geotransform[1] + 0.5) \
263            - tw_xoff
264        tw_ysize = int((tgw_lry - t_geotransform[3]) / t_geotransform[5] + 0.5) \
265            - tw_yoff
266
267        if tw_xsize < 1 or tw_ysize < 1:
268            return 1
269
270        # Compute source window in pixel coordinates.
271        sw_xoff = int((tgw_ulx - self.geotransform[0]) / self.geotransform[1] + 0.1)
272        sw_yoff = int((tgw_uly - self.geotransform[3]) / self.geotransform[5] + 0.1)
273        sw_xsize = int((tgw_lrx - self.geotransform[0]) /
274                       self.geotransform[1] + 0.5) - sw_xoff
275        sw_ysize = int((tgw_lry - self.geotransform[3]) /
276                       self.geotransform[5] + 0.5) - sw_yoff
277
278        if sw_xsize < 1 or sw_ysize < 1:
279            return 1
280
281        # Open the source file, and copy the selected region.
282        s_fh = gdal.Open(self.filename)
283
284        return raster_copy(s_fh, sw_xoff, sw_yoff, sw_xsize, sw_ysize, s_band,
285                           t_fh, tw_xoff, tw_yoff, tw_xsize, tw_ysize, t_band,
286                           nodata_arg, verbose)
287
288
289# =============================================================================
290def Usage():
291    print('Usage: gdal_merge.py [-o out_filename] [-of out_format] [-co NAME=VALUE]*')
292    print('                     [-ps pixelsize_x pixelsize_y] [-tap] [-separate] [-q] [-v] [-pct]')
293    print('                     [-ul_lr ulx uly lrx lry] [-init "value [value...]"]')
294    print('                     [-n nodata_value] [-a_nodata output_nodata_value]')
295    print('                     [-ot datatype] [-createonly] input_files')
296    print('                     [--help-general]')
297    print('')
298    return 1
299
300
301def main(argv=None):
302    verbose = 0
303    quiet = 0
304    names = []
305    frmt = None
306    out_file = 'out.tif'
307
308    ulx = None
309    psize_x = None
310    separate = 0
311    copy_pct = 0
312    nodata = None
313    a_nodata = None
314    create_options = []
315    pre_init = []
316    band_type = None
317    createonly = 0
318    bTargetAlignedPixels = False
319    start_time = time.time()
320
321    if argv is None:
322        argv = argv
323    argv = gdal.GeneralCmdLineProcessor(argv)
324    if argv is None:
325        return 0
326
327    # Parse command line arguments.
328    i = 1
329    while i < len(argv):
330        arg = argv[i]
331
332        if arg == '-o':
333            i = i + 1
334            out_file = argv[i]
335
336        elif arg == '-v':
337            verbose = 1
338
339        elif arg == '-q' or arg == '-quiet':
340            quiet = 1
341
342        elif arg == '-createonly':
343            createonly = 1
344
345        elif arg == '-separate':
346            separate = 1
347
348        elif arg == '-seperate':
349            separate = 1
350
351        elif arg == '-pct':
352            copy_pct = 1
353
354        elif arg == '-ot':
355            i = i + 1
356            band_type = gdal.GetDataTypeByName(argv[i])
357            if band_type == gdal.GDT_Unknown:
358                print('Unknown GDAL data type: %s' % argv[i])
359                return 1
360
361        elif arg == '-init':
362            i = i + 1
363            str_pre_init = argv[i].split()
364            for x in str_pre_init:
365                pre_init.append(float(x))
366
367        elif arg == '-n':
368            i = i + 1
369            nodata = float(argv[i])
370
371        elif arg == '-a_nodata':
372            i = i + 1
373            a_nodata = float(argv[i])
374
375        elif arg == '-f' or arg == '-of':
376            i = i + 1
377            frmt = argv[i]
378
379        elif arg == '-co':
380            i = i + 1
381            create_options.append(argv[i])
382
383        elif arg == '-ps':
384            psize_x = float(argv[i + 1])
385            psize_y = -1 * abs(float(argv[i + 2]))
386            i = i + 2
387
388        elif arg == '-tap':
389            bTargetAlignedPixels = True
390
391        elif arg == '-ul_lr':
392            ulx = float(argv[i + 1])
393            uly = float(argv[i + 2])
394            lrx = float(argv[i + 3])
395            lry = float(argv[i + 4])
396            i = i + 4
397
398        elif arg[:1] == '-':
399            print('Unrecognized command option: %s' % arg)
400            return Usage()
401
402        else:
403            names.append(arg)
404
405        i = i + 1
406
407    if not names:
408        print('No input files selected.')
409        return Usage()
410
411    if frmt is None:
412        frmt = GetOutputDriverFor(out_file)
413
414    Driver = gdal.GetDriverByName(frmt)
415    if Driver is None:
416        print('Format driver %s not found, pick a supported driver.' % frmt)
417        return 1
418
419    DriverMD = Driver.GetMetadata()
420    if 'DCAP_CREATE' not in DriverMD:
421        print('Format driver %s does not support creation and piecewise writing.\nPlease select a format that does, such as GTiff (the default) or HFA (Erdas Imagine).' % frmt)
422        return 1
423
424    # Collect information on all the source files.
425    file_infos = names_to_fileinfos(names)
426
427    if ulx is None:
428        ulx = file_infos[0].ulx
429        uly = file_infos[0].uly
430        lrx = file_infos[0].lrx
431        lry = file_infos[0].lry
432
433        for fi in file_infos:
434            ulx = min(ulx, fi.ulx)
435            uly = max(uly, fi.uly)
436            lrx = max(lrx, fi.lrx)
437            lry = min(lry, fi.lry)
438
439    if psize_x is None:
440        psize_x = file_infos[0].geotransform[1]
441        psize_y = file_infos[0].geotransform[5]
442
443    if band_type is None:
444        band_type = file_infos[0].band_type
445
446    # Try opening as an existing file.
447    gdal.PushErrorHandler('CPLQuietErrorHandler')
448    t_fh = gdal.Open(out_file, gdal.GA_Update)
449    gdal.PopErrorHandler()
450
451    # Create output file if it does not already exist.
452    if t_fh is None:
453
454        if bTargetAlignedPixels:
455            ulx = math.floor(ulx / psize_x) * psize_x
456            lrx = math.ceil(lrx / psize_x) * psize_x
457            lry = math.floor(lry / -psize_y) * -psize_y
458            uly = math.ceil(uly / -psize_y) * -psize_y
459
460        geotransform = [ulx, psize_x, 0, uly, 0, psize_y]
461
462        xsize = int((lrx - ulx) / geotransform[1] + 0.5)
463        ysize = int((lry - uly) / geotransform[5] + 0.5)
464
465        if separate != 0:
466            bands = 0
467
468            for fi in file_infos:
469                bands = bands + fi.bands
470        else:
471            bands = file_infos[0].bands
472
473        t_fh = Driver.Create(out_file, xsize, ysize, bands,
474                             band_type, create_options)
475        if t_fh is None:
476            print('Creation failed, terminating gdal_merge.')
477            return 1
478
479        t_fh.SetGeoTransform(geotransform)
480        t_fh.SetProjection(file_infos[0].projection)
481
482        if copy_pct:
483            t_fh.GetRasterBand(1).SetRasterColorTable(file_infos[0].ct)
484    else:
485        if separate != 0:
486            bands = 0
487            for fi in file_infos:
488                bands = bands + fi.bands
489            if t_fh.RasterCount < bands:
490                print('Existing output file has less bands than the input files. You should delete it before. Terminating gdal_merge.')
491                return 1
492        else:
493            bands = min(file_infos[0].bands, t_fh.RasterCount)
494
495    # Do we need to set nodata value ?
496    if a_nodata is not None:
497        for i in range(t_fh.RasterCount):
498            t_fh.GetRasterBand(i + 1).SetNoDataValue(a_nodata)
499
500    # Do we need to pre-initialize the whole mosaic file to some value?
501    if pre_init is not None:
502        if t_fh.RasterCount <= len(pre_init):
503            for i in range(t_fh.RasterCount):
504                t_fh.GetRasterBand(i + 1).Fill(pre_init[i])
505        elif len(pre_init) == 1:
506            for i in range(t_fh.RasterCount):
507                t_fh.GetRasterBand(i + 1).Fill(pre_init[0])
508
509    # Copy data from source files into output file.
510    t_band = 1
511
512    if quiet == 0 and verbose == 0:
513        progress(0.0)
514    fi_processed = 0
515
516    for fi in file_infos:
517        if createonly != 0:
518            continue
519
520        if verbose != 0:
521            print("")
522            print("Processing file %5d of %5d, %6.3f%% completed in %d minutes."
523                  % (fi_processed + 1, len(file_infos),
524                     fi_processed * 100.0 / len(file_infos),
525                     int(round((time.time() - start_time) / 60.0))))
526            fi.report()
527
528        if separate == 0:
529            for band in range(1, bands + 1):
530                fi.copy_into(t_fh, band, band, nodata, verbose)
531        else:
532            for band in range(1, fi.bands + 1):
533                fi.copy_into(t_fh, band, t_band, nodata, verbose)
534                t_band = t_band + 1
535
536        fi_processed = fi_processed + 1
537        if quiet == 0 and verbose == 0:
538            progress(fi_processed / float(len(file_infos)))
539
540    # Force file to be closed.
541    t_fh = None
542
543
544if __name__ == '__main__':
545    sys.exit(main(sys.argv))
546