gdal2cesium.py

#!/usr/bin/env python
#******************************************************************************
#  $Id: gdal2cesium.py 2014-10-01 12:01:23Z $
# 
# Project:  Cesium terrain generator for GDAL raster formats - S.I.T. Comune di Prato (Italy)
# Support:  Gis3w s.a.s. (http://gis3w.it)
# Purpose:  Convert a raster into a heightmap terrain for Cesium 3D Javascript library
#           - generate a global geodetic TMS data structure
#           - tiles are generated according to the Cesium heightmap binary format v1.0 (http://cesiumjs.org/data-and-assets/terrain/formats/heightmap-1.0.html)
#           - the max zoom level is calculated on the base of the raster horizontal resolution
#           - zoom levels up to the 0 level are always created to complete the parent-child relationships required by the Cesium format
# Author:   Giovanni Allegri (http://giovanniallegri.it, http://gis3w.it)
#
###############################################################################
# Copyright (c) 2014, S.I.T. Comune di Prato (Italy)

#  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.                                   
# 
#  Permission is hereby granted, free of charge, to any person obtaining a
#  copy of this software and associated documentation files (the "Software"),
#  to deal in the Software without restriction, including without limitation
#  the rights to use, copy, modify, merge, publish, distribute, sublicense,
#  and/or sell copies of the Software, and to permit persons to whom the
#  Software is furnished to do so, subject to the following conditions:
# 
#  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
#  OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
#  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
#  THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
#  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
#  FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
#  DEALINGS IN THE SOFTWARE.
#******************************************************************************

import os,sys,math,glob,struct,shutil
import subprocess

try:
    from osgeo import gdal
    from osgeo import osr
    from osgeo import gdal_array
except:
    import gdal
    print('You are using "old gen" bindings. gdal2cesium needs "new gen" bindings.')
    sys.exit(1)

from shapely.geometry import mapping, Polygon, LineString
from osgeo import ogr


try:
    import numpy
    import osgeo.gdal_array as gdalarray
except:
    print('gdal2cesium needs Numpy.')
    sys.exit(1)

MAXZOOMLEVEL = 32
resampling_list = ('average','near','bilinear','cubic','cubicspline','lanczos')

def makepoly(ulx,uly,lrx,lry):
    return Polygon([(ulx, uly), (lrx, uly), (lrx, lry), (ulx, lry), (ulx, uly)])
    
def makeline(ulx,uly,lrx,lry):
    return LineString([(ulx, uly), (lrx, uly), (lrx, lry), (ulx, lry), (ulx, uly)])
	
def splitpath(path):
    parts=[]
    (path, tail)=os.path.split( path)
    while path and tail:
         parts.append( tail)
         (path,tail)=os.path.split(path)
    parts.append( os.path.join(path,tail) )
    return map( os.path.normpath, parts)[::-1]

class GlobalGeodetic(object):
    def __init__(self, tileSize = 256):
        self.tileSize = tileSize

    def LatLonToPixels(self, lat, lon, zoom):
        "Converts lat/lon to pixel coordinates in given zoom of the EPSG:4326 pyramid"

        res = 180.0 / self.tileSize / 2**zoom
        px = (180 + lat) / res
        py = (90 + lon) / res
        return px, py

    def PixelsToTile(self, px, py):
        "Returns coordinates of the tile covering region in pixel coordinates"

        tx = int( math.ceil( px / float(self.tileSize) ) - 1 )
        ty = int( math.ceil( py / float(self.tileSize) ) - 1 )
        return tx, ty

    def LatLonToTile(self, lat, lon, zoom):
        "Returns the tile for zoom which covers given lat/lon coordinates"

        px, py = self.LatLonToPixels( lat, lon, zoom)
        return self.PixelsToTile(px,py)

    def Resolution(self, zoom ):
        "Resolution (arc/pixel) for given zoom level (measured at Equator)"

        return 180.0 / self.tileSize / 2**zoom
    def ZoomForPixelSize(self, pixelSize ):
        "Maximal scaledown zoom of the pyramid closest to the pixelSize."

        for i in range(MAXZOOMLEVEL):
            if pixelSize > self.Resolution(i):
                if i!=0:
                    return i-1
                else:
                    return 0 # We don't want to scale up

    def TileBounds(self, tx, ty, zoom):
        "Returns bounds of the given tile"
        res = 180.0 / self.tileSize / 2**zoom
        return (
            tx*self.tileSize*res - 180,
            ty*self.tileSize*res - 90,
            (tx+1)*self.tileSize*res - 180,
            (ty+1)*self.tileSize*res - 90
        )
        
    def TileBoundsForTileSize(self, tx, ty, zoom, extrapixels):
        res = 180.0 / self.tileSize / 2**zoom
        # we have to calculate a wider bound to consider the overlapping pixel according to Cesium format
        extrafactor = res*extrapixels
        return (
            tx*self.tileSize*res - 180,
            (ty*self.tileSize*res) - extrafactor - 90,
            ((tx+1)*self.tileSize*res - 180) + extrafactor,
            (ty+1)*self.tileSize*res - 90
        )

    def TileLatLonBounds(self, tx, ty, zoom):
        "Returns bounds of the given tile in the SWNE form"
        b = self.TileBounds(tx, ty, zoom)
        return (b[1],b[0],b[3],b[2])

