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

Revision 23, 27.8 KB checked in by astrid.emde, 12 years ago (diff)
Line 
1#!/usr/bin/env python
2###############################################################################
3#
4# Purpose:  Module for retiling (merging) tiles and building tiled pyramids
5# Author:   Christian Meuller, christian.mueller@nvoe.at
6#
7###############################################################################
8# Copyright (c) 2007, Christian Mueller
9#
10# Permission is hereby granted, free of charge, to any person obtaining a
11# copy of this software and associated documentation files (the "Software"),
12# to deal in the Software without restriction, including without limitation
13# the rights to use, copy, modify, merge, publish, distribute, sublicense,
14# and/or sell copies of the Software, and to permit persons to whom the
15# Software is furnished to do so, subject to the following conditions:
16#
17# The above copyright notice and this permission notice shall be included
18# in all copies or substantial portions of the Software.
19#
20# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
21# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
22# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
23# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
24# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
25# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
26# DEALINGS IN THE SOFTWARE.
27###############################################################################
28
29
30try:
31    from osgeo import gdal
32    from osgeo import ogr
33    from osgeo import osr
34    from osgeo.gdalconst import *
35except:
36    import gdal
37    import ogr
38    import osr
39    from gdalconst import *
40
41import sys
42import os
43import math
44
45class AffineTransformDecorator:
46    """ A class providing some usefull methods for affine Transformations """
47    def __init__(self, transform ):
48        self.geotransform=transform
49        self.scaleX=self.geotransform[1]
50        self.scaleY=self.geotransform[5]
51        if self.scaleY > 0:
52            self.scaleY*=-1
53        self.ulx = self.geotransform[0]
54        self.uly = self.geotransform[3]
55
56    def pointsFor(self,width,height):
57        xlist=[]
58        ylist=[]
59        w=self.scaleX * width;
60        h=self.scaleY * height;
61
62        xlist.append(self.ulx)
63        ylist.append(self.uly)
64        xlist.append(self.ulx+w)
65        ylist.append(self.uly)
66        xlist.append(self.ulx+w)
67        ylist.append(self.uly+h)
68        xlist.append(self.ulx)
69        ylist.append(self.uly+h)
70        return [ xlist, ylist]
71
72
73class DataSetCache:
74    """ A class for caching source tiles """
75    def __init__(self ):
76        self.cacheSize=8
77        self.queue=[]
78        self.dict={}
79    def get(self,name ):
80
81        if name in self.dict:
82            return self.dict[name]
83        result = gdal.Open(name)
84        if result is None:
85            print("Error openenig:",NameError)
86            sys.exit(1)
87        if len(self.queue)==self.cacheSize:
88            toRemove = self.queue.pop(0)
89            del self.dict[toRemove]
90        self.queue.append(name)
91        self.dict[name]=result
92        return result
93    def __del__(self):
94        for name, dataset in self.dict.items():
95            del dataset
96        del self.queue
97        del self.dict
98
99
100
101class tile_info:
102    """ A class holding info how to tile """
103    def __init__(self,xsize,ysize,tileWidth,tileHeight):
104        self.tileWidth=tileWidth
105        self.tileHeight=tileHeight
106        self.countTilesX= xsize / tileWidth
107        self.countTilesY= ysize / tileHeight
108        self.lastTileWidth = xsize - self.countTilesX *  tileWidth
109        self.lastTileHeight = ysize - self.countTilesY *  tileHeight
110
111        if (self.lastTileWidth > 0 ):
112            self.countTilesX=self.countTilesX+1
113        else:
114            self.lastTileWidth=tileWidth
115
116        if (self.lastTileHeight > 0 ):
117            self.countTilesY=self.countTilesY+1
118        else:
119            self.lastTileHeight=tileHeight
120
121
122
123    def report( self ):
124        print('tileWidth       %d' % self.tileWidth)
125        print('tileHeight      %d' % self.tileHeight)
126        print('countTilesX:    %d' % self.countTilesX)
127        print('countTilesY:    %d' % self.countTilesY)
128        print('lastTileWidth:  %d' % self.lastTileWidth)
129        print('lastTileHeight: %d' % self.lastTileHeight)
130
131
132
133class mosaic_info:
134    """A class holding information about a GDAL file or a GDAL fileset"""
135
136    def __init__(self, filename,inputDS ):
137        """
138        Initialize mosaic_info from filename
139
140        filename -- Name of file to read.
141
142        """
143        self.TempDriver=gdal.GetDriverByName("MEM")
144        self.filename = filename
145        self.cache = DataSetCache()
146        self.ogrTileIndexDS = inputDS
147
148        self.ogrTileIndexDS.GetLayer().ResetReading()
149        feature = self.ogrTileIndexDS.GetLayer().GetNextFeature()
150        imgLocation = feature.GetField(0)
151
152        fhInputTile = self.cache.get(imgLocation)
153
154        self.bands = fhInputTile.RasterCount
155        self.band_type = fhInputTile.GetRasterBand(1).DataType
156        self.projection = fhInputTile.GetProjection()
157
158
159        dec = AffineTransformDecorator(fhInputTile.GetGeoTransform())
160        self.scaleX=dec.scaleX
161        self.scaleY=dec.scaleY
162        ct = fhInputTile.GetRasterBand(1).GetRasterColorTable()
163        if ct is not None:
164           self.ct = ct.Clone()
165        else:
166           self.ct = None
167
168        extent = self.ogrTileIndexDS.GetLayer().GetExtent()
169        self.ulx = extent[0];
170        self.uly = extent[3]
171        self.lrx = extent[1]
172        self.lry = extent[2]
173
174        self.xsize = int(round((self.lrx-self.ulx) / self.scaleX))
175        self.ysize = abs(int(round((self.uly-self.lry) / self.scaleY)))
176
177
178    def __del__(self):
179        del self.cache
180        del self.ogrTileIndexDS
181
182    def getDataSet(self,minx,miny,maxx,maxy):
183
184        self.ogrTileIndexDS.GetLayer().ResetReading()
185        self.ogrTileIndexDS.GetLayer().SetSpatialFilterRect(minx,miny,maxx,maxy)
186        features = []
187        envelope = None
188        while True:
189            feature = self.ogrTileIndexDS.GetLayer().GetNextFeature();
190            if feature is None:
191                break
192            features.append(feature)
193            if envelope is None:
194                envelope=feature.GetGeometryRef().GetEnvelope()
195            else:
196                featureEnv = feature.GetGeometryRef().GetEnvelope()
197                envelope= ( min(featureEnv[0],envelope[0]),max(featureEnv[1],envelope[1]),
198                            min(featureEnv[2],envelope[2]),max(featureEnv[3],envelope[3]))
199
200        if envelope is None:
201            return None
202
203        #enlarge to query rect if necessairy
204        envelope= ( min(minx,envelope[0]),max(maxx,envelope[1]),
205                    min(miny,envelope[2]),max(maxy,envelope[3]))
206
207
208        self.ogrTileIndexDS.GetLayer().SetSpatialFilter(None)
209
210         # merge tiles
211
212
213        resultSizeX =int(math.ceil(((maxx-minx) / self.scaleX )))
214        resultSizeY =int(math.ceil(((miny-maxy) / self.scaleY )))
215
216        resultDS = self.TempDriver.Create( "TEMP", resultSizeX, resultSizeY, self.bands,self.band_type,[])
217        resultDS.SetGeoTransform( [minx,self.scaleX,0,maxy,0,self.scaleY] )
218
219
220        for feature in features:
221            featureName =  feature.GetField(0)
222            sourceDS=self.cache.get(featureName)
223            dec = AffineTransformDecorator(sourceDS.GetGeoTransform())
224            #calculate read and write offsets
225            readOffsetX =int(round((minx-dec.ulx) / self.scaleX))
226            readOffsetY =int(round((maxy-dec.uly) / self.scaleY))
227            writeOffsetX=0
228            if readOffsetX<0:
229                writeOffsetX=readOffsetX*-1;
230                readOffsetX=0
231            writeOffsetY=0
232            if readOffsetY<0:
233                writeOffsetY=readOffsetY*-1;
234                readOffsetY=0
235            #calculate read and write dimensions
236            readX=min(resultSizeX,sourceDS.RasterXSize-readOffsetX,resultSizeX-writeOffsetX)
237            if readX<=0:
238                continue
239            readY=min(resultSizeY,sourceDS.RasterYSize-readOffsetY,resultSizeY-writeOffsetY)
240            if readY<=0:
241                continue
242
243#            print "READ",readOffsetX,readOffsetY,readX,readY
244
245            for bandNr in range(1,self.bands+1):
246                s_band = sourceDS.GetRasterBand( bandNr )
247                t_band = resultDS.GetRasterBand( bandNr )
248                if self.ct is not None:
249                    t_band.SetRasterColorTable(self.ct)
250                data = s_band.ReadRaster( readOffsetX,readOffsetY,readX,readY, readX,readY, self.band_type )
251                t_band.WriteRaster(writeOffsetX,writeOffsetY,readX,readY,data )
252
253        return resultDS
254
255    def closeDataSet(self, memDS):
256        del memDS
257        #self.TempDriver.Delete("TEMP")
258
259
260    def report( self ):
261        print('Filename: '+ self.filename)
262        print('File Size: %dx%dx%d' \
263              % (self.xsize, self.ysize, self.bands))
264        print('Pixel Size: %f x %f' \
265              % (self.scaleX,self.scaleY))
266        print('UL:(%f,%f)   LR:(%f,%f)' \
267              % (self.ulx,self.uly,self.lrx,self.lry))
268
269
270def getTileIndexFromFiles( inputTiles, driverTyp):
271
272    if Verbose:
273        print("Building internal Index for %d tile(s) ..." % len(inputTiles), end=' ')
274
275    ogrTileIndexDS = createTileIndex("TileIndex",TileIndexFieldName,None,driverTyp);
276    for inputTile in inputTiles:
277
278        fhInputTile = gdal.Open(inputTile)
279        if fhInputTile is None:
280             return None
281
282        dec = AffineTransformDecorator(fhInputTile.GetGeoTransform())
283        points = dec.pointsFor(fhInputTile.RasterXSize, fhInputTile.RasterYSize)
284
285        addFeature(ogrTileIndexDS,inputTile,points[0],points[1])
286        del fhInputTile
287
288    if Verbose:
289        print("finished")
290    #ogrTileIndexDS.GetLayer().SyncToDisk()
291    return ogrTileIndexDS
292
293
294
295
296def getTargetDir (level = -1):
297    if level==-1:
298        return TargetDir
299    else:
300        return TargetDir+str(level)+os.sep
301
302
303
304
305def tileImage(minfo, ti ):
306    """
307
308    Tile image in mosaicinfo minfo  based on tileinfo ti
309
310    returns list of created tiles
311
312    """
313
314
315    OGRDS=createTileIndex("TileResult_0", TileIndexFieldName, Source_SRS,TileIndexDriverTyp)
316
317
318    yRange = list(range(1,ti.countTilesY+1))
319    xRange = list(range(1,ti.countTilesX+1))
320
321    for yIndex in yRange:
322        for xIndex in xRange:
323            offsetY=(yIndex-1)* ti.tileHeight
324            offsetX=(xIndex-1)* ti.tileWidth
325            if yIndex==ti.countTilesY:
326                height=ti.lastTileHeight
327            else:
328                height=ti.tileHeight
329
330            if xIndex==ti.countTilesX:
331                width=ti.lastTileWidth
332            else:
333                width=ti.tileWidth
334            tilename=getTileName(minfo,ti, xIndex, yIndex)
335            createTile(minfo, offsetX, offsetY, width, height,tilename,OGRDS)
336
337
338    if TileIndexName is not None:
339        shapeName=getTargetDir()+TileIndexName
340        copyTileIndexToDisk(OGRDS,shapeName)
341
342    if CsvFileName is not None:
343        csvName=getTargetDir()+CsvFileName
344        copyTileIndexToCSV(OGRDS,csvName)
345
346
347    return OGRDS
348
349def copyTileIndexToDisk(OGRDS, fileName):
350    SHAPEDS = createTileIndex(fileName, TileIndexFieldName, OGRDS.GetLayer().GetSpatialRef(), "ESRI Shapefile")
351    OGRDS.GetLayer().ResetReading()
352    while True:
353      feature = OGRDS.GetLayer().GetNextFeature()
354      if feature is None:
355          break
356      newFeature = feature.Clone()
357      basename = os.path.basename(feature.GetField(0))
358      newFeature.SetField(0,basename)
359      SHAPEDS.GetLayer().CreateFeature(newFeature)
360    closeTileIndex(SHAPEDS)
361
362def copyTileIndexToCSV(OGRDS, fileName):
363    csvfile = open(fileName, 'w')
364    OGRDS.GetLayer().ResetReading()
365    while True:
366      feature = OGRDS.GetLayer().GetNextFeature()
367      if feature is None:
368          break
369      basename = os.path.basename(feature.GetField(0))
370      csvfile.write(basename);
371      geom = feature.GetGeometryRef()
372      coords = geom.GetEnvelope();
373
374      for i in range(len(coords)):
375          csvfile.write(CsvDelimiter)
376          csvfile.write("%f" % coords[i])
377      csvfile.write("\n");
378
379    csvfile.close()
380
381
382
383def createPyramidTile(levelMosaicInfo, offsetX, offsetY, width, height,tileName,OGRDS):
384
385    sx= levelMosaicInfo.scaleX*2
386    sy= levelMosaicInfo.scaleY*2
387
388    dec = AffineTransformDecorator([levelMosaicInfo.ulx+offsetX*sx,sx,0,
389                                    levelMosaicInfo.uly+offsetY*sy,0,sy])
390
391
392
393    s_fh = levelMosaicInfo.getDataSet(dec.ulx,dec.uly+height*dec.scaleY,
394                         dec.ulx+width*dec.scaleX,dec.uly)
395    if s_fh is None:
396        return
397
398
399    if OGRDS is not None:
400        points = dec.pointsFor(width, height)
401        addFeature(OGRDS, tileName, points[0], points[1])
402
403
404    if BandType is None:
405        bt=levelMosaicInfo.band_type
406    else:
407        bt=BandType
408
409    geotransform = [dec.ulx, dec.scaleX, 0,dec.uly,0,dec.scaleY]
410
411
412    bands = levelMosaicInfo.bands
413
414    if MemDriver is None:
415        t_fh = Driver.Create( tileName, width, height, bands,bt,CreateOptions)
416    else:
417        t_fh = MemDriver.Create( tileName, width, height, bands,bt,CreateOptions)
418
419    if t_fh is None:
420        print('Creation failed, terminating gdal_tile.')
421        sys.exit( 1 )
422
423
424    t_fh.SetGeoTransform( geotransform )
425    t_fh.SetProjection( levelMosaicInfo.projection)
426    if levelMosaicInfo.ct is not None:
427        for band in range(1,bands+1):
428            t_band = t_fh.GetRasterBand( band )
429            t_band.SetRasterColorTable(levelMosaicInfo.ct)
430
431    res = gdal.ReprojectImage(s_fh,t_fh,None,None,ResamplingMethod)
432    if  res!=0:
433        print("Reprojection failed for %s, error %d" % (tileName,res))
434        sys.exit( 1 )
435
436
437    levelMosaicInfo.closeDataSet(s_fh);
438
439    if MemDriver is not None:
440        tt_fh = Driver.CreateCopy( tileName, t_fh, 0, CreateOptions )
441
442    if Verbose:
443        print(tileName + " : " + str(offsetX)+"|"+str(offsetY)+"-->"+str(width)+"-"+str(height))
444
445
446
447
448def createTile( minfo, offsetX,offsetY,width,height, tilename,OGRDS):
449    """
450
451    Create tile
452    return name of created tile
453
454    """
455
456    if BandType is None:
457        bt=minfo.band_type
458    else:
459        bt=BandType
460
461
462    dec = AffineTransformDecorator([minfo.ulx,minfo.scaleX,0,minfo.uly,0,minfo.scaleY])
463
464
465    s_fh = minfo.getDataSet(dec.ulx+offsetX*dec.scaleX,dec.uly+offsetY*dec.scaleY+height*dec.scaleY,
466                         dec.ulx+offsetX*dec.scaleX+width*dec.scaleX,
467                         dec.uly+offsetY*dec.scaleY)
468    if s_fh is None:
469        return;
470
471
472    geotransform = [dec.ulx+offsetX*dec.scaleX, dec.scaleX, 0,
473                    dec.uly+offsetY*dec.scaleY,  0,dec.scaleY]
474
475    if OGRDS is not None:
476        dec2 = AffineTransformDecorator(geotransform)
477        points = dec2.pointsFor(width, height)
478        addFeature(OGRDS, tilename, points[0], points[1])
479
480
481
482    bands = minfo.bands
483
484    if MemDriver is None:
485        t_fh = Driver.Create( tilename, width, height, bands,bt,CreateOptions)
486    else:
487        t_fh = MemDriver.Create( tilename, width, height, bands,bt,CreateOptions)
488
489    if t_fh is None:
490        print('Creation failed, terminating gdal_tile.')
491        sys.exit( 1 )
492
493    t_fh.SetGeoTransform( geotransform )
494    if Source_SRS is not None:
495        t_fh.SetProjection( Source_SRS.ExportToWkt())
496
497    readX=min(s_fh.RasterXSize,width)
498    readY=min(s_fh.RasterYSize,height)
499    for band in range(1,bands+1):
500        s_band = s_fh.GetRasterBand( band )
501        t_band = t_fh.GetRasterBand( band )
502        if minfo.ct is not None:
503            t_band.SetRasterColorTable(minfo.ct)
504
505#        data = s_band.ReadRaster( offsetX,offsetY,width,height,width,height, t_band.DataType )
506        data = s_band.ReadRaster( 0,0,readX,readY,readX,readY,  t_band.DataType )
507        t_band.WriteRaster( 0,0,readX,readY, data,readX,readY, t_band.DataType )
508
509    minfo.closeDataSet(s_fh);
510
511    if MemDriver is not None:
512        tt_fh = Driver.CreateCopy( tilename, t_fh, 0, CreateOptions )
513
514    if Verbose:
515        print(tilename + " : " + str(offsetX)+"|"+str(offsetY)+"-->"+str(width)+"-"+str(height))
516
517
518
519def createTileIndex(dsName,fieldName,srs,driverName):
520
521    OGRDriver = ogr.GetDriverByName(driverName);
522    if OGRDriver is None:
523        print('ESRI Shapefile driver not found')
524        sys.exit( 1 )
525
526    OGRDataSource=OGRDriver.Open(dsName)
527    if OGRDataSource is not None:
528        OGRDataSource.Destroy()
529        OGRDriver.DeleteDataSource(dsName)
530        if Verbose:
531            print('truncating index '+ dsName)
532
533    OGRDataSource=OGRDriver.CreateDataSource(dsName)
534    if OGRDataSource is None:
535        print('Could not open datasource '+dsName)
536        sys.exit( 1 )
537
538    OGRLayer = OGRDataSource.CreateLayer("index", srs, ogr.wkbPolygon)
539    if OGRLayer is None:
540        print('Could not create Layer')
541        sys.exit( 1 )
542
543    OGRFieldDefn = ogr.FieldDefn(fieldName,ogr.OFTString)
544    if OGRFieldDefn is None:
545        print('Could not create FieldDefn for '+fieldName)
546        sys.exit( 1 )
547
548    OGRFieldDefn.SetWidth(256)
549    if OGRLayer.CreateField(OGRFieldDefn) != 0:
550        print('Could not create Field for '+fieldName)
551        sys.exit( 1 )
552
553    return OGRDataSource
554
555def addFeature(OGRDataSource,location,xlist,ylist):
556
557    OGRLayer=OGRDataSource.GetLayer();
558    OGRFeature = ogr.Feature(OGRLayer.GetLayerDefn())
559    if OGRFeature is None:
560        print('Could not create Feature')
561        sys.exit( 1 )
562
563    OGRFeature.SetField(TileIndexFieldName,location);
564    wkt = 'POLYGON ((%f %f,%f %f,%f %f,%f %f,%f %f ))' % (xlist[0],ylist[0],
565            xlist[1],ylist[1],xlist[2],ylist[2],xlist[3],ylist[3],xlist[0],ylist[0])
566    OGRGeometry=ogr.CreateGeometryFromWkt(wkt,OGRLayer.GetSpatialRef())
567    if (OGRGeometry is None):
568        print('Could not create Geometry')
569        sys.exit( 1 )
570
571    OGRFeature.SetGeometryDirectly(OGRGeometry)
572
573    OGRLayer.CreateFeature(OGRFeature)
574    OGRFeature.Destroy()
575
576def closeTileIndex(OGRDataSource):
577    OGRDataSource.Destroy()
578
579
580def buildPyramid(minfo,createdTileIndexDS,tileWidth, tileHeight):
581
582    inputDS=createdTileIndexDS
583    for level in range(1,Levels+1):
584        levelMosaicInfo = mosaic_info(minfo.filename,inputDS)
585        levelOutputTileInfo = tile_info(levelMosaicInfo.xsize/2,levelMosaicInfo.ysize/2,tileWidth,tileHeight)
586        inputDS=buildPyramidLevel(levelMosaicInfo,levelOutputTileInfo,level)
587
588
589def buildPyramidLevel(levelMosaicInfo,levelOutputTileInfo, level):
590    yRange = list(range(1,levelOutputTileInfo.countTilesY+1))
591    xRange = list(range(1,levelOutputTileInfo.countTilesX+1))
592
593    OGRDS=createTileIndex("TileResult_"+str(level), TileIndexFieldName, Source_SRS,TileIndexDriverTyp)
594
595    for yIndex in yRange:
596        for xIndex in xRange:
597            offsetY=(yIndex-1)* levelOutputTileInfo.tileHeight
598            offsetX=(xIndex-1)* levelOutputTileInfo.tileWidth
599            if yIndex==levelOutputTileInfo.countTilesY:
600                height=levelOutputTileInfo.lastTileHeight
601            else:
602                height=levelOutputTileInfo.tileHeight
603
604            if xIndex==levelOutputTileInfo.countTilesX:
605                width=levelOutputTileInfo.lastTileWidth
606            else:
607                width=levelOutputTileInfo.tileWidth
608            tilename=getTileName(levelMosaicInfo,levelOutputTileInfo, xIndex, yIndex,level)
609            createPyramidTile(levelMosaicInfo, offsetX, offsetY, width, height,tilename,OGRDS)
610
611
612    if TileIndexName is not None:
613        shapeName=getTargetDir(level)+TileIndexName
614        copyTileIndexToDisk(OGRDS,shapeName)
615
616    if CsvFileName is not None:
617        csvName=getTargetDir(level)+CsvFileName
618        copyTileIndexToCSV(OGRDS,csvName)
619
620
621    return OGRDS
622
623def getTileName(minfo,ti,xIndex,yIndex,level = -1):
624    """
625    creates the tile file name
626    """
627    max = ti.countTilesX
628    if (ti.countTilesY > max):
629        max=ti.countTilesY
630    countDigits= len(str(max))
631    parts=os.path.splitext(os.path.basename(minfo.filename))
632    if parts[0][0]=="@" : #remove possible leading "@"
633       parts = ( parts[0][1:len(parts[0])], parts[1])
634
635    if Extension is None:
636        format=getTargetDir(level)+parts[0]+"_%0"+str(countDigits)+"i"+"_%0"+str(countDigits)+"i"+parts[1]
637    else:
638        format=getTargetDir(level)+parts[0]+"_%0"+str(countDigits)+"i"+"_%0"+str(countDigits)+"i"+"."+Extension
639    return format % (yIndex,xIndex)
640
641def UsageFormat():
642    print('Valid formats:')
643    count = gdal.GetDriverCount()
644    for index in range(count):
645       driver= gdal.GetDriver(index)
646       print(driver.ShortName)
647
648# =============================================================================
649def Usage():
650     print('Usage: gdal_retile.py ')
651     print('        [-v] [-co NAME=VALUE]* [-of out_format]')
652     print('        [-ps pixelWidth pixelHeight]')
653     print('        [-ot  {Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/')
654     print('               CInt16/CInt32/CFloat32/CFloat64}]')
655     print('        [ -tileIndex tileIndexName [-tileIndexField fieldName]]')
656     print('        [ -csv fileName [-csvDelim delimiter]]')
657     print('        [-s_srs srs_def]  [-pyramidOnly] -levels numberoflevels')
658     print('        [-r {near/bilinear/cubic/cubicspline/lanczos}]')
659     print('        -targetDir TileDirectory input_files')
660
661# =============================================================================
662
663
664# =============================================================================
665#
666# Program mainline.
667#
668
669def main(args):
670
671    global Verbose
672    global CreateOptions
673    global Names
674    global TileWidth
675    global TileHeight
676    global Format
677    global BandType
678    global Driver
679    global Extension
680    global MemDriver
681    global TileIndexFieldName
682    global TileIndexName
683    global CsvDelimiter
684    global CsvFileName
685
686    global TileIndexDriverTyp
687    global Source_SRS
688    global TargetDir
689    global ResamplingMethod
690    global Levels
691    global PyramidOnly
692
693    gdal.AllRegister()
694
695    argv = gdal.GeneralCmdLineProcessor( args )
696    if argv is None:
697        return 1
698
699    # Parse command line arguments.
700    i = 1
701    while i < len(argv):
702        arg = argv[i]
703
704        if arg == '-of':
705            i+=1
706            Format = argv[i]
707        elif arg == '-ot':
708            i+=1
709            BandType = gdal.GetDataTypeByName( argv[i] )
710            if BandType == gdal.GDT_Unknown:
711                print('Unknown GDAL data type: ', argv[i])
712                return 1
713        elif arg == '-co':
714            i+=1
715            CreateOptions.append( argv[i] )
716
717
718        elif arg == '-v':
719            Verbose = True
720
721        elif arg == '-targetDir':
722            i+=1
723            TargetDir=argv[i]
724
725            if os.path.exists(TargetDir)==False:
726                print("TargetDir " + TargetDir + " does not exist")
727                return 1
728            if TargetDir[len(TargetDir)-1:] != os.sep:
729                TargetDir =  TargetDir+os.sep
730
731        elif arg == '-ps':
732            i+=1
733            TileWidth=int(argv[i])
734            i+=1
735            TileHeight=int(argv[i])
736
737        elif arg == '-r':
738            i+=1
739            ResamplingMethodString=argv[i]
740            if ResamplingMethodString=="near":
741                ResamplingMethod=GRA_NearestNeighbour
742            elif ResamplingMethodString=="bilinear":
743                 ResamplingMethod=GRA_Bilinear
744            elif ResamplingMethodString=="cubic":
745                 ResamplingMethod=GRA_Cubic
746            elif ResamplingMethodString=="cubicspline":
747                 ResamplingMethod=GRA_CubicSpline
748            elif ResamplingMethodString=="lanczos":
749                ResamplingMethod=GRA_Lanczos
750            else:
751                print("Unknown resampling method:" ,ResamplingMethodString)
752                return 1
753        elif arg == '-levels':
754            i+=1
755            Levels=int(argv[i])
756            if Levels<1:
757                print("Invalid number of levels", Levels)
758                return 1
759        elif arg == '-s_srs':
760            i+=1
761            Source_SRS = osr.SpatialReference()
762            if Source_SRS.SetFromUserInput( argv[i] ) != 0:
763                print('invalid -s_srs: ' + argv[i]);
764                return 1;
765
766        elif arg ==  "-pyramidOnly":
767            PyramidOnly=True
768        elif arg == '-tileIndex':
769            i+=1
770            TileIndexName=argv[i]
771            parts=os.path.splitext(TileIndexName)
772            if len(parts[1])==0:
773                TileIndexName+=".shp"
774
775        elif arg == '-tileIndexField':
776            i+=1
777            TileIndexFieldName=argv[i]
778        elif arg == '-csv':
779            i+=1
780            CsvFileName=argv[i]
781            parts=os.path.splitext(CsvFileName)
782            if len(parts[1])==0:
783                CsvFileName+=".csv"
784        elif arg == '-csvDelim':
785            i+=1
786            CsvDelimiter=argv[i]
787        elif arg[:1] == '-':
788            print('Unrecognised command option: ', arg)
789            Usage()
790            return 1
791
792        else:
793            Names.append( arg )
794        i+=1
795
796    if len(Names) == 0:
797        print('No input files selected.')
798        Usage()
799        return 1
800
801    if (TileWidth==0 or TileHeight==0):
802        print("Invalid tile dimension %d,%d" % (TileWidth,TileHeight))
803        return 1
804
805    if (TargetDir is None):
806        print("Missing Directory for Tiles -targetDir")
807        Usage()
808        return 1
809
810    if Levels > 0:    #prepare Dirs for pyramid
811        for levelIndx in range (1,Levels+1):
812            leveldir=TargetDir+str(levelIndx)+os.sep
813            if (os.path.exists(leveldir)):
814                continue
815            os.mkdir(leveldir)
816            if (os.path.exists(leveldir)==False):
817                print("Cannot create level dir:", leveldir)
818                return 1
819            if Verbose :
820                print("Created level dir: ",leveldir)
821
822
823    Driver = gdal.GetDriverByName(Format)
824    if Driver is None:
825        print('Format driver %s not found, pick a supported driver.' % Format)
826        UsageFormat()
827        return 1
828
829
830
831
832    DriverMD = Driver.GetMetadata()
833    Extension=DriverMD.get(DMD_EXTENSION);
834    if 'DCAP_CREATE' not in DriverMD:
835        MemDriver=gdal.GetDriverByName("MEM")
836
837
838    tileIndexDS=getTileIndexFromFiles(Names,TileIndexDriverTyp)
839    if tileIndexDS is None:
840        print("Error building tile index")
841        return 1;
842    minfo = mosaic_info(Names[0],tileIndexDS)
843    ti=tile_info(minfo.xsize,minfo.ysize, TileWidth, TileHeight)
844
845    if Source_SRS is None and len(minfo.projection) > 0 :
846       Source_SRS = osr.SpatialReference()
847       if Source_SRS.SetFromUserInput( minfo.projection ) != 0:
848           print('invalid projection  ' + minfo.projection);
849           return 1
850
851    if Verbose:
852        minfo.report()
853        ti.report()
854
855
856    if PyramidOnly==False:
857       dsCreatedTileIndex = tileImage(minfo,ti)
858       tileIndexDS.Destroy()
859    else:
860       dsCreatedTileIndex=tileIndexDS
861
862    if Levels>0:
863       buildPyramid(minfo,dsCreatedTileIndex,TileWidth, TileHeight)
864
865    if Verbose:
866        print("FINISHED")
867    return 0
868
869def initGlobals():
870    """ Only used for unit tests """
871    global Verbose
872    global CreateOptions
873    global Names
874    global TileWidth
875    global TileHeight
876    global Format
877    global BandType
878    global Driver
879    global Extension
880    global MemDriver
881    global TileIndexFieldName
882    global TileIndexName
883    global TileIndexDriverTyp
884    global CsvDelimiter
885    global CsvFileName
886    global Source_SRS
887    global TargetDir
888    global ResamplingMethod
889    global Levels
890    global PyramidOnly
891
892    Verbose=False
893    CreateOptions = []
894    Names=[]
895    TileWidth=256
896    TileHeight=256
897    Format='GTiff'
898    BandType = None
899    Driver=None
900    Extension=None
901    MemDriver=None
902    TileIndexFieldName='location'
903    TileIndexName=None
904    TileIndexDriverTyp="Memory"
905    CsvDelimiter=";"
906    CsvFileName=None
907
908    Source_SRS=None
909    TargetDir=None
910    ResamplingMethod=GRA_NearestNeighbour
911    Levels=0
912    PyramidOnly=False
913
914
915#global vars
916Verbose=False
917CreateOptions = []
918Names=[]
919TileWidth=256
920TileHeight=256
921Format='GTiff'
922BandType = None
923Driver=None
924Extension=None
925MemDriver=None
926TileIndexFieldName='location'
927TileIndexName=None
928TileIndexDriverTyp="Memory"
929CsvDelimiter=";"
930CsvFileName=None
931Source_SRS=None
932TargetDir=None
933ResamplingMethod=GRA_NearestNeighbour
934Levels=0
935PyramidOnly=False
936
937
938
939if __name__ == '__main__':
940    sys.exit(main(sys.argv))
941
Note: See TracBrowser for help on using the repository browser.