#!/usr/bin/env python # -*- coding: utf-8 -*- # Copyright (c) 2000, Atlantis Scientific Inc. (www.atlsci.com) # Copyright (C) 2005 Gabriel Ebner <ge@gabrielebner.at> # # 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 try: from osgeo import gdal except ImportError: import gdal import sys # ============================================================================= 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: """A class holding information about a GDAL file.""" 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() != 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 ('Unrecognised command option: ', arg) Usage() sys.exit( 1 ) else: names.append( arg ) i = i + 1 if len(names) == 0: print ('No input files selected.') Usage() sys.exit( 1 ) # Collect information on all the source files. file_infos = names_to_fileinfos( names ) if len(file_infos) == 0: 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 len(projection) > 0: 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 accout' % (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 != 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')