import fnmatch
class PostProcessor(object):
    def __init__(self,outPath = "./",tmpOutPath = "./tmp"):
        self.rootPath = tmpOutPath
        self.pattern = '*.terrain'
        self.processedPath = outPath
        self.rtype = numpy.dtype('int16')
    
    def walk_tiles(self,folder = '.'):
        for root, _, files in os.walk(self.rootPath):
            for filename in fnmatch.filter(files, self.pattern):
                yield( os.path.join(root, filename))
        
    def get_tiles(self):
        terrains = []
        for terrain in self.walk_tiles():
            terrains.append(terrain)
        return terrains
    
    def extract_data(self,fin,rb):
        data = numpy.fromfile(fin,dtype=self.rtype,count=4096)
        data_mat = data.reshape((64,64))
        if rb == 'r':
            data_slice = data_mat[:,0] # left col
        else:
            data_slice = data_mat[0] # top row
        return data_slice
        
    def augment_tile(self,fin,data_slice_r,data_slice_b,pixel_br):
        data_complete = numpy.fromfile(fin,dtype=self.rtype,count=4097)
        data = data_complete[:4096]
        maskbytes = data_complete[4096]
        data_mat = data.reshape((64,64))
    
        data_slice_b_br = numpy.c_[[data_slice_b],[pixel_br]] # add bottom right pixel to bottom row as new col
        data_mat_r = numpy.c_[data_mat,data_slice_r] # add right col to data
        data_mat_brpr = numpy.r_[data_mat_r,data_slice_b_br] # add bottom row to data
        
        return data_mat_brpr,maskbytes
    
    def write_tile(self,tilename,new_tile,maskbytes):
        tilepath = os.path.join(self.processedPath,tilename)
        if not os.path.exists(os.path.dirname(tilepath)):
            os.makedirs(os.path.dirname(tilepath))
         
         
        tilearrayint = new_tile.astype(numpy.int16)
        data = tilearrayint.flatten()
        data_with_mask = numpy.append(data,maskbytes)
        data_with_mask.tofile(tilepath)
        
    def run(self):
        for terrain in self.get_tiles():
            pathparts = splitpath(terrain)
            idx = len(pathparts)
            root = os.path.join(*pathparts[:idx-3])
            y = int(pathparts[idx-1].split(".")[0])
            x = int(pathparts[idx-2])
            z = int(pathparts[idx-3])
            right_tile = os.path.join(root,str(z),str(x+1),"%s.terrain" % y)
            bottom_tile = os.path.join(root,str(z),str(x),"%s.terrain" % str(y-1))
            bottom_right_tile = os.path.join(root,str(z),str(x+1),"%s.terrain" % str(y-1))
            
            if os.path.exists(right_tile):
                with open(right_tile, 'rb') as right_tile_f:
                    data_slice_r = self.extract_data(right_tile_f,'r')
            else:
                data_slice_r = numpy.empty(64)
                data_slice_r.fill(5000)
            
            if os.path.exists(bottom_tile):
                with open(bottom_tile, 'rb') as bottom_tile_f:
                    data_slice_b = self.extract_data(bottom_tile_f,'t')
            else:
                data_slice_b = numpy.empty(64)
                data_slice_b.fill(5000)
            
            if os.path.exists(bottom_right_tile):
                with open(bottom_right_tile, 'rb') as bottom_right_tile_f:
                    data = numpy.fromfile(bottom_right_tile_f,dtype=self.rtype,count=1)
                    pixel_br = data[0]
            else:
                pixel_br = 5000
        
            with open(terrain, 'rb') as terrain_f:
                new_tile,maskbytes = self.augment_tile(terrain_f,data_slice_r,data_slice_b,pixel_br)
                tilename = os.path.join(*pathparts[idx-3:idx])
                self.write_tile(tilename,new_tile,maskbytes)

