1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31 from __future__ import print_function
32 import os
33 import sys
34 import shutil
35
36 from osgeo import gdal
37 from osgeo import ogr
38 from osgeo import osr
39
40 progress = gdal.TermProgress_nocb
41
42
72
73
75 """ A class for caching source tiles """
76
78 self.cacheSize = 8
79 self.queue = []
80 self.dict = {}
81
82 - def get(self, name):
83
84 if name in self.dict:
85 return self.dict[name]
86 result = gdal.Open(name)
87 if result is None:
88 print("Error opening: %s" % NameError)
89 sys.exit(1)
90 if len(self.queue) == self.cacheSize:
91 toRemove = self.queue.pop(0)
92 del self.dict[toRemove]
93 self.queue.append(name)
94 self.dict[name] = result
95 return result
96
98 for dataset in self.dict.values():
99 del dataset
100 del self.queue
101 del self.dict
102
103
105 """ A class holding info how to tile """
106
107 - def __init__(self, xsize, ysize, tileWidth, tileHeight, overlap):
108 self.width = xsize
109 self.height = ysize
110 self.tileWidth = tileWidth
111 self.tileHeight = tileHeight
112 self.countTilesX = 1
113 if xsize > tileWidth:
114 self.countTilesX += int((xsize - tileWidth + (tileWidth - overlap) - 1) / (tileWidth - overlap))
115 self.countTilesY = 1
116 if ysize > tileHeight:
117 self.countTilesY += int((ysize - tileHeight + (tileHeight - overlap) - 1) / (tileHeight - overlap))
118 self.overlap = overlap
119
121 print('tileWidth: %d' % self.tileWidth)
122 print('tileHeight: %d' % self.tileHeight)
123 print('countTilesX: %d' % self.countTilesX)
124 print('countTilesY: %d' % self.countTilesY)
125 print('overlap: %d' % self.overlap)
126
127
129 """A class holding information about a GDAL file or a GDAL fileset"""
130
174
176 del self.cache
177 del self.ogrTileIndexDS
178
180
181 self.ogrTileIndexDS.GetLayer().ResetReading()
182 self.ogrTileIndexDS.GetLayer().SetSpatialFilterRect(minx, miny, maxx, maxy)
183 features = []
184 envelope = None
185 while True:
186 feature = self.ogrTileIndexDS.GetLayer().GetNextFeature()
187 if feature is None:
188 break
189 features.append(feature)
190 if envelope is None:
191 envelope = feature.GetGeometryRef().GetEnvelope()
192 else:
193 featureEnv = feature.GetGeometryRef().GetEnvelope()
194 envelope = (min(featureEnv[0], envelope[0]), max(featureEnv[1], envelope[1]),
195 min(featureEnv[2], envelope[2]), max(featureEnv[3], envelope[3]))
196
197 if envelope is None:
198 return None
199
200
201 envelope = (min(minx, envelope[0]), max(maxx, envelope[1]),
202 min(miny, envelope[2]), max(maxy, envelope[3]))
203
204 self.ogrTileIndexDS.GetLayer().SetSpatialFilter(None)
205
206
207
208 resultSizeX = int((maxx - minx) / self.scaleX + 0.5)
209 resultSizeY = int((miny - maxy) / self.scaleY + 0.5)
210
211 resultDS = self.TempDriver.Create("TEMP", resultSizeX, resultSizeY, self.bands, self.band_type, [])
212 resultDS.SetGeoTransform([minx, self.scaleX, 0, maxy, 0, self.scaleY])
213
214 for bandNr in range(1, self.bands + 1):
215 t_band = resultDS.GetRasterBand(bandNr)
216 if self.nodata is not None:
217 t_band.Fill(self.nodata)
218 t_band.SetNoDataValue(self.nodata)
219
220 for feature in features:
221 featureName = feature.GetField(0)
222 sourceDS = self.cache.get(featureName)
223 dec = AffineTransformDecorator(sourceDS.GetGeoTransform())
224
225 dec.lrx = dec.ulx + sourceDS.RasterXSize * dec.scaleX
226 dec.lry = dec.uly + sourceDS.RasterYSize * dec.scaleY
227
228
229 tgw_ulx = max(dec.ulx, minx)
230 tgw_lrx = min(dec.lrx, maxx)
231 if self.scaleY < 0:
232 tgw_uly = min(dec.uly, maxy)
233 tgw_lry = max(dec.lry, miny)
234 else:
235 tgw_uly = max(dec.uly, maxy)
236 tgw_lry = min(dec.lry, miny)
237
238
239 sw_xoff = int((tgw_ulx - dec.ulx) / dec.scaleX + 0.5)
240 sw_yoff = int((tgw_uly - dec.uly) / dec.scaleY + 0.5)
241 sw_xsize = min(sourceDS.RasterXSize, int((tgw_lrx - dec.ulx) / dec.scaleX + 0.5)) - sw_xoff
242 sw_ysize = min(sourceDS.RasterYSize, int((tgw_lry - dec.uly) / dec.scaleY + 0.5)) - sw_yoff
243 if sw_xsize <= 0 or sw_ysize <= 0:
244 continue
245
246
247 tw_xoff = int((tgw_ulx - minx) / self.scaleX + 0.5)
248 tw_yoff = int((tgw_uly - maxy) / self.scaleY + 0.5)
249 tw_xsize = min(resultDS.RasterXSize, int((tgw_lrx - minx) / self.scaleX + 0.5)) - tw_xoff
250 tw_ysize = min(resultDS.RasterYSize, int((tgw_lry - maxy) / self.scaleY + 0.5)) - tw_yoff
251 if tw_xsize <= 0 or tw_ysize <= 0:
252 continue
253
254 assert tw_xoff >= 0
255 assert tw_yoff >= 0
256 assert sw_xoff >= 0
257 assert sw_yoff >= 0
258
259 for bandNr in range(1, self.bands + 1):
260 s_band = sourceDS.GetRasterBand(bandNr)
261 t_band = resultDS.GetRasterBand(bandNr)
262 if self.ct is not None:
263 t_band.SetRasterColorTable(self.ct)
264 t_band.SetRasterColorInterpretation(self.ci[bandNr - 1])
265
266 data = s_band.ReadRaster(sw_xoff, sw_yoff, sw_xsize, sw_ysize, tw_xsize, tw_ysize, self.band_type)
267 if data is None:
268 print(gdal.GetLastErrorMsg())
269
270 t_band.WriteRaster(tw_xoff, tw_yoff, tw_xsize, tw_ysize, data)
271
272 return resultDS
273
276
277
279 print('Filename: ' + self.filename)
280 print('File Size: %dx%dx%d'
281 % (self.xsize, self.ysize, self.bands))
282 print('Pixel Size: %f x %f'
283 % (self.scaleX, self.scaleY))
284 print('UL:(%f,%f) LR:(%f,%f)'
285 % (self.ulx, self.uly, self.lrx, self.lry))
286
287
289 if g.Verbose:
290 print("Building internal Index for %d tile(s) ..." % len(g.Names), end=" ")
291
292 ogrTileIndexDS = createTileIndex(g.Verbose, "TileIndex", g.TileIndexFieldName, None, g.TileIndexDriverTyp)
293 for inputTile in g.Names:
294
295 fhInputTile = gdal.Open(inputTile)
296 if fhInputTile is None:
297 return None
298
299 dec = AffineTransformDecorator(fhInputTile.GetGeoTransform())
300 points = dec.pointsFor(fhInputTile.RasterXSize, fhInputTile.RasterYSize)
301
302 addFeature(g.TileIndexFieldName, ogrTileIndexDS, inputTile, points[0], points[1])
303 del fhInputTile
304
305 if g.Verbose:
306 print("finished")
307
308 return ogrTileIndexDS
309
310
312 if level == -1:
313 return g.TargetDir
314 return g.TargetDir + str(level) + os.sep
315
316
318 """
319
320 Tile image in mosaicinfo minfo based on tileinfo ti
321
322 returns list of created tiles
323
324 """
325
326 g.LastRowIndx = -1
327 OGRDS = createTileIndex(g.Verbose, "TileResult_0", g.TileIndexFieldName, g.Source_SRS, g.TileIndexDriverTyp)
328
329 yRange = list(range(1, ti.countTilesY + 1))
330 xRange = list(range(1, ti.countTilesX + 1))
331
332 if not g.Quiet and not g.Verbose:
333 progress(0.0)
334 processed = 0
335 total = len(xRange) * len(yRange)
336
337 for yIndex in yRange:
338 for xIndex in xRange:
339 offsetY = (yIndex - 1) * (ti.tileHeight - ti.overlap)
340 offsetX = (xIndex - 1) * (ti.tileWidth - ti.overlap)
341 height = ti.tileHeight
342 width = ti.tileWidth
343 if g.UseDirForEachRow:
344 tilename = getTileName(g, minfo, ti, xIndex, yIndex, 0)
345 else:
346 tilename = getTileName(g, minfo, ti, xIndex, yIndex)
347
348 if offsetX + width > ti.width:
349 width = ti.width - offsetX
350 if offsetY + height > ti.height:
351 height = ti.height - offsetY
352
353 feature_only = g.Resume and os.path.exists(tilename)
354 createTile(g, minfo, offsetX, offsetY, width, height, tilename, OGRDS, feature_only)
355
356 if not g.Quiet and not g.Verbose:
357 processed += 1
358 progress(processed / float(total))
359
360 if g.TileIndexName is not None:
361 if g.UseDirForEachRow and not g.PyramidOnly:
362 shapeName = getTargetDir(g, 0) + g.TileIndexName
363 else:
364 shapeName = getTargetDir(g) + g.TileIndexName
365 copyTileIndexToDisk(g, OGRDS, shapeName)
366
367 if g.CsvFileName is not None:
368 if g.UseDirForEachRow and not g.PyramidOnly:
369 csvName = getTargetDir(g, 0) + g.CsvFileName
370 else:
371 csvName = getTargetDir(g) + g.CsvFileName
372 copyTileIndexToCSV(g, OGRDS, csvName)
373
374 return OGRDS
375
376
392
393
395 csvfile = open(fileName, 'w')
396 OGRDS.GetLayer().ResetReading()
397 while True:
398 feature = OGRDS.GetLayer().GetNextFeature()
399 if feature is None:
400 break
401 basename = os.path.basename(feature.GetField(0))
402 if g.UseDirForEachRow:
403 t = os.path.split(os.path.dirname(feature.GetField(0)))
404 basename = t[1] + "/" + basename
405 csvfile.write(basename)
406 geom = feature.GetGeometryRef()
407 coords = geom.GetEnvelope()
408
409 for coord in coords:
410 csvfile.write(g.CsvDelimiter)
411 csvfile.write("%f" % coord)
412 csvfile.write("\n")
413
414 csvfile.close()
415
416
417 -def createPyramidTile(g, levelMosaicInfo, offsetX, offsetY, width, height, tileName, OGRDS, feature_only):
418
419 temp_tilename = tileName + '.tmp'
420
421 sx = levelMosaicInfo.scaleX * 2
422 sy = levelMosaicInfo.scaleY * 2
423
424 dec = AffineTransformDecorator([levelMosaicInfo.ulx + offsetX * sx, sx, 0,
425 levelMosaicInfo.uly + offsetY * sy, 0, sy])
426
427 if OGRDS is not None:
428 points = dec.pointsFor(width, height)
429 addFeature(g.TileIndexFieldName, OGRDS, tileName, points[0], points[1])
430
431 if feature_only:
432 return
433
434 s_fh = levelMosaicInfo.getDataSet(dec.ulx, dec.uly + height * dec.scaleY,
435 dec.ulx + width * dec.scaleX, dec.uly)
436 if s_fh is None:
437 return
438
439 if g.BandType is None:
440 bt = levelMosaicInfo.band_type
441 else:
442 bt = g.BandType
443
444 geotransform = [dec.ulx, dec.scaleX, 0, dec.uly, 0, dec.scaleY]
445
446 bands = levelMosaicInfo.bands
447
448 if g.MemDriver is None:
449 t_fh = g.Driver.Create(temp_tilename, width, height, bands, bt, g.CreateOptions)
450 else:
451 t_fh = g.MemDriver.Create('', width, height, bands, bt)
452
453 if t_fh is None:
454 print('Creation failed, terminating gdal_tile.')
455 sys.exit(1)
456
457 t_fh.SetGeoTransform(geotransform)
458 t_fh.SetProjection(levelMosaicInfo.projection)
459 for band in range(1, bands + 1):
460 t_band = t_fh.GetRasterBand(band)
461 if levelMosaicInfo.ct is not None:
462 t_band.SetRasterColorTable(levelMosaicInfo.ct)
463 t_band.SetRasterColorInterpretation(levelMosaicInfo.ci[band - 1])
464
465 if levelMosaicInfo.nodata is not None:
466 for band in range(1, bands + 1):
467 t_band = t_fh.GetRasterBand(band)
468 t_band.Fill(levelMosaicInfo.nodata)
469 t_band.SetNoDataValue(levelMosaicInfo.nodata)
470
471 res = gdal.ReprojectImage(s_fh, t_fh, None, None, g.ResamplingMethod)
472 if res != 0:
473 print("Reprojection failed for %s, error %d" % (temp_tilename, res))
474 sys.exit(1)
475
476 levelMosaicInfo.closeDataSet(s_fh)
477
478 if g.MemDriver is None:
479 t_fh.FlushCache()
480 else:
481 tt_fh = g.Driver.CreateCopy(temp_tilename, t_fh, 0, g.CreateOptions)
482 tt_fh.FlushCache()
483 tt_fh = None
484
485 t_fh = None
486
487 if os.path.exists(tileName):
488 os.remove(tileName)
489 shutil.move(temp_tilename, tileName)
490
491 if g.Verbose:
492 print(tileName + " : " + str(offsetX) + "|" + str(offsetY) + "-->" + str(width) + "-" + str(height))
493
494
495 -def createTile(g, minfo, offsetX, offsetY, width, height, tilename, OGRDS, feature_only):
496 """
497
498 Create tile
499
500 """
501 temp_tilename = tilename + '.tmp'
502
503 if g.BandType is None:
504 bt = minfo.band_type
505 else:
506 bt = g.BandType
507
508 dec = AffineTransformDecorator([minfo.ulx, minfo.scaleX, 0, minfo.uly, 0, minfo.scaleY])
509
510 geotransform = [dec.ulx + offsetX * dec.scaleX, dec.scaleX, 0,
511 dec.uly + offsetY * dec.scaleY, 0, dec.scaleY]
512
513 if OGRDS is not None:
514 dec2 = AffineTransformDecorator(geotransform)
515 points = dec2.pointsFor(width, height)
516 addFeature(g.TileIndexFieldName, OGRDS, tilename, points[0], points[1])
517
518 if feature_only:
519 return
520
521 s_fh = minfo.getDataSet(dec.ulx + offsetX * dec.scaleX, dec.uly + offsetY * dec.scaleY + height * dec.scaleY,
522 dec.ulx + offsetX * dec.scaleX + width * dec.scaleX,
523 dec.uly + offsetY * dec.scaleY)
524 if s_fh is None:
525 return
526
527 bands = minfo.bands
528
529 if g.MemDriver is None:
530 t_fh = g.Driver.Create(temp_tilename, width, height, bands, bt, g.CreateOptions)
531 else:
532 t_fh = g.MemDriver.Create('', width, height, bands, bt)
533
534 if t_fh is None:
535 print('Creation failed, terminating gdal_tile.')
536 sys.exit(1)
537
538 t_fh.SetGeoTransform(geotransform)
539 if g.Source_SRS is not None:
540 t_fh.SetProjection(g.Source_SRS.ExportToWkt())
541
542 readX = min(s_fh.RasterXSize, width)
543 readY = min(s_fh.RasterYSize, height)
544 for band in range(1, bands + 1):
545 s_band = s_fh.GetRasterBand(band)
546 t_band = t_fh.GetRasterBand(band)
547 if minfo.ct is not None:
548 t_band.SetRasterColorTable(minfo.ct)
549 if minfo.nodata is not None:
550 t_band.Fill(minfo.nodata)
551 t_band.SetNoDataValue(minfo.nodata)
552
553 data = s_band.ReadRaster(0, 0, readX, readY, readX, readY, t_band.DataType)
554 t_band.WriteRaster(0, 0, readX, readY, data, readX, readY, t_band.DataType)
555
556 minfo.closeDataSet(s_fh)
557
558 if g.MemDriver is None:
559 t_fh.FlushCache()
560 else:
561 tt_fh = g.Driver.CreateCopy(temp_tilename, t_fh, 0, g.CreateOptions)
562 tt_fh.FlushCache()
563 tt_fh = None
564
565 t_fh = None
566
567 if os.path.exists(tilename):
568 os.remove(tilename)
569 shutil.move(temp_tilename, tilename)
570
571 if g.Verbose:
572 print(tilename + " : " + str(offsetX) + "|" + str(offsetY) + "-->" + str(width) + "-" + str(height))
573
574
576 OGRDriver = ogr.GetDriverByName(driverName)
577 if OGRDriver is None:
578 print('ESRI Shapefile driver not found')
579 sys.exit(1)
580
581 OGRDataSource = OGRDriver.Open(dsName)
582 if OGRDataSource is not None:
583 OGRDataSource.Destroy()
584 OGRDriver.DeleteDataSource(dsName)
585 if Verbose:
586 print('truncating index ' + dsName)
587
588 OGRDataSource = OGRDriver.CreateDataSource(dsName)
589 if OGRDataSource is None:
590 print('Could not open datasource ' + dsName)
591 sys.exit(1)
592
593 OGRLayer = OGRDataSource.CreateLayer("index", srs, ogr.wkbPolygon)
594 if OGRLayer is None:
595 print('Could not create Layer')
596 sys.exit(1)
597
598 OGRFieldDefn = ogr.FieldDefn(fieldName, ogr.OFTString)
599 if OGRFieldDefn is None:
600 print('Could not create FieldDefn for ' + fieldName)
601 sys.exit(1)
602
603 OGRFieldDefn.SetWidth(256)
604 if OGRLayer.CreateField(OGRFieldDefn) != 0:
605 print('Could not create Field for ' + fieldName)
606 sys.exit(1)
607
608 return OGRDataSource
609
610
611 -def addFeature(TileIndexFieldName, OGRDataSource, location, xlist, ylist):
612 OGRLayer = OGRDataSource.GetLayer()
613 OGRFeature = ogr.Feature(OGRLayer.GetLayerDefn())
614 if OGRFeature is None:
615 print('Could not create Feature')
616 sys.exit(1)
617
618 OGRFeature.SetField(TileIndexFieldName, location)
619 wkt = 'POLYGON ((%f %f,%f %f,%f %f,%f %f,%f %f ))' % (xlist[0], ylist[0],
620 xlist[1], ylist[1], xlist[2], ylist[2], xlist[3], ylist[3], xlist[0], ylist[0])
621 OGRGeometry = ogr.CreateGeometryFromWkt(wkt, OGRLayer.GetSpatialRef())
622 if OGRGeometry is None:
623 print('Could not create Geometry')
624 sys.exit(1)
625
626 OGRFeature.SetGeometryDirectly(OGRGeometry)
627
628 OGRLayer.CreateFeature(OGRFeature)
629 OGRFeature.Destroy()
630
631
634
635
636 -def buildPyramid(g, minfo, createdTileIndexDS, tileWidth, tileHeight, overlap):
637 inputDS = createdTileIndexDS
638 for level in range(1, g.Levels + 1):
639 g.LastRowIndx = -1
640 levelMosaicInfo = mosaic_info(minfo.filename, inputDS)
641 levelOutputTileInfo = tile_info(int(levelMosaicInfo.xsize / 2), int(levelMosaicInfo.ysize / 2), tileWidth, tileHeight, overlap)
642 inputDS = buildPyramidLevel(g, levelMosaicInfo, levelOutputTileInfo, level)
643
644
646 yRange = list(range(1, levelOutputTileInfo.countTilesY + 1))
647 xRange = list(range(1, levelOutputTileInfo.countTilesX + 1))
648
649 OGRDS = createTileIndex(g.Verbose, "TileResult_" + str(level), g.TileIndexFieldName, g.Source_SRS, g.TileIndexDriverTyp)
650
651 for yIndex in yRange:
652 for xIndex in xRange:
653 offsetY = (yIndex - 1) * (levelOutputTileInfo.tileHeight - levelOutputTileInfo.overlap)
654 offsetX = (xIndex - 1) * (levelOutputTileInfo.tileWidth - levelOutputTileInfo.overlap)
655 height = levelOutputTileInfo.tileHeight
656 width = levelOutputTileInfo.tileWidth
657
658 if offsetX + width > levelOutputTileInfo.width:
659 width = levelOutputTileInfo.width - offsetX
660 if offsetY + height > levelOutputTileInfo.height:
661 height = levelOutputTileInfo.height - offsetY
662
663 tilename = getTileName(g, levelMosaicInfo, levelOutputTileInfo, xIndex, yIndex, level)
664
665 feature_only = g.Resume and os.path.exists(tilename)
666 createPyramidTile(g, levelMosaicInfo, offsetX, offsetY, width, height, tilename, OGRDS, feature_only)
667
668 if g.TileIndexName is not None:
669 shapeName = getTargetDir(g, level) + g.TileIndexName
670 copyTileIndexToDisk(g, OGRDS, shapeName)
671
672 if g.CsvFileName is not None:
673 csvName = getTargetDir(g, level) + g.CsvFileName
674 copyTileIndexToCSV(g, OGRDS, csvName)
675
676 return OGRDS
677
678
679 -def getTileName(g, minfo, ti, xIndex, yIndex, level=-1):
680 """
681 creates the tile file name
682 """
683 maxim = ti.countTilesX
684 if ti.countTilesY > maxim:
685 maxim = ti.countTilesY
686 countDigits = len(str(maxim))
687 parts = os.path.splitext(os.path.basename(minfo.filename))
688 if parts[0][0] == "@":
689 parts = (parts[0][1:len(parts[0])], parts[1])
690
691 yIndex_str = ("%0" + str(countDigits) + "i") % (yIndex,)
692 xIndex_str = ("%0" + str(countDigits) + "i") % (xIndex,)
693
694 if g.UseDirForEachRow:
695 frmt = getTargetDir(g, level) + str(yIndex) + os.sep + parts[0] + "_" + yIndex_str + "_" + xIndex_str
696
697 if g.LastRowIndx < yIndex:
698 g.LastRowIndx = yIndex
699 if not os.path.exists(getTargetDir(g, level) + str(yIndex)):
700 os.mkdir(getTargetDir(g, level) + str(yIndex))
701 else:
702 frmt = getTargetDir(g, level) + parts[0] + "_" + yIndex_str + "_" + xIndex_str
703
704 if g.Extension is None:
705 frmt += parts[1]
706 else:
707 frmt += "." + g.Extension
708 return frmt
709
710
717
718
719
720
722 print('Usage: gdal_retile.py ')
723 print(' [-v] [-q] [-co NAME=VALUE]* [-of out_format]')
724 print(' [-ps pixelWidth pixelHeight]')
725 print(' [-overlap val_in_pixel]')
726 print(' [-ot {Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/')
727 print(' CInt16/CInt32/CFloat32/CFloat64}]')
728 print(' [ -tileIndex tileIndexName [-tileIndexField fieldName]]')
729 print(' [ -csv fileName [-csvDelim delimiter]]')
730 print(' [-s_srs srs_def] [-pyramidOnly] -levels numberoflevels')
731 print(' [-r {near/bilinear/cubic/cubicspline/lanczos}]')
732 print(' [-useDirForEachRow] [-resume]')
733 print(' -targetDir TileDirectory input_files')
734
735
736
737 -def main(args=None, g=None):
738 if g is None:
739 g = RetileGlobals()
740
741 gdal.AllRegister()
742
743 if args is None:
744 args = sys.argv
745 argv = gdal.GeneralCmdLineProcessor(args)
746 if argv is None:
747 return 1
748
749
750 i = 1
751 while i < len(argv):
752 arg = argv[i]
753
754 if arg == '-of' or arg == '-f':
755 i += 1
756 g.Format = argv[i]
757 elif arg == '-ot':
758 i += 1
759 g.BandType = gdal.GetDataTypeByName(argv[i])
760 if g.BandType == gdal.GDT_Unknown:
761 print('Unknown GDAL data type: %s' % argv[i])
762 return 1
763 elif arg == '-co':
764 i += 1
765 g.CreateOptions.append(argv[i])
766
767 elif arg == '-v':
768 g.Verbose = True
769 elif arg == '-q':
770 g.Quiet = True
771
772 elif arg == '-targetDir':
773 i += 1
774 g.TargetDir = argv[i]
775
776 if not os.path.exists(g.TargetDir):
777 print("TargetDir " + g.TargetDir + " does not exist")
778 return 1
779 if g.TargetDir[len(g.TargetDir) - 1:] != os.sep:
780 g.TargetDir = g.TargetDir + os.sep
781
782 elif arg == '-ps':
783 i += 1
784 g.TileWidth = int(argv[i])
785 i += 1
786 g.TileHeight = int(argv[i])
787
788 elif arg == '-overlap':
789 i += 1
790 g.Overlap = int(argv[i])
791
792 elif arg == '-r':
793 i += 1
794 ResamplingMethodString = argv[i]
795 if ResamplingMethodString == "near":
796 g.ResamplingMethod = gdal.GRA_NearestNeighbour
797 elif ResamplingMethodString == "bilinear":
798 g.ResamplingMethod = gdal.GRA_Bilinear
799 elif ResamplingMethodString == "cubic":
800 g.ResamplingMethod = gdal.GRA_Cubic
801 elif ResamplingMethodString == "cubicspline":
802 g.ResamplingMethod = gdal.GRA_CubicSpline
803 elif ResamplingMethodString == "lanczos":
804 g.ResamplingMethod = gdal.GRA_Lanczos
805 else:
806 print("Unknown resampling method: %s" % ResamplingMethodString)
807 return 1
808 elif arg == '-levels':
809 i += 1
810 g.Levels = int(argv[i])
811 if g.Levels < 1:
812 print("Invalid number of levels : %d" % g.Levels)
813 return 1
814 elif arg == '-s_srs':
815 i += 1
816 g.Source_SRS = osr.SpatialReference()
817 if g.Source_SRS.SetFromUserInput(argv[i]) != 0:
818 print('invalid -s_srs: ' + argv[i])
819 return 1
820
821 elif arg == "-pyramidOnly":
822 g.PyramidOnly = True
823 elif arg == '-tileIndex':
824 i += 1
825 g.TileIndexName = argv[i]
826 parts = os.path.splitext(g.TileIndexName)
827 if not parts[1]:
828 g.TileIndexName += ".shp"
829
830 elif arg == '-tileIndexField':
831 i += 1
832 g.TileIndexFieldName = argv[i]
833 elif arg == '-csv':
834 i += 1
835 g.CsvFileName = argv[i]
836 parts = os.path.splitext(g.CsvFileName)
837 if not parts[1]:
838 g.CsvFileName += ".csv"
839 elif arg == '-csvDelim':
840 i += 1
841 g.CsvDelimiter = argv[i]
842 elif arg == '-useDirForEachRow':
843 g.UseDirForEachRow = True
844 elif arg == "-resume":
845 g.Resume = True
846 elif arg[:1] == '-':
847 print('Unrecognized command option: %s' % arg)
848 Usage()
849 return 1
850
851 else:
852 g.Names.append(arg)
853 i += 1
854
855 if not g.Names:
856 print('No input files selected.')
857 Usage()
858 return 1
859
860 if (g.TileWidth == 0 or g.TileHeight == 0):
861 print("Invalid tile dimension %d,%d" % (g.TileWidth, g.TileHeight))
862 return 1
863 if (g.TileWidth - g.Overlap <= 0 or g.TileHeight - g.Overlap <= 0):
864 print("Overlap too big w.r.t tile height/width")
865 return 1
866
867 if g.TargetDir is None:
868 print("Missing Directory for Tiles -targetDir")
869 Usage()
870 return 1
871
872
873 if g.UseDirForEachRow and not g.PyramidOnly:
874 leveldir = g.TargetDir + str(0) + os.sep
875 if not os.path.exists(leveldir):
876 os.mkdir(leveldir)
877
878 if g.Levels > 0:
879 startIndx = 1
880 for levelIndx in range(startIndx, g.Levels + 1):
881 leveldir = g.TargetDir + str(levelIndx) + os.sep
882 if os.path.exists(leveldir):
883 continue
884 os.mkdir(leveldir)
885 if not os.path.exists(leveldir):
886 print("Cannot create level dir: %s" % leveldir)
887 return 1
888 if g.Verbose:
889 print("Created level dir: %s" % leveldir)
890
891 g.Driver = gdal.GetDriverByName(g.Format)
892 if g.Driver is None:
893 print('Format driver %s not found, pick a supported driver.' % g.Format)
894 UsageFormat()
895 return 1
896
897 DriverMD = g.Driver.GetMetadata()
898 g.Extension = DriverMD.get(gdal.DMD_EXTENSION)
899 if 'DCAP_CREATE' not in DriverMD:
900 g.MemDriver = gdal.GetDriverByName("MEM")
901
902 tileIndexDS = getTileIndexFromFiles(g)
903 if tileIndexDS is None:
904 print("Error building tile index")
905 return 1
906 minfo = mosaic_info(g.Names[0], tileIndexDS)
907 ti = tile_info(minfo.xsize, minfo.ysize, g.TileWidth, g.TileHeight, g.Overlap)
908
909 if g.Source_SRS is None and minfo.projection:
910 g.Source_SRS = osr.SpatialReference()
911 if g.Source_SRS.SetFromUserInput(minfo.projection) != 0:
912 print('invalid projection ' + minfo.projection)
913 return 1
914
915 if g.Verbose:
916 minfo.report()
917 ti.report()
918
919 if not g.PyramidOnly:
920 dsCreatedTileIndex = tileImage(g, minfo, ti)
921 tileIndexDS.Destroy()
922 else:
923 dsCreatedTileIndex = tileIndexDS
924
925 if g.Levels > 0:
926 buildPyramid(g, minfo, dsCreatedTileIndex, g.TileWidth, g.TileHeight, g.Overlap)
927
928 if g.Verbose:
929 print("FINISHED")
930 return 0
931
932
934 __slots__ = ['Verbose',
935 'Quiet',
936 'CreateOptions',
937 'Names',
938 'TileWidth',
939 'TileHeight',
940 'Overlap',
941 'Format',
942 'BandType',
943 'Driver',
944 'Extension',
945 'MemDriver',
946 'TileIndexFieldName',
947 'TileIndexName',
948 'TileIndexDriverTyp',
949 'CsvDelimiter',
950 'CsvFileName',
951 'Source_SRS',
952 'TargetDir',
953 'ResamplingMethod',
954 'Levels',
955 'PyramidOnly',
956 'LastRowIndx',
957 'UseDirForEachRow',
958 'Resume']
959
961 """ Only used for unit tests """
962 self.Verbose = False
963 self.Quiet = False
964 self.CreateOptions = []
965 self.Names = []
966 self.TileWidth = 256
967 self.TileHeight = 256
968 self.Overlap = 0
969 self.Format = 'GTiff'
970 self.BandType = None
971 self.Driver = None
972 self.Extension = None
973 self.MemDriver = None
974 self.TileIndexFieldName = 'location'
975 self.TileIndexName = None
976 self.TileIndexDriverTyp = "Memory"
977 self.CsvDelimiter = ";"
978 self.CsvFileName = None
979
980 self.Source_SRS = None
981 self.TargetDir = None
982 self.ResamplingMethod = gdal.GRA_NearestNeighbour
983 self.Levels = 0
984 self.PyramidOnly = False
985 self.LastRowIndx = -1
986 self.UseDirForEachRow = False
987 self.Resume = False
988
989
990 if __name__ == '__main__':
991 sys.exit(main(sys.argv))
992