source: branches/0.5/win/PostNAS-0.5/bin/gdal_merge.py @ 23

Revision 23, 14.9 KB checked in by astrid.emde, 12 years ago (diff)
Line 
1#!/usr/bin/env python
2###############################################################################
3# $Id: gdal_merge.py 18412 2009-12-30 12:53:37Z chaitanya $
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#
12# This library is free software; you can redistribute it and/or
13# modify it under the terms of the GNU Library General Public
14# License as published by the Free Software Foundation; either
15# version 2 of the License, or (at your option) any later version.
16#
17# This library is distributed in the hope that it will be useful,
18# but WITHOUT ANY WARRANTY; without even the implied warranty of
19# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
20# Library General Public License for more details.
21#
22# You should have received a copy of the GNU Library General Public
23# License along with this library; if not, write to the
24# Free Software Foundation, Inc., 59 Temple Place - Suite 330,
25# Boston, MA 02111-1307, USA.
26###############################################################################
27
28try:
29    from osgeo import gdal
30    gdal.TermProgress = gdal.TermProgress_nocb
31except ImportError:
32    import gdal
33
34import sys
35import glob
36
37__version__ = '$id$'[5:-1]
38verbose = 0
39quiet = 0
40
41
42# =============================================================================
43def raster_copy( s_fh, s_xoff, s_yoff, s_xsize, s_ysize, s_band_n,
44                 t_fh, t_xoff, t_yoff, t_xsize, t_ysize, t_band_n,
45                 nodata=None ):
46
47    if nodata is not None:
48        return raster_copy_with_nodata(
49            s_fh, s_xoff, s_yoff, s_xsize, s_ysize, s_band_n,
50            t_fh, t_xoff, t_yoff, t_xsize, t_ysize, t_band_n,
51            nodata )
52
53    if verbose != 0:
54        print('Copy %d,%d,%d,%d to %d,%d,%d,%d.' \
55              % (s_xoff, s_yoff, s_xsize, s_ysize,
56             t_xoff, t_yoff, t_xsize, t_ysize ))
57
58    s_band = s_fh.GetRasterBand( s_band_n )
59    t_band = t_fh.GetRasterBand( t_band_n )
60
61    data = s_band.ReadRaster( s_xoff, s_yoff, s_xsize, s_ysize,
62                              t_xsize, t_ysize, t_band.DataType )
63    t_band.WriteRaster( t_xoff, t_yoff, t_xsize, t_ysize,
64                        data, t_xsize, t_ysize, t_band.DataType )
65       
66
67    return 0
68   
69# =============================================================================
70def raster_copy_with_nodata( s_fh, s_xoff, s_yoff, s_xsize, s_ysize, s_band_n,
71                             t_fh, t_xoff, t_yoff, t_xsize, t_ysize, t_band_n,
72                             nodata ):
73    try:
74        import numpy as Numeric
75    except ImportError:
76        import Numeric
77   
78    if verbose != 0:
79        print('Copy %d,%d,%d,%d to %d,%d,%d,%d.' \
80              % (s_xoff, s_yoff, s_xsize, s_ysize,
81             t_xoff, t_yoff, t_xsize, t_ysize ))
82
83    s_band = s_fh.GetRasterBand( s_band_n )
84    t_band = t_fh.GetRasterBand( t_band_n )
85
86    data_src = s_band.ReadAsArray( s_xoff, s_yoff, s_xsize, s_ysize,
87                                   t_xsize, t_ysize )
88    data_dst = t_band.ReadAsArray( t_xoff, t_yoff, t_xsize, t_ysize )
89
90    nodata_test = Numeric.equal(data_src,nodata)
91    to_write = Numeric.choose( nodata_test, (data_src, data_dst) )
92                               
93    t_band.WriteArray( to_write, t_xoff, t_yoff )
94
95    return 0
96   
97# =============================================================================
98def names_to_fileinfos( names ):
99    """
100    Translate a list of GDAL filenames, into file_info objects.
101
102    names -- list of valid GDAL dataset names.
103
104    Returns a list of file_info objects.  There may be less file_info objects
105    than names if some of the names could not be opened as GDAL files.
106    """
107   
108    file_infos = []
109    for name in names:
110        fi = file_info()
111        if fi.init_from_name( name ) == 1:
112            file_infos.append( fi )
113
114    return file_infos
115
116# *****************************************************************************
117class file_info:
118    """A class holding information about a GDAL file."""
119
120    def init_from_name(self, filename):
121        """
122        Initialize file_info from filename
123
124        filename -- Name of file to read.
125
126        Returns 1 on success or 0 if the file can't be opened.
127        """
128        fh = gdal.Open( filename )
129        if fh is None:
130            return 0
131
132        self.filename = filename
133        self.bands = fh.RasterCount
134        self.xsize = fh.RasterXSize
135        self.ysize = fh.RasterYSize
136        self.band_type = fh.GetRasterBand(1).DataType
137        self.projection = fh.GetProjection()
138        self.geotransform = fh.GetGeoTransform()
139        self.ulx = self.geotransform[0]
140        self.uly = self.geotransform[3]
141        self.lrx = self.ulx + self.geotransform[1] * self.xsize
142        self.lry = self.uly + self.geotransform[5] * self.ysize
143
144        ct = fh.GetRasterBand(1).GetRasterColorTable()
145        if ct is not None:
146            self.ct = ct.Clone()
147        else:
148            self.ct = None
149
150        return 1
151
152    def report( self ):
153        print('Filename: '+ self.filename)
154        print('File Size: %dx%dx%d' \
155              % (self.xsize, self.ysize, self.bands))
156        print('Pixel Size: %f x %f' \
157              % (self.geotransform[1],self.geotransform[5]))
158        print('UL:(%f,%f)   LR:(%f,%f)' \
159              % (self.ulx,self.uly,self.lrx,self.lry))
160
161    def copy_into( self, t_fh, s_band = 1, t_band = 1, nodata_arg=None ):
162        """
163        Copy this files image into target file.
164
165        This method will compute the overlap area of the file_info objects
166        file, and the target gdal.Dataset object, and copy the image data
167        for the common window area.  It is assumed that the files are in
168        a compatible projection ... no checking or warping is done.  However,
169        if the destination file is a different resolution, or different
170        image pixel type, the appropriate resampling and conversions will
171        be done (using normal GDAL promotion/demotion rules).
172
173        t_fh -- gdal.Dataset object for the file into which some or all
174        of this file may be copied.
175
176        Returns 1 on success (or if nothing needs to be copied), and zero one
177        failure.
178        """
179        t_geotransform = t_fh.GetGeoTransform()
180        t_ulx = t_geotransform[0]
181        t_uly = t_geotransform[3]
182        t_lrx = t_geotransform[0] + t_fh.RasterXSize * t_geotransform[1]
183        t_lry = t_geotransform[3] + t_fh.RasterYSize * t_geotransform[5]
184
185        # figure out intersection region
186        tgw_ulx = max(t_ulx,self.ulx)
187        tgw_lrx = min(t_lrx,self.lrx)
188        if t_geotransform[5] < 0:
189            tgw_uly = min(t_uly,self.uly)
190            tgw_lry = max(t_lry,self.lry)
191        else:
192            tgw_uly = max(t_uly,self.uly)
193            tgw_lry = min(t_lry,self.lry)
194       
195        # do they even intersect?
196        if tgw_ulx >= tgw_lrx:
197            return 1
198        if t_geotransform[5] < 0 and tgw_uly <= tgw_lry:
199            return 1
200        if t_geotransform[5] > 0 and tgw_uly >= tgw_lry:
201            return 1
202           
203        # compute target window in pixel coordinates.
204        tw_xoff = int((tgw_ulx - t_geotransform[0]) / t_geotransform[1] + 0.1)
205        tw_yoff = int((tgw_uly - t_geotransform[3]) / t_geotransform[5] + 0.1)
206        tw_xsize = int((tgw_lrx - t_geotransform[0])/t_geotransform[1] + 0.5) \
207                   - tw_xoff
208        tw_ysize = int((tgw_lry - t_geotransform[3])/t_geotransform[5] + 0.5) \
209                   - tw_yoff
210
211        if tw_xsize < 1 or tw_ysize < 1:
212            return 1
213
214        # Compute source window in pixel coordinates.
215        sw_xoff = int((tgw_ulx - self.geotransform[0]) / self.geotransform[1])
216        sw_yoff = int((tgw_uly - self.geotransform[3]) / self.geotransform[5])
217        sw_xsize = int((tgw_lrx - self.geotransform[0]) \
218                       / self.geotransform[1] + 0.5) - sw_xoff
219        sw_ysize = int((tgw_lry - self.geotransform[3]) \
220                       / self.geotransform[5] + 0.5) - sw_yoff
221
222        if sw_xsize < 1 or sw_ysize < 1:
223            return 1
224
225        # Open the source file, and copy the selected region.
226        s_fh = gdal.Open( self.filename )
227
228        return \
229            raster_copy( s_fh, sw_xoff, sw_yoff, sw_xsize, sw_ysize, s_band,
230                         t_fh, tw_xoff, tw_yoff, tw_xsize, tw_ysize, t_band,
231                         nodata_arg )
232
233
234# =============================================================================
235def Usage():
236    print('Usage: gdal_merge.py [-o out_filename] [-of out_format] [-co NAME=VALUE]*')
237    print('                     [-ps pixelsize_x pixelsize_y] [-separate] [-q] [-v] [-pct]')
238    print('                     [-ul_lr ulx uly lrx lry] [-n nodata_value] [-init value]')
239    print('                     [-ot datatype] [-createonly] input_files')
240    print('                     [--help-general]')
241    print('')
242
243# =============================================================================
244#
245# Program mainline.
246#
247
248def main( argv=None ):
249
250    names = []
251    format = 'GTiff'
252    out_file = 'out.tif'
253
254    ulx = None
255    psize_x = None
256    separate = 0
257    copy_pct = 0
258    nodata = None
259    create_options = []
260    pre_init = None
261    band_type = None
262    createonly = 0
263   
264    gdal.AllRegister()
265    if argv is None:
266        argv = sys.argv
267    argv = gdal.GeneralCmdLineProcessor( argv )
268    if argv is None:
269        sys.exit( 0 )
270
271    # Parse command line arguments.
272    i = 1
273    while i < len(argv):
274        arg = argv[i]
275
276        if arg == '-o':
277            i = i + 1
278            out_file = argv[i]
279
280        elif arg == '-v':
281            verbose = 1
282
283        elif arg == '-q' or arg == '-quiet':
284            quiet = 1
285
286        elif arg == '-createonly':
287            createonly = 1
288
289        elif arg == '-separate':
290            separate = 1
291
292        elif arg == '-seperate':
293            separate = 1
294
295        elif arg == '-pct':
296            copy_pct = 1
297
298        elif arg == '-ot':
299            i = i + 1
300            band_type = gdal.GetDataTypeByName( argv[i] )
301            if band_type == gdal.GDT_Unknown:
302                print('Unknown GDAL data type: ', argv[i])
303                sys.exit( 1 )
304
305        elif arg == '-init':
306            i = i + 1
307            pre_init = float(argv[i])
308
309        elif arg == '-n':
310            i = i + 1
311            nodata = float(argv[i])
312
313        elif arg == '-f':
314            # for backward compatibility.
315            i = i + 1
316            format = argv[i]
317
318        elif arg == '-of':
319            i = i + 1
320            format = argv[i]
321
322        elif arg == '-co':
323            i = i + 1
324            create_options.append( argv[i] )
325
326        elif arg == '-ps':
327            psize_x = float(argv[i+1])
328            psize_y = -1 * abs(float(argv[i+2]))
329            i = i + 2
330
331        elif arg == '-ul_lr':
332            ulx = float(argv[i+1])
333            uly = float(argv[i+2])
334            lrx = float(argv[i+3])
335            lry = float(argv[i+4])
336            i = i + 4
337
338        elif arg[:1] == '-':
339            print('Unrecognised command option: ', arg)
340            Usage()
341            sys.exit( 1 )
342
343        else:
344            # Expand any possible wildcards from command line arguments
345            f = glob.glob( arg )
346            if len(f) == 0:
347                print('File not found: "%s"' % (str( arg )))
348            names += f # append 1 or more files
349        i = i + 1
350
351    if len(names) == 0:
352        print('No input files selected.')
353        Usage()
354        sys.exit( 1 )
355
356    Driver = gdal.GetDriverByName(format)
357    if Driver is None:
358        print('Format driver %s not found, pick a supported driver.' % format)
359        sys.exit( 1 )
360
361    DriverMD = Driver.GetMetadata()
362    if 'DCAP_CREATE' not in DriverMD:
363        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).' % format)
364        sys.exit( 1 )
365
366    # Collect information on all the source files.
367    file_infos = names_to_fileinfos( names )
368
369    if ulx is None:
370        ulx = file_infos[0].ulx
371        uly = file_infos[0].uly
372        lrx = file_infos[0].lrx
373        lry = file_infos[0].lry
374       
375        for fi in file_infos:
376            ulx = min(ulx, fi.ulx)
377            uly = max(uly, fi.uly)
378            lrx = max(lrx, fi.lrx)
379            lry = min(lry, fi.lry)
380
381    if psize_x is None:
382        psize_x = file_infos[0].geotransform[1]
383        psize_y = file_infos[0].geotransform[5]
384
385    if band_type is None:
386        band_type = file_infos[0].band_type
387
388    # Try opening as an existing file.
389    gdal.PushErrorHandler( 'CPLQuietErrorHandler' )
390    t_fh = gdal.Open( out_file, gdal.GA_Update )
391    gdal.PopErrorHandler()
392   
393    # Create output file if it does not already exist.
394    if t_fh is None:
395        geotransform = [ulx, psize_x, 0, uly, 0, psize_y]
396
397        xsize = int((lrx - ulx) / geotransform[1] + 0.5)
398        ysize = int((lry - uly) / geotransform[5] + 0.5)
399
400        if separate != 0:
401            bands = len(file_infos)
402        else:
403            bands = file_infos[0].bands
404
405        t_fh = Driver.Create( out_file, xsize, ysize, bands,
406                              band_type, create_options )
407        if t_fh is None:
408            print('Creation failed, terminating gdal_merge.')
409            sys.exit( 1 )
410           
411        t_fh.SetGeoTransform( geotransform )
412        t_fh.SetProjection( file_infos[0].projection )
413
414        if copy_pct:
415            t_fh.GetRasterBand(1).SetRasterColorTable(file_infos[0].ct)
416    else:
417        if separate != 0:
418            bands = len(file_infos)
419            if t_fh.RasterCount < bands :
420                print('Existing output file has less bands than the number of input files. You should delete it before. Terminating gdal_merge.')
421                sys.exit( 1 )
422        else:
423            bands = min(file_infos[0].bands,t_fh.RasterCount)
424
425    # Do we need to pre-initialize the whole mosaic file to some value?
426    if pre_init is not None:
427        for i in range(t_fh.RasterCount):
428            t_fh.GetRasterBand(i+1).Fill( pre_init )
429
430    # Copy data from source files into output file.
431    t_band = 1
432
433    if quiet == 0 and verbose == 0:
434        gdal.TermProgress( 0.0 )
435    fi_processed = 0
436   
437    for fi in file_infos:
438        if createonly != 0:
439            continue
440       
441        if verbose != 0:
442            print("")
443            print("Processing file %5d of %5d, %6.3f%% completed." \
444                  % (fi_processed+1,len(file_infos),
445                     fi_processed * 100.0 / len(file_infos)) )
446            fi.report()
447
448        if separate == 0 :
449            for band in range(1, bands+1):
450                fi.copy_into( t_fh, band, band, nodata )
451        else:
452            fi.copy_into( t_fh, 1, t_band, nodata )
453            t_band = t_band+1
454           
455        fi_processed = fi_processed+1
456        if quiet == 0 and verbose == 0:
457            gdal.TermProgress( fi_processed / float(len(file_infos))  )
458   
459    # Force file to be closed.
460    t_fh = None
461
462if __name__ == '__main__':
463    sys.exit(main())
Note: See TracBrowser for help on using the repository browser.