class GDAL2Cesium(object):
    # -------------------------------------------------------------------------    
    def process(self):
        for inumpyut_file in self.inumpyuts:
            self.inumpyut = inumpyut_file
            self.pre_process_inumpyut(inumpyut_file)
        
        self.merge_inumpyuts_data()
        
        self.make_tiles()
        
        print """Running post processing"""
        pp = PostProcessor(self.output,self.tmpoutput)
        pp.run()
        print """Post processing terminated"""
        shutil.rmtree(self.tmpoutput)
        
        
    def merge_inumpyuts_data(self):
        # Merge tminmax. We will use the extent containing all the layers for the lower zooms and only the higher resolution layer for the highest zooms
        global_tminmax = []
        for _inumpyut,inumpyut_data in self.inumpyuts_data.iteritems():
            #print "Inumpyut: %s" % _inumpyut
            minz = inumpyut_data[0]
            maxz = inumpyut_data[1]
            tminmax = inumpyut_data[2]
            for tz,tminmax_values in enumerate(tminmax):
                if (self.user_tminz is not None and tz < self.user_tminz) or (self.user_tmaxz is not None and tz > self.user_tmaxz):
                    continue
                if tz <= maxz:
                    #print "  tz: %s, tminmax: %s" % (tz,tminmax_values)
                    if len(global_tminmax)<=tz:
                        global_tminmax.append(list(tminmax_values))
                    else:
                        tminx = tminmax_values[0]
                        tminy = tminmax_values[1]
                        tmaxx = tminmax_values[2]
                        tmaxy = tminmax_values[3]
                        if tminx < global_tminmax[tz][0]:
                            global_tminmax[tz][0] = tminx
                        if tminy < global_tminmax[tz][1]:
                            global_tminmax[tz][1] = tminy
                        if tmaxx > global_tminmax[tz][2]:
                            global_tminmax[tz][2] = tmaxx
                        if tmaxy > global_tminmax[tz][3]:
                            global_tminmax[tz][3] = tmaxy
                            
        self.tminmax = global_tminmax
        
        # Split zooms in zoom ranges based on resolutions (to build the related vrt files)
        for _inumpyut,inumpyut_data in self.inumpyuts_data.iteritems():
            minz = inumpyut_data[0]
            maxz = inumpyut_data[1]
            if self.tminz is None or minz < self.tminz:
                self.tminz = minz
            if self.tmaxz is None or maxz > self.tmaxz:
                self.tmaxz = maxz
            for zoom in range(minz,maxz+1):
                if (self.user_tminz is not None and tz < self.user_tminz) or (self.user_tmaxz is not None and tz > self.user_tmaxz):
                    continue
                if self.zoom_resolutions.get(zoom) is None:
                    self.zoom_resolutions[zoom] = (inumpyut_data[3],inumpyut_data[4])
                else:
                    # the worst resolution is assigned to the common zoom levels (we check only resx, because resy will be consequently correlated)
                    if self.zoom_resolutions[zoom][0] < inumpyut_data[3]:
                        self.zoom_resolutions[zoom] = (inumpyut_data[3],inumpyut_data[4])
        '''print "MERGED"
        for tz,tminmax_values in enumerate(self.global_tminmax): 
            print "  tz: %s, tminmax: %s" % (tz,tminmax_values)
        '''
        
    # -------------------------------------------------------------------------
    def error(self, msg, details = "" ):
        """Print an error message and stop the processing"""
        if details:
            self.parser.error(msg + "\n\n" + details)
        else:
            self.parser.error(msg)
        exit(1)

    # -------------------------------------------------------------------------
    def progressbar(self, complete = 0.0):
        """Print progressbar for float value 0..1"""

        gdal.TermProgress_nocb(complete)

    # -------------------------------------------------------------------------
    def stop(self):
        """Stop the rendering immediately"""
        self.stopped = True

    # -------------------------------------------------------------------------
    def __init__(self, arguments ):
        """Constructor function - initialization"""
        try:
            subprocess.call(["gdalbuildvrt","--help"])
        except:
            print "gdalbuildvrt is required to run gdal2cesium in multi inumpyuts mode"
            exit(1)
        
        self.stopped = False
        self.multi_suffix = ''
        self.inumpyut = None
        self.default_base_output = 'tiles'
        self.min_tile_tz = None
        self.inumpyuts_data = {}
        self.inumpyuts_files_or_vrt = []
        self.vrts = {}
        self.tminmax = None
        self.zoom_resolutions = {}
        self.tminz = None
        self.tmaxz = None
        
        gdal.AllRegister()
        self.mem_drv = gdal.GetDriverByName( 'MEM' )
        self.geodetic = GlobalGeodetic()

        # Tile format
        self.tilesize = 64
        self.tileext = 'terrain'
        
        self.epsg4326 = "EPSG:4326"
        
        self.tilelayer = None

        self.scaledquery = True
        # How big should be query window be for scaling down
        # Later on reset according the chosen resampling algorightm
        self.querysize = 4 * self.tilesize
        
        # pixel overlap between tiles according to Ceiusm heightmap format
        self.extrapixels = 0 

        # RUN THE ARGUMENT PARSER:
        self.optparse_init()
        self.options, self.args = self.parser.parse_args(args=arguments)
        self.options.srcnodata = None
        if not self.args:
            self.error("No inumpyut file specified")

        # POSTPROCESSING OF PARSED ARGUMENTS:
        # Workaround for old versions of GDAL
        try:
            if (self.options.verbose and self.options.resampling == 'near') or gdal.TermProgress_nocb:
                pass
        except:
            self.error("This version of GDAL is not supported. Please upgrade to 1.6+.")
            #,"You can try run crippled version of gdal2tiles with parameters: -v -r 'near'")
        
        self.inumpyuts = [i for i in self.args]

        # Default values for not given options
        if self.options.output:
            self.output = self.options.output
        else:
            if len(self.inumpyuts)>0:
                self.multi_suffix = '_multi'
            self.output = os.path.join(self.default_base_output,os.path.basename( self.inumpyuts[0] ).split('.')[0]+self.multi_suffix)
            self.options.title = os.path.basename( self.inumpyuts[0]+self.multi_suffix )
        self.tmpoutput = os.path.join(self.output,'tmp')

        # Supported options
        self.resampling = None

        if self.options.resampling == 'average':
            try:
                if gdal.RegenerateOverview:
                    pass
            except:
                self.error("'average' resampling algorithm is not available.", "Please use -r 'near' argument or upgrade to newer version of GDAL.")
        elif self.options.resampling == 'near':
            self.resampling = gdal.GRA_NearestNeighbour
            self.querysize = self.tilesize
        elif self.options.resampling == 'bilinear':
            self.resampling = gdal.GRA_Bilinear
            self.querysize = self.tilesize * 2
        elif self.options.resampling == 'cubic':
            self.resampling = gdal.GRA_Cubic
        elif self.options.resampling == 'cubicspline':
            self.resampling = gdal.GRA_CubicSpline
        elif self.options.resampling == 'lanczos':
            self.resampling = gdal.GRA_Lanczos

        # User specified zoom levels
        self.user_tminz = None
        self.user_tmaxz = None
        if self.options.zoom:
            minmax = self.options.zoom.split('-',1)
            minmax.extend([''])
            min, max = minmax[:2]
            self.user_tminz = int(min)
            if max:
                self.user_tmaxz = int(max)
            else:
                self.user_tmaxz = int(min) 

        # Output the results
        if self.options.verbose:
            print("Options:", self.options)
            print("Inumpyut:", self.inumpyuts[0]+self.multi_suffix)
            print("Output:", self.output)
            print("Cache: %s MB" % (gdal.GetCacheMax() / 1024 / 1024))
            print('')

    # -------------------------------------------------------------------------
    def optparse_init(self):
        """Prepare the option parser for inumpyut (argv)"""
        from optparse import OptionParser, OptionGroup
        usage = "Usage: %prog [options] inumpyut_file(s)"
        p = OptionParser(usage, version="%prog ")
        
        p.add_option("-s", "--s_srs", dest="s_srs",
                          help="Define inumpyut raster CRS (eg EPSG:3003)")
        p.add_option('-z', '--zoom', dest="zoom",
                          help="Zoom levels to render (format:'2-5' or '10').")
        p.add_option("-r", "--resampling", dest="resampling", type='choice', choices=resampling_list,
                        help="Resampling method (%s) - default 'average'" % ",".join(resampling_list))
        p.add_option('-e', '--resume', dest="resume", action="store_true",
                          help="Resume mode. Generate only missing files.")        
        p.add_option("-v", "--verbose",
                          action="store_true", dest="verbose",
                          help="Print status messages to stdout")
        p.add_option("-o", "--o_dir",dest="output",
                          help="Root output directory")
        p.add_option("-i", "--index",dest="createtileindexshp",action="store_true",default=False,
                          help="Create the shapefile of tiles index (True or False)")
        p.add_option("-k", "--keep",dest="keepfiles",action="store_true",default=False,
                          help="Keep temporary files reated by gdal2cesium")
              
        p.set_defaults(resume=False,verbose=False,resampling='average')
        self.parser = p
                
    # -------------------------------------------------------------------------
    def pre_process_inumpyut(self,_inumpyut):
        """Initialization of the inumpyut raster, reprojection if necessary"""
        print "Processing: %s" % _inumpyut
        
        inumpyut_or_vrt = _inumpyut
        
        if not self.mem_drv:
            raise Exception("The 'MEM' driver was not found, is it available in this GDAL build?")

        # Open the inumpyut file
        if self.inumpyut:
            in_ds = gdal.Open(_inumpyut, gdal.GA_ReadOnly)
        else:
            raise Exception("No inumpyut file was specified")

        if self.options.verbose:
            print("Inumpyut file:", "( %sP x %sL - %s bands)" % (self.in_ds.RasterXSize, self.in_ds.RasterYSize, self.in_ds.RasterCount))

        if not in_ds:
            # Note: GDAL prints the ERROR message too
            self.error("It is not possible to open the inumpyut file '%s'." % _inumpyut )

        # Read metadata from the inumpyut file
        if in_ds.RasterCount == 0:
            self.error( "Inumpyut file '%s' has no raster band" % _inumpyut )

        if in_ds.GetRasterBand(1).GetRasterColorTable():
            # TODO: Process directly paletted dataset by generating VRT in memory
            self.error( "Please convert this file to RGB/RGBA and run gdal2tiles on the result.",
            """From paletted file you can create RGBA file (temp.vrt) by:
gdal_translate -of vrt -expand rgba %s temp.vrt
then run:
gdal2tiles temp.vrt""" % _inumpyut )

        # Get NODATA value
        in_nodata = []
        for i in range(1, in_ds.RasterCount+1):
            if in_ds.GetRasterBand(i).GetNoDataValue() != None:
                ndata = in_ds.GetRasterBand(i).GetNoDataValue()
                if math.isnan(ndata):
                    ndata = 'none'
                in_nodata.append( ndata )
        if self.options.srcnodata:
            nds = list(map( float, self.options.srcnodata.split(',')))
            if len(nds) < in_ds.RasterCount:
                in_nodata = (nds * in_ds.RasterCount)[:in_ds.RasterCount]
            else:
                in_nodata = nds

        if self.options.verbose:
            print("NODATA: %s" % in_nodata)

        #
        # Here we should have RGBA inumpyut dataset opened in in_ds
        #
        if self.options.verbose:
            print("Preprocessed file:", "( %sP x %sL - %s bands)" % (in_ds.RasterXSize, in_ds.RasterYSize, in_ds.RasterCount))

        # Spatial Reference System of the inumpyut raster
        self.in_srs = None

        if self.options.s_srs:
            self.in_srs = osr.SpatialReference()
            self.in_srs.SetFromUserInumpyut(self.options.s_srs)
            self.in_srs_wkt = self.in_srs.ExportToWkt()
        else:
            self.in_srs_wkt = in_ds.GetProjection()
            if not self.in_srs_wkt and in_ds.GetGCPCount() != 0:
                self.in_srs_wkt = in_ds.GetGCPProjection()
            if self.in_srs_wkt:
                self.in_srs = osr.SpatialReference()
                self.in_srs.ImportFromWkt(self.in_srs_wkt)

        # Spatial Reference System of tiles
        self.out_srs = osr.SpatialReference()
        self.out_srs.ImportFromEPSG(4326)

        # Are the reference systems the same? Reproject if necessary.
        out_ds = None
                              
        if (in_ds.GetGeoTransform() == (0.0, 1.0, 0.0, 0.0, 0.0, 1.0)) and (in_ds.GetGCPCount() == 0):
            self.error("There is no georeference - neither affine transformation (worldfile) nor GCPs. You can generate only 'raster' profile tiles.",
            "Either gdal2tiles with parameter -p 'raster' or use another GIS software for georeference e.g. gdal_transform -gcp / -a_ullr / -a_srs")
        
        in_srs_code = self.in_srs.GetAttrValue("AUTHORITY", 0)
        in_ds_srs = osr.SpatialReference()
        res = in_ds_srs.ImportFromWkt(in_ds.GetProjection())
        
        if res != 0 and in_srs_code is None:
                print "ERROR! The inumpyut file %s has no SRS associated and no SRS has been defined in inumpyut (-s parameter)" % _inumpyut
                exit(1)
 
        if self.in_srs:
            if in_ds_srs.ExportToProj4() != self.out_srs.ExportToProj4():
                if (self.in_srs.ExportToProj4() != self.out_srs.ExportToProj4()) or (in_ds.GetGCPCount() != 0):
                    print "WARNING! Inumpyut file %s has a SR different from EPSG:4326 (WGS84). This can make the processing significantly slow." % _inumpyut
                    # Generation of VRT dataset in tile projection, default 'nearest neighbour' warping
                    out_ds = gdal.AutoCreateWarpedVRT( in_ds, self.in_srs_wkt, self.out_srs.ExportToWkt() )

                    # TODO: HIGH PRIORITY: Correction of AutoCreateWarpedVRT according the max zoomlevel for correct direct warping!!!

                    if self.options.verbose:
                        print("Warping of the raster by AutoCreateWarpedVRT (result saved into 'tiles.vrt')")
                    out_ds.GetDriver().CreateCopy("%s.vrt" % _inumpyut, out_ds)
                    inumpyut_or_vrt = "%s.vrt" % _inumpyut

                    # Note: self.in_srs and self.in_srs_wkt contain still the non-warped reference system!!!

        else:
            self.error("Inumpyut file has unknown SRS.", "Use --s_srs ESPG:xyz (or similar) to provide source reference system." )

        if out_ds and self.options.verbose:
            print("Projected file:", "tiles.vrt", "( %sP x %sL - %s bands)" % (out_ds.RasterXSize, out_ds.RasterYSize, out_ds.RasterCount))
        
        if not out_ds:
            out_ds = in_ds

        #
        # Here we should have a raster (out_ds) in the correct Spatial Reference system
        #

        # Get alpha band (either directly or from NODATA value)
        alphaband = out_ds.GetRasterBand(1).GetMaskBand()
        if (alphaband.GetMaskFlags() & gdal.GMF_ALPHA) or out_ds.RasterCount==4 or out_ds.RasterCount==2:
            # TODO: Better test for alpha band in the dataset
            dataBandsCount = out_ds.RasterCount - 1
        else:
            dataBandsCount = out_ds.RasterCount
        
        # Read the georeference 
        out_gt = out_ds.GetGeoTransform()

        # Report error in case rotation/skew is in geotransform (possible only in 'raster' profile)
        if (out_gt[2], out_gt[4]) != (0,0):
            self.error("Georeference of the raster contains rotation or skew. Such raster is not supported. Please use gdalwarp first.")
            # TODO: Do the warping in this case automaticaly

        #
        # Here we expect: pixel is square, no rotation on the raster
        #

        # Output Bounds - coordinates in the output SRS
        ominx = out_gt[0]
        omaxx = out_gt[0]+out_ds.RasterXSize*out_gt[1]
        omaxy = out_gt[3]
        ominy = out_gt[3]-out_ds.RasterYSize*out_gt[1]
        # Note: maybe round(x, 14) to avoid the gdal_translate behaviour, when 0 becomes -1e-15

        if self.options.verbose:
            print("Bounds (output srs):", round(ominx, 13), ominy, omaxx, omaxy)

        #
        # Calculating ranges for tiles in different zoom levels
        #

        geodetic = GlobalGeodetic() # from globalmaptiles.py

        # Generate table with min max tile coordinates for all zoomlevels
        tminmax = list(range(0,32))
        for tz in range(0, 32):
            tminx, tminy = geodetic.LatLonToTile( ominx, ominy, tz )
            tmaxx, tmaxy = geodetic.LatLonToTile( omaxx, omaxy, tz )
            # crop tiles extending world limits (+-180,+-90)
            tminx, tminy = max(0, tminx), max(0, tminy)
            tmaxx, tmaxy = min(2**(tz+1)-1, tmaxx), min(2**tz-1, tmaxy)
            tminmax[tz] = (tminx, tminy, tmaxx, tmaxy)

        # Get the maximal zoom level (closest possible zoom level up on the resolution of raster)
        tminz = geodetic.ZoomForPixelSize( out_gt[1] * max( out_ds.RasterXSize, out_ds.RasterYSize) / float(self.tilesize) )
        if self.options.verbose:
            print ('Min Zoom: %s' % tminz)

        # Get the maximal zoom level (closest possible zoom level up on the resolution of raster)
        tmaxz = geodetic.ZoomForPixelSize( out_gt[1] )
        if self.options.verbose:
            print ('Max Zoom: %s' % tmaxz)
        
        self.inumpyuts_data[_inumpyut] = [tminz,tmaxz,tminmax,out_gt[1],out_gt[5]]
        
        self.inumpyuts_files_or_vrt.append(inumpyut_or_vrt)

        if self.options.verbose:
            print("Bounds (latlong):", ominx, ominy, omaxx, omaxy)

    def make_vrt(self,resx,resy,i):
        inumpyuts = " ".join(self.inumpyuts_files_or_vrt)
        if self.options.verbose:
            print "Building VRT file cesium_%s.vrt" % s
        try:
            res = subprocess.check_output("gdalbuildvrt -srcnodata 0 -resolution user -tr %s %s cesium_%s.vrt %s" % (abs(resx),abs(resy),i,inumpyuts), shell=True)
        except:
            exit(1)
    
    def make_tiles(self):
        # Generate the vrt files for zoom ranges
        i = 0
        tmp_res = -1
        vrt_file = None
        for tz in range(self.tminz,self.tmaxz+1):
            res = self.zoom_resolutions[tz][0]  # I check only with resx, because resy will be positively correlated
            if res != tmp_res:
		if i>0:
			self.vrts[vrt_file][1] = tz-1
                tmp_res = res
                resx = self.zoom_resolutions[tz][0]
                resy = self.zoom_resolutions[tz][1]
                self.make_vrt(resx,resy,i)
                vrt_file = "cesium_%s.vrt" % i
                self.vrts[vrt_file] = [tz,None]
                i += 1
            if tz == self.tmaxz:
                self.vrts[vrt_file][1] = tz
        
        self.ti_cum = 0
        if self.options.createtileindexshp and self.tilelayer is None:
            driver = ogr.GetDriverByName('Esri Shapefile')
            shptileindexfile = os.path.join(self.output,'tilesindex.shp')
            if os.path.exists(shptileindexfile):
                for f in glob.glob(self.output+'/tilesindex.*'):
                    os.remove(f)
            shptileindex = driver.CreateDataSource(shptileindexfile)
            self.tilelayer = shptileindex.CreateLayer('tiles', None, ogr.wkbLineString)
            self.tilelayer.CreateField(ogr.FieldDefn('id', ogr.OFTInteger))
            self.tilelayer.CreateField(ogr.FieldDefn('zoom', ogr.OFTInteger))
            self.tilelayer.CreateField(ogr.FieldDefn('tile', ogr.OFTString))
            self.tilelayer.CreateField(ogr.FieldDefn('children', ogr.OFTInteger))
            
        # Generate parent tiles
        self.generate_fake_parent_tiles()
        
        # For each vrt (i.e. zoom range) generate the tiles
        self.steps = len(self.vrts)
        self.step = 1
        for vrt in self.vrts.keys():
            self.process_vrt(vrt)
            if not self.options.keepfiles:
                try:
                    os.remove(vrt)
                except:
                    pass
            self.step += 1
        
        self.create_layerjsonfile()
        
        if self.options.createtileindexshp and self.tilelayer is not None:
            shptileindex.Destroy()
            shptileindex = self.tilelayer = feat = geom = None
            
        print """Processing finished. Tiles written to "%s".""" % self.output
    
    def process_vrt(self,vrt):
        self.open_inumpyut(vrt)
        self.generate_tiles(vrt)

        
    def open_inumpyut(self,vrt):
        if vrt:
            self.in_ds = gdal.Open(vrt, gdal.GA_ReadOnly)
        else:
            raise Exception("No vrt file was specified")
            
        if self.options.verbose:
            print("Inumpyut file:", "( %sP x %sL - %s bands)" % (self.in_ds.RasterXSize, self.in_ds.RasterYSize, self.in_ds.RasterCount))

        if not self.in_ds:
            # Note: GDAL prints the ERROR message too
            self.error("It is not possible to open the inumpyut file '%s'." % vrt )
            
        if self.in_ds.RasterCount == 0:
            self.error( "Inumpyut file '%s' has no raster band" % vrt )
        
        self.out_ds = self.in_ds
        
        # Get alpha band (either directly or from NODATA value)
        self.alphaband = self.out_ds.GetRasterBand(1).GetMaskBand()
        if (self.alphaband.GetMaskFlags() & gdal.GMF_ALPHA) or self.out_ds.RasterCount==4 or self.out_ds.RasterCount==2:
            self.dataBandsCount = self.out_ds.RasterCount - 1
        else:
            self.dataBandsCount = self.out_ds.RasterCount

    # -------------------------------------------------------------------------
    def make_child_flags(self,N,S,E,W):
        # Cesium format neighbor tiles flags
        HAS_SW = 0x01
        HAS_SE = 0x02
        HAS_NW = 0x04
        HAS_NE = 0x08
        
        NB_FLAGS = 0x00
        
        if N & W:
            NB_FLAGS = NB_FLAGS | HAS_NW
        if N & E:
            NB_FLAGS = NB_FLAGS | HAS_NE
        if S & W:
            NB_FLAGS = NB_FLAGS | HAS_SW
        if S & E:
            NB_FLAGS = NB_FLAGS | HAS_SE
        
        return NB_FLAGS
        
    def generate_fake_parent_tiles(self):
        tx = None
        for tz in range(self.tminz-1,-1,-1):
            tminx, tminy, tmaxx, tmaxy = self.tminmax[tz]
            tminx_c, tminy_c, tmaxx_c, tmaxy_c = self.tminmax[tz+1]
            for ty in range(tmaxy, tminy-1, -1):
                for tx in range(tminx, tmaxx+1):
                    tminx_cpot = tx * 2
                    tmaxx_cpot = tminx_cpot + 1
                    tminy_cpot = ty * 2
                    tmaxy_cpot = tminy_cpot + 1
                    
                    N = S = E = W = False
                    if tminx_cpot >= tminx_c:
                        W = True
                    if tmaxx_cpot <= tmaxx_c:
                        E = True
                    if tminy_cpot >= tminy_c:
                        S = True              
                    if tmaxy_cpot <= tmaxy_c:
                        N = True
                    
                    NB_FLAGS = self.make_child_flags(N,S,E,W)
                    if self.options.verbose:
                        print "Fake tile %s,%s,%s" % (tz,tx,ty)
                    self.write_fake_tile(tz,tx,ty,NB_FLAGS)
   
        # Write missing zero level tile with no children, tx 0 in case the zero level parent tileX is 1, 1 otherwise
        if tx:
            tx = 1-tx
        if tx is None:
            tx = 0

        self.write_fake_tile(0,tx,0,0x00)
            
    def write_fake_tile(self,tz,tx,ty,NB_FLAGS):
        tilefilename = os.path.join(self.tmpoutput, str(tz), str(tx), "%s.%s" % (ty, self.tileext))
        # Create directories for the tile
        if not os.path.exists(os.path.dirname(tilefilename)):
            os.makedirs(os.path.dirname(tilefilename))
            
        if self.options.createtileindexshp and self.tilelayer is not None:
            self.ti_cum += 1
            tilelayerdefn = self.tilelayer.GetLayerDefn()
            feat = ogr.Feature(tilelayerdefn)
            feat.SetField('id', self.ti_cum)
            feat.SetField('zoom', tz)
            feat.SetField('tile', "%s_%s_%s" % (tz, tx, ty))
            feat.SetField('children', NB_FLAGS)
            b = self.geodetic.TileBounds(tx, ty, tz)
            geom = ogr.CreateGeometryFromWkb(makeline(b[0], b[3], b[2], b[1]).wkb)
            feat.SetGeometry(geom)
            self.tilelayer.CreateFeature(feat)
            feat = geom = None
        
        # convert to integer representation of heightmap accordind to Cesium format and append children flags byte
        tilearrayint = (numpy.zeros(4096,numpy.dtype('int16')) + 1000) * 5
        tilearrayint.tofile(tilefilename)
        child_water_bytes = struct.pack('<BB',NB_FLAGS,0x00)
        with open(tilefilename,'ab') as outfile:
            outfile.write(child_water_bytes)
    
    def generate_tiles(self,vrt):
        """Generation of the Csium tiles from the inumpyut raster"""
        print("Generating Tiles (round %s of %s):" % (self.step,self.steps))

        # Cesium format neighbor tiles flags
        HAS_SW = 0x01
        HAS_SE = 0x02
        HAS_NW = 0x04
        HAS_NE = 0x08
        
        tminz = self.vrts[vrt][0]
        tmaxz = self.vrts[vrt][1]
        
        tcount = 0
        for tz in range(tmaxz, tminz-1, -1):
            tminx, tminy, tmaxx, tmaxy = self.tminmax[tz]
            tcount += (1+abs(tmaxx-tminx)) * (1+abs(tmaxy-tminy))
         
        ti = 0
        for tz in range(tmaxz, tminz-1, -1):
            # do not overwrite any real tile with successive inputs' fake tiles
            tminx, tminy, tmaxx, tmaxy = self.tminmax[tz]
            if tz < self.tmaxz:
                tminx_c, tminy_c, tmaxx_c, tmaxy_c = self.tminmax[tz+1]
            
            if self.options.verbose:
                tcount_zoom = (1+abs(tmaxx-tminx)) * (1+abs(tmaxy-tminy))
                print ("Tminx - Tmax: %s-%s" % (tminx,tmaxx))
                print ("Tminy - Tmaxy: %s-%s" % (tminy,tmaxy))
                print("Tile count for zoom %s: %s" % (tz,tcount_zoom))
            
            for ty in range(tmaxy, tminy-1, -1):
                for tx in range(tminx, tmaxx+1):
                    if self.options.resume and os.path.exists(os.path.join(self.tmpoutput, str(tz), str(tx), "%s.%s" % (ty, self.tileext))):
                        continue
                    # By Default the children flags are set to 0, which means no childern tiles (Cesium format)
                    NB_FLAGS = 0x00
                    # Child flags are calculated for all the zooms except the higher one (which has not children tiles)
                    if tz < self.tmaxz:
                        tminx_cpot = tx * 2
                        tmaxx_cpot = tminx_cpot + 1
                        tminy_cpot = ty * 2
                        tmaxy_cpot = tminy_cpot + 1
                        
                        N = S = E = W = False
                        if tminx_cpot >= tminx_c and tminx_cpot <= tmaxx_c:
                            W = True
                        if tmaxx_cpot >= tminx_c and tmaxx_cpot <= tmaxx_c:
                            E = True
                        if tminy_cpot >= tminy_c and tminy_cpot <= tmaxy_c:
                            S = True              
                        if tmaxy_cpot >= tminy_c and tmaxy_cpot <= tmaxy_c:
                            N = True
                            
                        NB_FLAGS = self.make_child_flags(N,S,E,W)
                    if self.stopped:
                        break
                    ti += 1
                    self.ti_cum += 1
                    tilearray  = self.process_tile(tz,tx,ty,ti,NB_FLAGS)
                    self.write_tile(tilearray,tz,tx,ty,NB_FLAGS)
                    if not self.options.verbose:
                        self.progressbar( ti / float(tcount) )


    def process_tile(self,tz,tx,ty,ti,NB_FLAGS):
        ds = self.out_ds
        tilebands = self.dataBandsCount
        querysize = self.querysize

        b = self.geodetic.TileBounds(tx, ty, tz)
        tilesize_aug = self.tilesize + self.extrapixels
        b_aug = self.geodetic.TileBoundsForTileSize(tx, ty, tz, self.extrapixels)
        if self.options.verbose:
            print "Tile bounds %s,%s,%s,%s" % (b[0], b[1], b[2], b[3])
            print "Tile bounds augomented %s,%s,%s,%s" % (b_aug[0], b_aug[1], b_aug[2], b_aug[3])
            #print "Tile poly: %s" % makepoly(b_aug[0], b_aug[1], b_aug[2], b_aug[3]).wkt

        if self.options.createtileindexshp and self.tilelayer is not None:
            '''
            shptileindex.write({
                'geometry': mapping(makepoly(b[0], b[3], b[2], b[1])),
                'properties': {'id': 123},
            })
            '''
            tilelayerdefn = self.tilelayer.GetLayerDefn()
            feat = ogr.Feature(tilelayerdefn)
            feat.SetField('id', self.ti_cum)
            feat.SetField('zoom', tz)
            feat.SetField('tile', "%s_%s_%s" % (tz, tx, ty))
            feat.SetField('children', NB_FLAGS)
            geom = ogr.CreateGeometryFromWkb(makeline(b[0], b[3], b[2], b[1]).wkb)
            feat.SetGeometry(geom)
            self.tilelayer.CreateFeature(feat)
            feat = geom = None
        rb, wb = self.geo_query( ds, b_aug[0], b_aug[3], b_aug[2], b_aug[1])
        nativesize = wb[0]+wb[2] # Pixel size in the raster covering query geo extent
        if self.options.verbose:
            print("\tNative Extent (querysize",nativesize,"): ", rb, wb)

        # Tile bounds in raster coordinates for ReadRaster query with extrapixels for Cesium tiles overlap
        querysize = self.querysize + ((self.querysize/self.tilesize) * self.extrapixels)
        rb, wb = self.geo_query( ds, b_aug[0], b_aug[3], b_aug[2], b_aug[1], querysize=querysize)

        rx, ry, rxsize, rysize = rb
        wx, wy, wxsize, wysize = wb
        if wxsize == 0:
            wxsize = 1
        if wysize == 0:
            wysize = 1

        if self.options.verbose:
            print("\tReadRaster Extent: ", (rx, ry, rxsize, rysize), (wx, wy, wxsize, wysize))

        # Query is in 'nearest neighbour' but can be bigger in then the tilesize
        # We scale down the query to the tilesize by supplied algorithm.
        # Tile dataset in memory
        dstile = self.mem_drv.Create('', tilesize_aug, tilesize_aug, tilebands, gdal.GDT_Float32)
        data = ds.ReadRaster(rx, ry, rxsize, rysize, wxsize, wysize, band_list=list(range(1,self.dataBandsCount+1)))
        datatype = gdal_array.GDALTypeCodeToNumericTypeCode(ds.GetRasterBand(1).DataType)
        if datatype != numpy.float32:
            data = numpy.frombuffer(data, dtype=datatype).astype(numpy.float32).tostring()
        
        if tilesize_aug == querysize:
            # Use the ReadRaster result directly in tiles ('nearest neighbour' query)
            dstile.WriteRaster(wx, wy, wxsize, wysize, data, band_list=list(range(1,self.dataBandsCount+1)))
        else:
            # Big ReadRaster query in memory scaled to the tilesize - all but 'near' algo
            dsquery = self.mem_drv.Create('', querysize, querysize, tilebands, gdal.GDT_Float32)
            dsquery.WriteRaster(wx, wy, wxsize, wysize, data, band_list=list(range(1,self.dataBandsCount+1)))
            self.scale_query_to_tile(dsquery, dstile)
            del dsquery
        del data
        
        tilearray = numpy.array(dstile.ReadAsArray())
        del dstile
        return tilearray
        #return None
    
    def write_tile(self,tilearray,tz,tx,ty,NB_FLAGS,WATER_MASK=0):
        tilefilename = os.path.join(self.tmpoutput, str(tz), str(tx), "%s.%s" % (ty, self.tileext))
        # Create directories for the tile
        if not os.path.exists(os.path.dirname(tilefilename)):
            os.makedirs(os.path.dirname(tilefilename))
        
        # convert to integer representation of heightmap accordind to Cesium format and append children flags byte       
        tilearray = (tilearray+1000) * 5
        tilearrayint = tilearray.astype(numpy.int16)
        tilearrayint.tofile(tilefilename)
        child_water_bytes = struct.pack('<BB',NB_FLAGS,WATER_MASK)
        with open(tilefilename,'ab') as outfile:
            outfile.write(child_water_bytes)
            
    def create_layerjsonfile(self):
        with open(os.path.join(self.output,'layer.json'),'w') as lj:
            lj.write("""{
              "tilejson": "2.1.0",
              "format": "heightmap-1.0",
              "version": "1.0.0",
              "scheme": "tms",
              "tiles": ["{z}/{x}/{y}.terrain"]
            }""")

    # -----------------------------------------------------------------------
    def geo_query(self, ds, ulx, uly, lrx, lry, querysize = 0):
        """For given dataset and query in cartographic coordinates
        returns parameters for ReadRaster() in raster coordinates and
        x/y shifts (for border tiles). If the querysize is not given, the
        extent is returned in the native resolution of dataset ds."""

        geotran = ds.GetGeoTransform()
        rx= int((ulx - geotran[0]) / geotran[1] + 0.001)
        ry= int((uly - geotran[3]) / geotran[5] + 0.001)
        rxsize= int((lrx - ulx) / geotran[1] + 0.5)
        rysize= int((lry - uly) / geotran[5] + 0.5)

        if not querysize:
            wxsize, wysize = rxsize, rysize
        else:
            wxsize, wysize = querysize, querysize

        # Coordinates should not go out of the bounds of the raster
        wx = 0
        if rx < 0:
            rxshift = abs(rx)
            wx = int( wxsize * (float(rxshift) / rxsize) )
            wxsize = wxsize - wx
            rxsize = rxsize - int( rxsize * (float(rxshift) / rxsize) )
            rx = 0
        if rx+rxsize > ds.RasterXSize:
            wxsize = int( wxsize * (float(ds.RasterXSize - rx) / rxsize) )
            rxsize = ds.RasterXSize - rx

        wy = 0
        if ry < 0:
            ryshift = abs(ry)
            wy = int( wysize * (float(ryshift) / rysize) )
            wysize = wysize - wy
            rysize = rysize - int( rysize * (float(ryshift) / rysize) )
            ry = 0
        if ry+rysize > ds.RasterYSize:
            wysize = int( wysize * (float(ds.RasterYSize - ry) / rysize) )
            rysize = ds.RasterYSize - ry

        return (rx, ry, rxsize, rysize), (wx, wy, wxsize, wysize)

    # -------------------------------------------------------------------------
    def scale_query_to_tile(self, dsquery, dstile):
        """Scales down query dataset to the tile dataset"""

        querysize = dsquery.RasterXSize
        tilesize = dstile.RasterXSize
        tilebands = dstile.RasterCount

        if self.options.resampling == 'average':
            for i in range(1,tilebands+1):
                res = gdal.RegenerateOverview( dsquery.GetRasterBand(i),
                    dstile.GetRasterBand(i), 'average' )
                if res != 0:
                    self.error("RegenerateOverview() failed")
        else:
            # Other algorithms are implemented by gdal.ReprojectImage().
            dsquery.SetGeoTransform( (0.0, tilesize / float(querysize), 0.0, 0.0, 0.0, tilesize / float(querysize)) )
            dstile.SetGeoTransform( (0.0, 1.0, 0.0, 0.0, 0.0, 1.0) )
            
            res = gdal.ReprojectImage(dsquery, dstile, None, None, self.resampling)    
            if res != 0:
                self.error("ReprojectImage() failed on %s, error %d" % (tilefilename, res))

if __name__=='__main__':
    argv = gdal.GeneralCmdLineProcessor( sys.argv )
    if argv:
        gdal2cesium = GDAL2Cesium( argv[1:] )
        gdal2cesium.process()