(self)
Initialization of the input raster, reprojection if necessary
Initialization of the input raster, reprojection if necessary
Initialization of the input raster reprojection if necessary
def open_input(self):
        """Initialization of the input raster, reprojection if necessary"""

        gdal.AllRegister()

        # Initialize necessary GDAL drivers

        self.out_drv = gdal.GetDriverByName(self.tiledriver)
        self.mem_drv = gdal.GetDriverByName('MEM')

        if not self.out_drv:
            raise Exception("The '%s' driver was not found, is it available in this GDAL build?"
                            , self.tiledriver)
        if not self.mem_drv:
            raise Exception("The 'MEM' driver was not found, is it available in this GDAL build?"
                            )

        # Open the input file

        if self.input:
            self.in_ds = gdal.Open(self.input, gdal.GA_ReadOnly)
        else:
            raise Exception('No input file was specified')

        if self.options.verbose:
            print('Input 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 input file '%s'."
                        % self.input)

        # Read metadata from the input file

        if self.in_ds.RasterCount == 0:
            self.error("Input file '%s' has no raster band"
                       % self.input)

        if self.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"""
                       % self.input)

        # Get NODATA value

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

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

        #
        # Here we should have RGBA input dataset opened in self.in_ds
        #

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

        # Spatial Reference System of the input raster

        self.in_srs = None

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

            # elif self.options.profile != 'raster':
            #   self.error("There is no spatial reference system info included in the input file.","You should run gdal2tiles with --s_srs EPSG:XXXX or similar.")

        # Spatial Reference System of tiles

        self.out_srs = osr.SpatialReference()

        if self.options.profile == 'mercator':
            self.out_srs.ImportFromEPSG(900913)
        elif self.options.profile == 'geodetic':
            self.out_srs.ImportFromEPSG(4326)
        else:
            self.out_srs = self.in_srs

        # Are the reference systems the same? Reproject if necessary.

        self.out_ds = None

        if self.options.profile in ('mercator', 'geodetic'):

            if self.in_ds.GetGeoTransform() == (
                0.0,
                1.0,
                0.0,
                0.0,
                0.0,
                1.0,
                ) and self.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"
                           )

            if self.in_srs:

                if self.in_srs.ExportToProj4() \
                    != self.out_srs.ExportToProj4() \
                    or self.in_ds.GetGCPCount() != 0:

                    # Generation of VRT dataset in tile projection, default 'nearest neighbour' warping

                    self.out_ds = gdal.AutoCreateWarpedVRT(self.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')")
                        self.out_ds.GetDriver().CreateCopy('tiles.vrt',
                                self.out_ds)

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

                    # Correction of AutoCreateWarpedVRT for NODATA values

                    if self.in_nodata != []:
                        import tempfile
                        tempfilename = tempfile.mktemp('-gdal2tiles.vrt'
                                )
                        self.out_ds.GetDriver().CreateCopy(tempfilename,
                                self.out_ds)

                        # open as a text file

                        s = open(tempfilename).read()

                        # Add the warping options

                        s = s.replace("""<GDALWarpOptions>""",
                                """<GDALWarpOptions>
      <Option name="INIT_DEST">NO_DATA</Option>
      <Option name="UNIFIED_SRC_NODATA">YES</Option>""")

                        # replace BandMapping tag for NODATA bands....

                        for i in range(len(self.in_nodata)):
                            s = \
                                s.replace("""<BandMapping src="%i" dst="%i"/>"""
                                     % (i + 1, i + 1),
                                    """<BandMapping src="%i" dst="%i">
          <SrcNoDataReal>%i</SrcNoDataReal>
          <SrcNoDataImag>0</SrcNoDataImag>
          <DstNoDataReal>%i</DstNoDataReal>
          <DstNoDataImag>0</DstNoDataImag>
        </BandMapping>"""
                                    % (i + 1, i + 1, self.in_nodata[i],
                                    self.in_nodata[i]))  # Or rewrite to white by: , 255 ))

                        # save the corrected VRT

                        open(tempfilename, 'w').write(s)

                        # open by GDAL as self.out_ds

                        self.out_ds = gdal.Open(tempfilename)  # , gdal.GA_ReadOnly)

                        # delete the temporary file

                        os.unlink(tempfilename)

                        # set NODATA_VALUE metadata

                        self.out_ds.SetMetadataItem('NODATA_VALUES',
                                '%i %i %i' % (self.in_nodata[0],
                                self.in_nodata[1], self.in_nodata[2]))

                        if self.options.verbose:
                            print("Modified warping result saved into 'tiles1.vrt'")
                            open('tiles1.vrt', 'w').write(s)

                    # -----------------------------------
                    # Correction of AutoCreateWarpedVRT for Mono (1 band) and RGB (3 bands) files without NODATA:
                    # equivalent of gdalwarp -dstalpha

                    if self.in_nodata == [] and self.out_ds.RasterCount \
                        in [1, 3]:
                        import tempfile
                        tempfilename = tempfile.mktemp('-gdal2tiles.vrt'
                                )
                        self.out_ds.GetDriver().CreateCopy(tempfilename,
                                self.out_ds)

                        # open as a text file

                        s = open(tempfilename).read()

                        # Add the warping options

                        s = s.replace("""<BlockXSize>""",
                                """<VRTRasterBand dataType="Byte" band="%i" subClass="VRTWarpedRasterBand">
    <ColorInterp>Alpha</ColorInterp>
  </VRTRasterBand>
  <BlockXSize>"""
                                % (self.out_ds.RasterCount + 1))
                        s = s.replace("""</GDALWarpOptions>""",
                                """<DstAlphaBand>%i</DstAlphaBand>
  </GDALWarpOptions>"""
                                 % (self.out_ds.RasterCount + 1))
                        s = s.replace("""</WorkingDataType>""",
                                """</WorkingDataType>
    <Option name="INIT_DEST">0</Option>"""
                                )

                        # save the corrected VRT

                        open(tempfilename, 'w').write(s)

                        # open by GDAL as self.out_ds

                        self.out_ds = gdal.Open(tempfilename)  # , gdal.GA_ReadOnly)

                        # delete the temporary file

                        os.unlink(tempfilename)

                        if self.options.verbose:
                            print("Modified -dstalpha warping result saved into 'tiles1.vrt'")
                            open('tiles1.vrt', 'w').write(s)
                    s = '''
                    '''
            else:

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

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

        if not self.out_ds:
            self.out_ds = self.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)

        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:

            # TODO: Better test for alpha band in the dataset

            self.dataBandsCount = self.out_ds.RasterCount - 1
        else:
            self.dataBandsCount = self.out_ds.RasterCount

        # KML test

        self.isepsg4326 = False
        srs4326 = osr.SpatialReference()
        srs4326.ImportFromEPSG(4326)
        if self.out_srs and srs4326.ExportToProj4() \
            == self.out_srs.ExportToProj4():
            self.kml = True
            self.isepsg4326 = True
            if self.options.verbose:
                print('KML autotest OK!')

        # Read the georeference

        self.out_gt = self.out_ds.GetGeoTransform()

        # originX, originY = self.out_gt[0], self.out_gt[3]
        # pixelSize = self.out_gt[1] # = self.out_gt[5]

        # Test the size of the pixel

        # MAPTILER - COMMENTED
        # if self.out_gt[1] != (-1 * self.out_gt[5]) and self.options.profile != 'raster':
            # TODO: Process corectly coordinates with are have swichted Y axis (display in OpenLayers too)
            # self.error("Size of the pixel in the output differ for X and Y axes.")

        # Report error in case rotation/skew is in geotransform (possible only in 'raster' profile)

        if (self.out_gt[2], self.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

        self.ominx = self.out_gt[0]
        self.omaxx = self.out_gt[0] + self.out_ds.RasterXSize \
            * self.out_gt[1]
        self.omaxy = self.out_gt[3]
        self.ominy = self.out_gt[3] - self.out_ds.RasterYSize \
            * self.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(self.ominx, 13),
                   self.ominy, self.omaxx, self.omaxy)

        #
        # Calculating ranges for tiles in different zoom levels
        #

        if self.options.profile == 'mercator':

            self.mercator = GlobalMercator()  # from globalmaptiles.py

            # Function which generates SWNE in LatLong for given tile

            self.tileswne = self.mercator.TileLatLonBounds

            # Generate table with min max tile coordinates for all zoomlevels

            self.tminmax = list(range(0, 32))
            for tz in range(0, 32):
                (tminx, tminy) = self.mercator.MetersToTile(self.ominx,
                        self.ominy, tz)
                (tmaxx, tmaxy) = self.mercator.MetersToTile(self.omaxx,
                        self.omaxy, tz)

                # crop tiles extending world limits (+-180,+-90)

                (tminx, tminy) = (max(0, tminx), max(0, tminy))
                (tmaxx, tmaxy) = (min(2 ** tz - 1, tmaxx), min(2 ** tz
                                  - 1, tmaxy))
                self.tminmax[tz] = (tminx, tminy, tmaxx, tmaxy)

            # TODO: Maps crossing 180E (Alaska?)

            # Get the minimal zoom level (map covers area equivalent to one tile)

            if self.tminz == None:
                self.tminz = \
                    self.mercator.ZoomForPixelSize(self.out_gt[1]
                        * max(self.out_ds.RasterXSize,
                        self.out_ds.RasterYSize) / float(self.tilesize))

            # Get the maximal zoom level (closest possible zoom level up on the resolution of raster)

            if self.tmaxz == None:
                self.tmaxz = \
                    self.mercator.ZoomForPixelSize(self.out_gt[1])

            if self.options.verbose:
                print ('Bounds (latlong):',
                       self.mercator.MetersToLatLon(self.ominx,
                       self.ominy),
                       self.mercator.MetersToLatLon(self.omaxx,
                       self.omaxy))
                print ('MinZoomLevel:', self.tminz)
                print ('MaxZoomLevel:', self.tmaxz, '(',
                       self.mercator.Resolution(self.tmaxz), ')')

        if self.options.profile == 'geodetic':

            self.geodetic = GlobalGeodetic(self.options.tmscompatible)  # from globalmaptiles.py

            # Function which generates SWNE in LatLong for given tile

            self.tileswne = self.geodetic.TileLatLonBounds

            # Generate table with min max tile coordinates for all zoomlevels

            self.tminmax = list(range(0, 32))
            for tz in range(0, 32):
                (tminx, tminy) = self.geodetic.LonLatToTile(self.ominx,
                        self.ominy, tz)
                (tmaxx, tmaxy) = self.geodetic.LonLatToTile(self.omaxx,
                        self.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))
                self.tminmax[tz] = (tminx, tminy, tmaxx, tmaxy)

            # TODO: Maps crossing 180E (Alaska?)

            # Get the maximal zoom level (closest possible zoom level up on the resolution of raster)

            if self.tminz == None:
                self.tminz = \
                    self.geodetic.ZoomForPixelSize(self.out_gt[1]
                        * max(self.out_ds.RasterXSize,
                        self.out_ds.RasterYSize) / float(self.tilesize))

            # Get the maximal zoom level (closest possible zoom level up on the resolution of raster)

            if self.tmaxz == None:
                self.tmaxz = \
                    self.geodetic.ZoomForPixelSize(self.out_gt[1])

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

        if self.options.profile == 'raster':

            log2 = lambda x: math.log10(x) / math.log10(2)  # log2 (base 2 logarithm)

            self.nativezoom = \
                int(max(math.ceil(log2(self.out_ds.RasterXSize
                    / float(self.tilesize))),
                    math.ceil(log2(self.out_ds.RasterYSize
                    / float(self.tilesize)))))

            if int(self.tmaxz or 0) < self.nativezoom:
                self.tmaxz = self.nativezoom

            if self.options.verbose:
                print ('Native zoom of the raster:', self.nativezoom)

            # Get the minimal zoom level (whole raster in one tile)

            if self.tminz == None:
                self.tminz = 0

            # Get the maximal zoom level (native resolution of the raster)

            if self.tmaxz == None:
                self.tmaxz = self.nativezoom

            # Generate table with min max tile coordinates for all zoomlevels

            self.tminmax = list(range(0, self.tmaxz + 1))
            self.tsize = list(range(0, self.tmaxz + 1))
            for tz in range(0, self.tmaxz + 1):
                tsize = 2.0 ** (self.nativezoom - tz) * self.tilesize
                (tminx, tminy) = (0, 0)
                tmaxx = int(math.ceil(self.out_ds.RasterXSize / tsize)) \
                    - 1
                tmaxy = int(math.ceil(self.out_ds.RasterYSize / tsize)) \
                    - 1
                self.tsize[tz] = math.ceil(tsize)
                self.tminmax[tz] = (tminx, tminy, tmaxx, tmaxy)

            # Function which generates SWNE in LatLong for given tile

            if self.kml and self.in_srs_wkt:
                self.ct = osr.CoordinateTransformation(self.in_srs,
                        srs4326)

                def rastertileswne(x, y, z):
                    pixelsizex = 2 ** (self.tmaxz - z) * self.out_gt[1]  # X-pixel size in level
                    pixelsizey = 2 ** (self.tmaxz - z) * self.out_gt[1]  # Y-pixel size in level (usually -1*pixelsizex)
                    west = self.out_gt[0] + x * self.tilesize \
                        * pixelsizex
                    east = west + self.tilesize * pixelsizex
                    south = self.ominy + y * self.tilesize * pixelsizex
                    north = south + self.tilesize * pixelsizex
                    if not self.isepsg4326:

                        # Transformation to EPSG:4326 (WGS84 datum)

                        (west, south) = self.ct.TransformPoint(west,
                                south)[:2]
                        (east, north) = self.ct.TransformPoint(east,
                                north)[:2]
                    return (south, west, north, east)

                self.tileswne = rastertileswne
            else:
                self.tileswne = lambda x, y, z: (0, 0, 0, 0)
def open_input ( self ) : gdal . AllRegister ( ) # Initialize necessary GDAL drivers self . out_drv = gdal . GetDriverByName ( self . tiledriver ) self . mem_drv = gdal . GetDriverByName ( 'MEM' ) if not self . out_drv : raise Exception ( "The '%s' driver was not found, is it available in this GDAL build?" , self . tiledriver ) if not self . mem_drv : raise Exception ( "The 'MEM' driver was not found, is it available in this GDAL build?" ) # Open the input file if self . input : self . in_ds = gdal . Open ( self . input , gdal . GA_ReadOnly ) else : raise Exception ( 'No input file was specified' ) if self . options . verbose : print ( 'Input 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 input file '%s'." % self . input ) # Read metadata from the input file if self . in_ds . RasterCount == 0 : self . error ( "Input file '%s' has no raster band" % self . input ) if self . 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""" % self . input ) # Get NODATA value self . in_nodata = [ ] for i in range ( 1 , self . in_ds . RasterCount + 1 ) : if self . in_ds . GetRasterBand ( i ) . GetNoDataValue ( ) != None : self . in_nodata . append ( self . in_ds . GetRasterBand ( i ) . GetNoDataValue ( ) ) if self . options . srcnodata : nds = list ( map ( float , self . options . srcnodata . split ( ',' ) ) ) if len ( nds ) < self . in_ds . RasterCount : self . in_nodata = ( nds * self . in_ds . RasterCount ) [ : self . in_ds . RasterCount ] else : self . in_nodata = nds if self . options . verbose : print ( 'NODATA: %s' % self . in_nodata ) # # Here we should have RGBA input dataset opened in self.in_ds # if self . options . verbose : print ( 'Preprocessed file:' , '( %sP x %sL - %s bands)' % ( self . in_ds . RasterXSize , self . in_ds . RasterYSize , self . in_ds . RasterCount ) ) # Spatial Reference System of the input raster self . in_srs = None if self . options . s_srs : self . in_srs = osr . SpatialReference ( ) self . in_srs . SetFromUserInput ( self . options . s_srs ) self . in_srs_wkt = self . in_srs . ExportToWkt ( ) else : self . in_srs_wkt = self . in_ds . GetProjection ( ) if not self . in_srs_wkt and self . in_ds . GetGCPCount ( ) != 0 : self . in_srs_wkt = self . in_ds . GetGCPProjection ( ) if self . in_srs_wkt : self . in_srs = osr . SpatialReference ( ) self . in_srs . ImportFromWkt ( self . in_srs_wkt ) # elif self.options.profile != 'raster': #   self.error("There is no spatial reference system info included in the input file.","You should run gdal2tiles with --s_srs EPSG:XXXX or similar.") # Spatial Reference System of tiles self . out_srs = osr . SpatialReference ( ) if self . options . profile == 'mercator' : self . out_srs . ImportFromEPSG ( 900913 ) elif self . options . profile == 'geodetic' : self . out_srs . ImportFromEPSG ( 4326 ) else : self . out_srs = self . in_srs # Are the reference systems the same? Reproject if necessary. self . out_ds = None if self . options . profile in ( 'mercator' , 'geodetic' ) : if self . in_ds . GetGeoTransform ( ) == ( 0.0 , 1.0 , 0.0 , 0.0 , 0.0 , 1.0 , ) and self . 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" ) if self . in_srs : if self . in_srs . ExportToProj4 ( ) != self . out_srs . ExportToProj4 ( ) or self . in_ds . GetGCPCount ( ) != 0 : # Generation of VRT dataset in tile projection, default 'nearest neighbour' warping self . out_ds = gdal . AutoCreateWarpedVRT ( self . 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')" ) self . out_ds . GetDriver ( ) . CreateCopy ( 'tiles.vrt' , self . out_ds ) # Note: self.in_srs and self.in_srs_wkt contain still the non-warped reference system!!! # Correction of AutoCreateWarpedVRT for NODATA values if self . in_nodata != [ ] : import tempfile tempfilename = tempfile . mktemp ( '-gdal2tiles.vrt' ) self . out_ds . GetDriver ( ) . CreateCopy ( tempfilename , self . out_ds ) # open as a text file s = open ( tempfilename ) . read ( ) # Add the warping options s = s . replace ( """<GDALWarpOptions>""" , """<GDALWarpOptions>
      <Option name="INIT_DEST">NO_DATA</Option>
      <Option name="UNIFIED_SRC_NODATA">YES</Option>""" ) # replace BandMapping tag for NODATA bands.... for i in range ( len ( self . in_nodata ) ) : s = s . replace ( """<BandMapping src="%i" dst="%i"/>""" % ( i + 1 , i + 1 ) , """<BandMapping src="%i" dst="%i">
          <SrcNoDataReal>%i</SrcNoDataReal>
          <SrcNoDataImag>0</SrcNoDataImag>
          <DstNoDataReal>%i</DstNoDataReal>
          <DstNoDataImag>0</DstNoDataImag>
        </BandMapping>""" % ( i + 1 , i + 1 , self . in_nodata [ i ] , self . in_nodata [ i ] ) ) # Or rewrite to white by: , 255 )) # save the corrected VRT open ( tempfilename , 'w' ) . write ( s ) # open by GDAL as self.out_ds self . out_ds = gdal . Open ( tempfilename ) # , gdal.GA_ReadOnly) # delete the temporary file os . unlink ( tempfilename ) # set NODATA_VALUE metadata self . out_ds . SetMetadataItem ( 'NODATA_VALUES' , '%i %i %i' % ( self . in_nodata [ 0 ] , self . in_nodata [ 1 ] , self . in_nodata [ 2 ] ) ) if self . options . verbose : print ( "Modified warping result saved into 'tiles1.vrt'" ) open ( 'tiles1.vrt' , 'w' ) . write ( s ) # ----------------------------------- # Correction of AutoCreateWarpedVRT for Mono (1 band) and RGB (3 bands) files without NODATA: # equivalent of gdalwarp -dstalpha if self . in_nodata == [ ] and self . out_ds . RasterCount in [ 1 , 3 ] : import tempfile tempfilename = tempfile . mktemp ( '-gdal2tiles.vrt' ) self . out_ds . GetDriver ( ) . CreateCopy ( tempfilename , self . out_ds ) # open as a text file s = open ( tempfilename ) . read ( ) # Add the warping options s = s . replace ( """<BlockXSize>""" , """<VRTRasterBand dataType="Byte" band="%i" subClass="VRTWarpedRasterBand">
    <ColorInterp>Alpha</ColorInterp>
  </VRTRasterBand>
  <BlockXSize>""" % ( self . out_ds . RasterCount + 1 ) ) s = s . replace ( """</GDALWarpOptions>""" , """<DstAlphaBand>%i</DstAlphaBand>
  </GDALWarpOptions>""" % ( self . out_ds . RasterCount + 1 ) ) s = s . replace ( """</WorkingDataType>""" , """</WorkingDataType>
    <Option name="INIT_DEST">0</Option>""" ) # save the corrected VRT open ( tempfilename , 'w' ) . write ( s ) # open by GDAL as self.out_ds self . out_ds = gdal . Open ( tempfilename ) # , gdal.GA_ReadOnly) # delete the temporary file os . unlink ( tempfilename ) if self . options . verbose : print ( "Modified -dstalpha warping result saved into 'tiles1.vrt'" ) open ( 'tiles1.vrt' , 'w' ) . write ( s ) s = '''
                    ''' else : self . error ( 'Input file has unknown SRS.' , 'Use --s_srs ESPG:xyz (or similar) to provide source reference system.' ) if self . out_ds and self . options . verbose : print ( 'Projected file:' , 'tiles.vrt' , '( %sP x %sL - %s bands)' % ( self . out_ds . RasterXSize , self . out_ds . RasterYSize , self . out_ds . RasterCount ) ) if not self . out_ds : self . out_ds = self . 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) 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 : # TODO: Better test for alpha band in the dataset self . dataBandsCount = self . out_ds . RasterCount - 1 else : self . dataBandsCount = self . out_ds . RasterCount # KML test self . isepsg4326 = False srs4326 = osr . SpatialReference ( ) srs4326 . ImportFromEPSG ( 4326 ) if self . out_srs and srs4326 . ExportToProj4 ( ) == self . out_srs . ExportToProj4 ( ) : self . kml = True self . isepsg4326 = True if self . options . verbose : print ( 'KML autotest OK!' ) # Read the georeference self . out_gt = self . out_ds . GetGeoTransform ( ) # originX, originY = self.out_gt[0], self.out_gt[3] # pixelSize = self.out_gt[1] # = self.out_gt[5] # Test the size of the pixel # MAPTILER - COMMENTED # if self.out_gt[1] != (-1 * self.out_gt[5]) and self.options.profile != 'raster': # TODO: Process corectly coordinates with are have swichted Y axis (display in OpenLayers too) # self.error("Size of the pixel in the output differ for X and Y axes.") # Report error in case rotation/skew is in geotransform (possible only in 'raster' profile) if ( self . out_gt [ 2 ] , self . 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 self . ominx = self . out_gt [ 0 ] self . omaxx = self . out_gt [ 0 ] + self . out_ds . RasterXSize * self . out_gt [ 1 ] self . omaxy = self . out_gt [ 3 ] self . ominy = self . out_gt [ 3 ] - self . out_ds . RasterYSize * self . 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 ( self . ominx , 13 ) , self . ominy , self . omaxx , self . omaxy ) # # Calculating ranges for tiles in different zoom levels # if self . options . profile == 'mercator' : self . mercator = GlobalMercator ( ) # from globalmaptiles.py # Function which generates SWNE in LatLong for given tile self . tileswne = self . mercator . TileLatLonBounds # Generate table with min max tile coordinates for all zoomlevels self . tminmax = list ( range ( 0 , 32 ) ) for tz in range ( 0 , 32 ) : ( tminx , tminy ) = self . mercator . MetersToTile ( self . ominx , self . ominy , tz ) ( tmaxx , tmaxy ) = self . mercator . MetersToTile ( self . omaxx , self . omaxy , tz ) # crop tiles extending world limits (+-180,+-90) ( tminx , tminy ) = ( max ( 0 , tminx ) , max ( 0 , tminy ) ) ( tmaxx , tmaxy ) = ( min ( 2 ** tz - 1 , tmaxx ) , min ( 2 ** tz - 1 , tmaxy ) ) self . tminmax [ tz ] = ( tminx , tminy , tmaxx , tmaxy ) # TODO: Maps crossing 180E (Alaska?) # Get the minimal zoom level (map covers area equivalent to one tile) if self . tminz == None : self . tminz = self . mercator . ZoomForPixelSize ( self . out_gt [ 1 ] * max ( self . out_ds . RasterXSize , self . out_ds . RasterYSize ) / float ( self . tilesize ) ) # Get the maximal zoom level (closest possible zoom level up on the resolution of raster) if self . tmaxz == None : self . tmaxz = self . mercator . ZoomForPixelSize ( self . out_gt [ 1 ] ) if self . options . verbose : print ( 'Bounds (latlong):' , self . mercator . MetersToLatLon ( self . ominx , self . ominy ) , self . mercator . MetersToLatLon ( self . omaxx , self . omaxy ) ) print ( 'MinZoomLevel:' , self . tminz ) print ( 'MaxZoomLevel:' , self . tmaxz , '(' , self . mercator . Resolution ( self . tmaxz ) , ')' ) if self . options . profile == 'geodetic' : self . geodetic = GlobalGeodetic ( self . options . tmscompatible ) # from globalmaptiles.py # Function which generates SWNE in LatLong for given tile self . tileswne = self . geodetic . TileLatLonBounds # Generate table with min max tile coordinates for all zoomlevels self . tminmax = list ( range ( 0 , 32 ) ) for tz in range ( 0 , 32 ) : ( tminx , tminy ) = self . geodetic . LonLatToTile ( self . ominx , self . ominy , tz ) ( tmaxx , tmaxy ) = self . geodetic . LonLatToTile ( self . omaxx , self . 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 ) ) self . tminmax [ tz ] = ( tminx , tminy , tmaxx , tmaxy ) # TODO: Maps crossing 180E (Alaska?) # Get the maximal zoom level (closest possible zoom level up on the resolution of raster) if self . tminz == None : self . tminz = self . geodetic . ZoomForPixelSize ( self . out_gt [ 1 ] * max ( self . out_ds . RasterXSize , self . out_ds . RasterYSize ) / float ( self . tilesize ) ) # Get the maximal zoom level (closest possible zoom level up on the resolution of raster) if self . tmaxz == None : self . tmaxz = self . geodetic . ZoomForPixelSize ( self . out_gt [ 1 ] ) if self . options . verbose : print ( 'Bounds (latlong):' , self . ominx , self . ominy , self . omaxx , self . omaxy ) if self . options . profile == 'raster' : log2 = lambda x : math . log10 ( x ) / math . log10 ( 2 ) # log2 (base 2 logarithm) self . nativezoom = int ( max ( math . ceil ( log2 ( self . out_ds . RasterXSize / float ( self . tilesize ) ) ) , math . ceil ( log2 ( self . out_ds . RasterYSize / float ( self . tilesize ) ) ) ) ) if int ( self . tmaxz or 0 ) < self . nativezoom : self . tmaxz = self . nativezoom if self . options . verbose : print ( 'Native zoom of the raster:' , self . nativezoom ) # Get the minimal zoom level (whole raster in one tile) if self . tminz == None : self . tminz = 0 # Get the maximal zoom level (native resolution of the raster) if self . tmaxz == None : self . tmaxz = self . nativezoom # Generate table with min max tile coordinates for all zoomlevels self . tminmax = list ( range ( 0 , self . tmaxz + 1 ) ) self . tsize = list ( range ( 0 , self . tmaxz + 1 ) ) for tz in range ( 0 , self . tmaxz + 1 ) : tsize = 2.0 ** ( self . nativezoom - tz ) * self . tilesize ( tminx , tminy ) = ( 0 , 0 ) tmaxx = int ( math . ceil ( self . out_ds . RasterXSize / tsize ) ) - 1 tmaxy = int ( math . ceil ( self . out_ds . RasterYSize / tsize ) ) - 1 self . tsize [ tz ] = math . ceil ( tsize ) self . tminmax [ tz ] = ( tminx , tminy , tmaxx , tmaxy ) # Function which generates SWNE in LatLong for given tile if self . kml and self . in_srs_wkt : self . ct = osr . CoordinateTransformation ( self . in_srs , srs4326 ) def rastertileswne ( x , y , z ) : pixelsizex = 2 ** ( self . tmaxz - z ) * self . out_gt [ 1 ] # X-pixel size in level pixelsizey = 2 ** ( self . tmaxz - z ) * self . out_gt [ 1 ] # Y-pixel size in level (usually -1*pixelsizex) west = self . out_gt [ 0 ] + x * self . tilesize * pixelsizex east = west + self . tilesize * pixelsizex south = self . ominy + y * self . tilesize * pixelsizex north = south + self . tilesize * pixelsizex if not self . isepsg4326 : # Transformation to EPSG:4326 (WGS84 datum) ( west , south ) = self . ct . TransformPoint ( west , south ) [ : 2 ] ( east , north ) = self . ct . TransformPoint ( east , north ) [ : 2 ] return ( south , west , north , east ) self . tileswne = rastertileswne else : self . tileswne = lambda x , y , z : ( 0 , 0 , 0 , 0 )