gdal_merge.ipynb

{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#!/usr/bin/env python\n",
      "###############################################################################\n",
      "# $Id$\n",
      "#\n",
      "# Project:  InSAR Peppers\n",
      "# Purpose:  Module to extract data from many rasters into one output.\n",
      "# Author:   Frank Warmerdam, warmerdam@pobox.com\n",
      "#\n",
      "###############################################################################\n",
      "# Copyright (c) 2000, Atlantis Scientific Inc. (www.atlsci.com)\n",
      "# Copyright (c) 2009-2011, Even Rouault <even dot rouault at mines-paris dot org>\n",
      "# \n",
      "# This library is free software; you can redistribute it and/or\n",
      "# modify it under the terms of the GNU Library General Public\n",
      "# License as published by the Free Software Foundation; either\n",
      "# version 2 of the License, or (at your option) any later version.\n",
      "# \n",
      "# This library is distributed in the hope that it will be useful,\n",
      "# but WITHOUT ANY WARRANTY; without even the implied warranty of\n",
      "# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU\n",
      "# Library General Public License for more details.\n",
      "# \n",
      "# You should have received a copy of the GNU Library General Public\n",
      "# License along with this library; if not, write to the\n",
      "# Free Software Foundation, Inc., 59 Temple Place - Suite 330,\n",
      "# Boston, MA 02111-1307, USA.\n",
      "###############################################################################\n",
      "# changes 29Apr2011\n",
      "# if the input image is a multi-band one, use all \n",
      "# the channels in building the stack\n",
      "# anssi.pekkarinen@fao.org\n",
      "\n",
      "\n",
      "try:\n",
      "    from osgeo import gdal\n",
      "except ImportError:\n",
      "    import gdal\n",
      "\n",
      "try:\n",
      "    progress = gdal.TermProgress_nocb\n",
      "except:\n",
      "    progress = gdal.TermProgress\n",
      "\n",
      "\n",
      "import sys\n",
      "import math\n",
      "\n",
      "__version__ = '$id$'[5:-1]\n",
      "verbose = 0\n",
      "quiet = 0\n",
      "\n",
      "\n",
      "# =============================================================================\n",
      "def raster_copy( s_fh, s_xoff, s_yoff, s_xsize, s_ysize, s_band_n,\n",
      "                 t_fh, t_xoff, t_yoff, t_xsize, t_ysize, t_band_n,\n",
      "                 nodata=None ):\n",
      "\n",
      "    if nodata is not None:\n",
      "        return raster_copy_with_nodata(\n",
      "            s_fh, s_xoff, s_yoff, s_xsize, s_ysize, s_band_n,\n",
      "            t_fh, t_xoff, t_yoff, t_xsize, t_ysize, t_band_n,\n",
      "            nodata )\n",
      "\n",
      "    if verbose != 0:\n",
      "        print('Copy %d,%d,%d,%d to %d,%d,%d,%d.' \\\n",
      "              % (s_xoff, s_yoff, s_xsize, s_ysize,\n",
      "             t_xoff, t_yoff, t_xsize, t_ysize ))\n",
      "\n",
      "    s_band = s_fh.GetRasterBand( s_band_n )\n",
      "    t_band = t_fh.GetRasterBand( t_band_n )\n",
      "\n",
      "    data = s_band.ReadRaster( s_xoff, s_yoff, s_xsize, s_ysize,\n",
      "                             t_xsize, t_ysize, t_band.DataType )\n",
      "    t_band.WriteRaster( t_xoff, t_yoff, t_xsize, t_ysize,\n",
      "                        data, t_xsize, t_ysize, t_band.DataType )\n",
      "        \n",
      "\n",
      "    return 0\n",
      "    \n",
      "# =============================================================================\n",
      "def raster_copy_with_nodata( s_fh, s_xoff, s_yoff, s_xsize, s_ysize, s_band_n,\n",
      "                             t_fh, t_xoff, t_yoff, t_xsize, t_ysize, t_band_n,\n",
      "                             nodata ):\n",
      "    try:\n",
      "        import numpy as Numeric\n",
      "    except ImportError:\n",
      "        import Numeric\n",
      "    \n",
      "    if verbose != 0:\n",
      "        print('Copy %d,%d,%d,%d to %d,%d,%d,%d.' \\\n",
      "              % (s_xoff, s_yoff, s_xsize, s_ysize,\n",
      "             t_poff, t_yoff, t_xsize, t_ysize ))\n",
      "\n",
      "    s_band = s_fh.GetRasterBand( s_band_n )\n",
      "    t_band = t_fh.GetRasterBand( t_band_n )\n",
      "\n",
      "    data_src = s_band.ReadAsArray( s_xoff, s_yoff, s_xsize, s_ysize,\n",
      "                                   t_xsize, t_ysize )\n",
      "    data_dst = t_band.ReadAsArray( t_xoff, t_yoff, t_xsize, t_ysize )\n",
      "\n",
      "    nodata_test = Numeric.equal(data_src,nodata)\n",
      "    to_write = Numeric.choose( nodata_test, (data_src, data_dst) )\n",
      "\n",
      "    t_band.WriteArray( to_write, t_xoff, t_yoff )\n",
      "\n",
      "    return 0\n",
      "    \n",
      "# =============================================================================\n",
      "def names_to_fileinfos( names ):\n",
      "    \"\"\"\n",
      "    Translate a list of GDAL filenames, into file_info objects.\n",
      "\n",
      "    names -- list of valid GDAL dataset names.\n",
      "\n",
      "    Returns a list of file_info objects.  There may be less file_info objects\n",
      "    than names if some of the names could not be opened as GDAL files.\n",
      "    \"\"\"\n",
      "    \n",
      "    file_infos = []\n",
      "    for name in names:\n",
      "        fi = file_info()\n",
      "        if fi.init_from_name( name ) == 1:\n",
      "            file_infos.append( fi )\n",
      "\n",
      "    return file_infos\n",
      "\n",
      "# *****************************************************************************\n",
      "class file_info:\n",
      "    \"\"\"A class holding information about a GDAL file.\"\"\"\n",
      "\n",
      "    def init_from_name(self, filename):\n",
      "        \"\"\"\n",
      "        Initialize file_info from filename\n",
      "\n",
      "        filename -- Name of file to read.\n",
      "\n",
      "        Returns 1 on success or 0 if the file can't be opened.\n",
      "        \"\"\"\n",
      "        fh = gdal.Open( filename )\n",
      "        if fh is None:\n",
      "            return 0\n",
      "\n",
      "        self.filename = filename\n",
      "        self.bands = fh.RasterCount\n",
      "        self.xsize = fh.RasterXSize\n",
      "        self.ysize = fh.RasterYSize\n",
      "        self.band_type = fh.GetRasterBand(1).DataType\n",
      "        self.projection = fh.GetProjection()\n",
      "        self.geotransform = fh.GetGeoTransform()\n",
      "        self.ulx = self.geotransform[0]\n",
      "        self.uly = self.geotransform[3]\n",
      "        self.lrx = self.ulx + self.geotransform[1] * self.xsize\n",
      "        self.lry = self.uly + self.geotransform[5] * self.ysize\n",
      "\n",
      "        ct = fh.GetRasterBand(1).GetRasterColorTable()\n",
      "        if ct is not None:\n",
      "            self.ct = ct.Clone()\n",
      "        else:\n",
      "            self.ct = None\n",
      "\n",
      "        return 1\n",
      "\n",
      "    def report( self ):\n",
      "        print('Filename: '+ self.filename)\n",
      "        print('File Size: %dx%dx%d' \\\n",
      "              % (self.xsize, self.ysize, self.bands))\n",
      "        print('Pixel Size: %f x %f' \\\n",
      "              % (self.geotransform[1],self.geotransform[5]))\n",
      "        print('UL:(%f,%f)   LR:(%f,%f)' \\\n",
      "              % (self.ulx,self.uly,self.lrx,self.lry))\n",
      "\n",
      "    def copy_into( self, t_fh, s_band = 1, t_band = 1, nodata_arg=None ):\n",
      "        \"\"\"\n",
      "        Copy this files image into target file.\n",
      "\n",
      "        This method will compute the overlap area of the file_info objects\n",
      "        file, and the target gdal.Dataset object, and copy the image data\n",
      "        for the common window area.  It is assumed that the files are in\n",
      "        a compatible projection ... no checking or warping is done.  However,\n",
      "        if the destination file is a different resolution, or different\n",
      "        image pixel type, the appropriate resampling and conversions will\n",
      "        be done (using normal GDAL promotion/demotion rules).\n",
      "\n",
      "        t_fh -- gdal.Dataset object for the file into which some or all\n",
      "        of this file may be copied.\n",
      "\n",
      "        Returns 1 on success (or if nothing needs to be copied), and zero one\n",
      "        failure.\n",
      "        \"\"\"\n",
      "        t_geotransform = t_fh.GetGeoTransform()\n",
      "        t_ulx = t_geotransform[0]\n",
      "        t_uly = t_geotransform[3]\n",
      "        t_lrx = t_geotransform[0] + t_fh.RasterXSize * t_geotransform[1]\n",
      "        t_lry = t_geotransform[3] + t_fh.RasterYSize * t_geotransform[5]\n",
      "\n",
      "        # figure out intersection region\n",
      "        tgw_ulx = max(t_ulx,self.ulx)\n",
      "        tgw_lrx = min(t_lrx,self.lrx)\n",
      "        if t_geotransform[5] < 0:\n",
      "            tgw_uly = min(t_uly,self.uly)\n",
      "            tgw_lry = max(t_lry,self.lry)\n",
      "        else:\n",
      "            tgw_uly = max(t_uly,self.uly)\n",
      "            tgw_lry = min(t_lry,self.lry)\n",
      "        \n",
      "        # do they even intersect?\n",
      "        if tgw_ulx >= tgw_lrx:\n",
      "            return 1\n",
      "        if t_geotransform[5] < 0 and tgw_uly <= tgw_lry:\n",
      "            return 1\n",
      "        if t_geotransform[5] > 0 and tgw_uly >= tgw_lry:\n",
      "            return 1\n",
      "            \n",
      "        # compute target window in pixel coordinates.\n",
      "        tw_xoff = int((tgw_ulx - t_geotransform[0]) / t_geotransform[1] + 0.1)\n",
      "        tw_yoff = int((tgw_uly - t_geotransform[3]) / t_geotransform[5] + 0.1)\n",
      "        tw_xsize = int((tgw_lrx - t_geotransform[0])/t_geotransform[1] + 0.5) \\\n",
      "                   - tw_xoff\n",
      "        tw_ysize = int((tgw_lry - t_geotransform[3])/t_geotransform[5] + 0.5) \\\n",
      "                   - tw_yoff\n",
      "\n",
      "        if tw_xsize < 1 or tw_ysize < 1:\n",
      "            return 1\n",
      "\n",
      "        # Compute source window in pixel coordinates.\n",
      "        sw_xoff = int((tgw_ulx - self.geotransform[0]) / self.geotransform[1])\n",
      "        sw_yoff = int((tgw_uly - self.geotransform[3]) / self.geotransform[5])\n",
      "        sw_xsize = int((tgw_lrx - self.geotransform[0]) \\\n",
      "                       / self.geotransform[1] + 0.5) - sw_xoff\n",
      "        sw_ysize = int((tgw_lry - self.geotransform[3]) \\\n",
      "                       / self.geotransform[5] + 0.5) - sw_yoff\n",
      "\n",
      "        if sw_xsize < 1 or sw_ysize < 1:\n",
      "            return 1\n",
      "\n",
      "        # Open the source file, and copy the selected region.\n",
      "        s_fh = gdal.Open( self.filename )\n",
      "\n",
      "        return \\\n",
      "            raster_copy( s_fh, sw_xoff, sw_yoff, sw_xsize, sw_ysize, s_band,\n",
      "                         t_fh, tw_xoff, tw_yoff, tw_xsize, tw_ysize, t_band,\n",
      "                         nodata_arg )\n",
      "\n",
      "\n",
      "# =============================================================================\n",
      "def Usage():\n",
      "    print('Usage: gdal_merge.py [-o out_filename] [-of out_format] [-co NAME=VALUE]*')\n",
      "    print('                     [-ps pixelsize_x pixelsize_y] [-tap] [-separate] [-q] [-v] [-pct]')\n",
      "    print('                     [-ul_lr ulx uly lrx lry] [-init \"value [value...]\"]')\n",
      "    print('                     [-n nodata_value] [-a_nodata output_nodata_value]')\n",
      "    print('                     [-ot datatype] [-createonly] input_files')\n",
      "    print('                     [--help-general]')\n",
      "    print('')\n",
      "\n",
      "# =============================================================================\n",
      "#\n",
      "# Program mainline.\n",
      "#\n",
      "\n",
      "def main( argv=None ):\n",
      "\n",
      "    global verbose, quiet\n",
      "    verbose = 0\n",
      "    quiet = 0\n",
      "    names = []\n",
      "    format = 'GTiff'\n",
      "    out_file = 'out.tif'\n",
      "\n",
      "    ulx = None\n",
      "    psize_x = None\n",
      "    separate = 0\n",
      "    copy_pct = 0\n",
      "    nodata = None\n",
      "    a_nodata = None\n",
      "    create_options = []\n",
      "    pre_init = []\n",
      "    band_type = None\n",
      "    createonly = 0\n",
      "    bTargetAlignedPixels = False\n",
      "    \n",
      "    gdal.AllRegister()\n",
      "    if argv is None:\n",
      "        argv = sys.argv\n",
      "    argv = gdal.GeneralCmdLineProcessor( argv )\n",
      "    if argv is None:\n",
      "        sys.exit( 0 )\n",
      "\n",
      "    # Parse command line arguments.\n",
      "    i = 1\n",
      "    while i < len(argv):\n",
      "        arg = argv[i]\n",
      "\n",
      "        if arg == '-o':\n",
      "            i = i + 1\n",
      "            out_file = argv[i]\n",
      "\n",
      "        elif arg == '-v':\n",
      "            verbose = 1\n",
      "\n",
      "        elif arg == '-q' or arg == '-quiet':\n",
      "            quiet = 1\n",
      "\n",
      "        elif arg == '-createonly':\n",
      "            createonly = 1\n",
      "\n",
      "        elif arg == '-separate':\n",
      "            separate = 1\n",
      "\n",
      "        elif arg == '-seperate':\n",
      "            separate = 1\n",
      "\n",
      "        elif arg == '-pct':\n",
      "            copy_pct = 1\n",
      "\n",
      "        elif arg == '-ot':\n",
      "            i = i + 1\n",
      "            band_type = gdal.GetDataTypeByName( argv[i] )\n",
      "            if band_type == gdal.GDT_Unknown:\n",
      "                print('Unknown GDAL data type: %s' % argv[i])\n",
      "                sys.exit( 1 )\n",
      "\n",
      "        elif arg == '-init':\n",
      "            i = i + 1\n",
      "            str_pre_init = argv[i].split()\n",
      "            for x in str_pre_init:\n",
      "                pre_init.append(float(x))\n",
      "\n",
      "        elif arg == '-n':\n",
      "            i = i + 1\n",
      "            nodata = float(argv[i])\n",
      "\n",
      "        elif arg == '-a_nodata':\n",
      "            i = i + 1\n",
      "            a_nodata = float(argv[i])\n",
      "\n",
      "        elif arg == '-f':\n",
      "            # for backward compatibility.\n",
      "            i = i + 1\n",
      "            format = argv[i]\n",
      "\n",
      "        elif arg == '-of':\n",
      "            i = i + 1\n",
      "            format = argv[i]\n",
      "\n",
      "        elif arg == '-co':\n",
      "            i = i + 1\n",
      "            create_options.append( argv[i] )\n",
      "\n",
      "        elif arg == '-ps':\n",
      "            psize_x = float(argv[i+1])\n",
      "            psize_y = -1 * abs(float(argv[i+2]))\n",
      "            i = i + 2\n",
      "\n",
      "        elif arg == '-tap':\n",
      "            bTargetAlignedPixels = True\n",
      "\n",
      "        elif arg == '-ul_lr':\n",
      "            ulx = float(argv[i+1])\n",
      "            uly = float(argv[i+2])\n",
      "            lrx = float(argv[i+3])\n",
      "            lry = float(argv[i+4])\n",
      "            i = i + 4\n",
      "\n",
      "        elif arg[:1] == '-':\n",
      "            print('Unrecognised command option: %s' % arg)\n",
      "            Usage()\n",
      "            sys.exit( 1 )\n",
      "\n",
      "        else:\n",
      "            names.append(arg)\n",
      "\n",
      "        i = i + 1\n",
      "\n",
      "    if len(names) == 0:\n",
      "        print('No input files selected.')\n",
      "        Usage()\n",
      "        sys.exit( 1 )\n",
      "\n",
      "    Driver = gdal.GetDriverByName(format)\n",
      "    if Driver is None:\n",
      "        print('Format driver %s not found, pick a supported driver.' % format)\n",
      "        sys.exit( 1 )\n",
      "\n",
      "    DriverMD = Driver.GetMetadata()\n",
      "    if 'DCAP_CREATE' not in DriverMD:\n",
      "        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)\n",
      "        sys.exit( 1 )\n",
      "\n",
      "    # Collect information on all the source files.\n",
      "    file_infos = names_to_fileinfos( names )\n",
      "\n",
      "    if ulx is None:\n",
      "        ulx = file_infos[0].ulx\n",
      "        uly = file_infos[0].uly\n",
      "        lrx = file_infos[0].lrx\n",
      "        lry = file_infos[0].lry\n",
      "        \n",
      "        for fi in file_infos:\n",
      "            ulx = min(ulx, fi.ulx)\n",
      "            uly = max(uly, fi.uly)\n",
      "            lrx = max(lrx, fi.lrx)\n",
      "            lry = min(lry, fi.lry)\n",
      "\n",
      "    if psize_x is None:\n",
      "        psize_x = file_infos[0].geotransform[1]\n",
      "        psize_y = file_infos[0].geotransform[5]\n",
      "\n",
      "    if band_type is None:\n",
      "        band_type = file_infos[0].band_type\n",
      "\n",
      "    # Try opening as an existing file.\n",
      "    gdal.PushErrorHandler( 'CPLQuietErrorHandler' )\n",
      "    t_fh = gdal.Open( out_file, gdal.GA_Update )\n",
      "    gdal.PopErrorHandler()\n",
      "    \n",
      "    # Create output file if it does not already exist.\n",
      "    if t_fh is None:\n",
      "    \n",
      "        if bTargetAlignedPixels:\n",
      "            ulx = math.floor(ulx / psize_x) * psize_x\n",
      "            lrx = math.ceil(lrx / psize_x) * psize_x\n",
      "            lry = math.floor(lry / -psize_y) * -psize_y\n",
      "            uly = math.ceil(uly / -psize_y) * -psize_y\n",
      "    \n",
      "        geotransform = [ulx, psize_x, 0, uly, 0, psize_y]\n",
      "\n",
      "        xsize = int((lrx - ulx) / geotransform[1] + 0.5)\n",
      "        ysize = int((lry - uly) / geotransform[5] + 0.5)\n",
      "\n",
      "\n",
      "        if separate != 0:\n",
      "            bands=0\n",
      "\n",
      "            for fi in file_infos:\n",
      "                bands=bands + fi.bands\n",
      "        else:\n",
      "            bands = file_infos[0].bands\n",
      "\n",
      "\n",
      "        t_fh = Driver.Create( out_file, xsize, ysize, bands,\n",
      "                              band_type, create_options )\n",
      "        if t_fh is None:\n",
      "            print('Creation failed, terminating gdal_merge.')\n",
      "            sys.exit( 1 )\n",
      "            \n",
      "        t_fh.SetGeoTransform( geotransform )\n",
      "        t_fh.SetProjection( file_infos[0].projection )\n",
      "\n",
      "        if copy_pct:\n",
      "            t_fh.GetRasterBand(1).SetRasterColorTable(file_infos[0].ct)\n",
      "    else:\n",
      "        if separate != 0:\n",
      "            bands=0\n",
      "            for fi in file_infos:\n",
      "                bands=bands + fi.bands            \n",
      "            if t_fh.RasterCount < bands :\n",
      "                print('Existing output file has less bands than the input files. You should delete it before. Terminating gdal_merge.')\n",
      "                sys.exit( 1 )\n",
      "        else:\n",
      "            bands = min(file_infos[0].bands,t_fh.RasterCount)\n",
      "\n",
      "    # Do we need to set nodata value ?\n",
      "    if a_nodata is not None:\n",
      "        for i in range(t_fh.RasterCount):\n",
      "            t_fh.GetRasterBand(i+1).SetNoDataValue(a_nodata)\n",
      "\n",
      "    # Do we need to pre-initialize the whole mosaic file to some value?\n",
      "    if pre_init is not None:\n",
      "        if t_fh.RasterCount <= len(pre_init):\n",
      "            for i in range(t_fh.RasterCount):\n",
      "                t_fh.GetRasterBand(i+1).Fill( pre_init[i] )\n",
      "        elif len(pre_init) == 1:\n",
      "            for i in range(t_fh.RasterCount):\n",
      "                t_fh.GetRasterBand(i+1).Fill( pre_init[0] )\n",
      "\n",
      "    # Copy data from source files into output file.\n",
      "    t_band = 1\n",
      "\n",
      "    if quiet == 0 and verbose == 0:\n",
      "        progress( 0.0 )\n",
      "    fi_processed = 0\n",
      "    \n",
      "    for fi in file_infos:\n",
      "        if createonly != 0:\n",
      "            continue\n",
      "        \n",
      "        if verbose != 0:\n",
      "            print(\"\")\n",
      "            print(\"Processing file %5d of %5d, %6.3f%% completed.\" \\\n",
      "                  % (fi_processed+1,len(file_infos),\n",
      "                     fi_processed * 100.0 / len(file_infos)) )\n",
      "            fi.report()\n",
      "\n",
      "        if separate == 0 :\n",
      "            for band in range(1, bands+1):\n",
      "                fi.copy_into( t_fh, band, band, nodata )\n",
      "        else:\n",
      "            for band in range(1, fi.bands+1):\n",
      "                fi.copy_into( t_fh, band, t_band, nodata )\n",
      "                t_band = t_band+1\n",
      "            \n",
      "        fi_processed = fi_processed+1\n",
      "        if quiet == 0 and verbose == 0:\n",
      "            progress( fi_processed / float(len(file_infos))  )\n",
      "    \n",
      "    # Force file to be closed.\n",
      "    t_fh = None\n",
      "\n",
      "if __name__ == '__main__':\n",
      "    sys.exit(main())"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}