Recently I was asked to import some geotiffs into a geoserver instance. I also needed to stitch the geotiffs together so I didn't have to create 100 some odd layers in geoserver.
I found a helpful python script here:
https://gist.github.com/philipn/1148693
But it was producing files that were way too big for my geoserver to handle. Inspired, I wrote the following to group a bunch of tiffs together in several composite geotiffs:
import os
import shutil
import sys
import time
import subprocess
from collections import OrderedDict
DIRS = ['/mnt/hgfs/f/imagery/va', '/mnt/hgfs/f/imagery/md']
GROUP_SIZE = 8
def runCmd( cmd ):
#print cmd
return subprocess.call( cmd, shell=True )
def tif(l) : return l.endswith( '.tif' ) and not l.startswith('group')
def sizeof_fmt(num):
for x in ['bytes', 'kb', 'mb', 'gb', 'tb']:
if num < 1024.0:
return "%3.2f %s" % (num, x)
num /= 1024.0
return "%3.2f%s" % (num, ' tb')
def human_time( secs ):
interval = OrderedDict( [("Y", 365*86400),
("M", 30*86400),
("W", 7*86400),
("D", 86400),
("h", 3600),
("m", 60),
("s", 1)])
secs = int(secs)
s = ""
for unit, value in interval.items():
subres = secs / value
if subres > 0:
secs = secs % subres
s += str(subres) + unit
return s
totalFilesProcessed = 0
totalFileCount = sum([len( filter( tif, sorted( os.listdir( dir ) ) ) ) for dir in DIRS])
pStart = time.time()
secsPerFile = 0;
for dir in DIRS:
dStart = time.time()
files = sorted( os.listdir( dir ) );
files = filter( tif, files );
i = 0
rel_i = i * GROUP_SIZE
dirFilesProcessed = 0
dirFileCount = len(files)
while( rel_i < dirFileCount ):
gStart = time.time()
tifsToProcess = [os.path.join(dir,x) for x in files[rel_i:rel_i + GROUP_SIZE]]
finalTif = os.path.join(dir, 'group-%i.tif' % i)
if os.path.exists( finalTif ):
print 'WARNING, %s exists, removing...' % finalTif
os.remove( finalTif )
print 'processing %i files in group %i inside %s:\n\t%s\n' % (len(tifsToProcess), i, dir, '\n\t'.join(tifsToProcess))
print 'uncompressing and nearblack\'ing...'
for t in tifsToProcess:
uc = '%s.uncompressed' % t
if not os.path.exists( uc ):
cmd = '/usr/local/bin/gdal_translate %s %s' % (t, uc)
if runCmd( cmd ) != 0:
sys.exit(1)
cmd = '/usr/local/bin/nearblack -near 20 %s' %uc
if runCmd( cmd ) != 0:
sys.exit(1)
print '\nstitching...'
cmd = '/usr/local/bin/gdalwarp -dstalpha -srcnodata 0 -wo "SKIP_NOSOURCE" --config "GDAL_CACHEMAX=500" -wm=500 %s.uncompressed %s' % ('.uncompressed '.join(tifsToProcess), finalTif)
if runCmd( cmd ) != 0:
sys.exit(1)
print '\nremoving temporary files...'
for t in tifsToProcess:
#print 'removing: %s.uncompressed' % t
os.remove( '%s.uncompressed' % t )
print '\nadding overviews...'
cmd = '/usr/local/bin/gdaladdo %s 2 4 8 16 32' % finalTif
if runCmd( cmd ) != 0:
sys.exit(1)
gStop = time.time()
dirFilesProcessed += len(tifsToProcess)
totalFilesProcessed += len(tifsToProcess)
print '\nproduced %s file in %s: %s' % (sizeof_fmt(os.path.getsize(finalTif)), human_time(gStop - gStart), finalTif)
print '%3.1f%% complete with directory %s. %i / %i files remain' % ((dirFilesProcessed * 100) / float(dirFileCount), dir, dirFileCount - dirFilesProcessed, dirFileCount )
print '%3.1f%% complete with total. %i / %i files remain' % ((totalFilesProcessed * 100) / float(totalFileCount), totalFileCount - totalFilesProcessed, totalFileCount)
secsPerFile = int((gStop - pStart) / (totalFilesProcessed))
print 'estimated time remaining for directory: %s' % (human_time( (dirFileCount - dirFilesProcessed) * secsPerFile ))
print 'estimated time remaining for all dirs: %s' % (human_time( (totalFileCount - totalFilesProcessed) * secsPerFile))
print '\n\n'
i += 1
rel_i = i * GROUP_SIZE
dStop = time.time()
print 'processed dir %s in %s' % (dir, human_time(dStop - dStart))
pStop = time.time()
print 'total runtime: %s' % (human_time(pStop - pStart))
Basically this script uncompresses and runs nearblack on each tiff in a group and then stitches them together using gdalwarp into a single image. You can define a set of directories to operate over and it will look through and create grouped tiffs in each.
KNOWN ISSUES
One issue with this script is it will blindly take the first 8 (in its posted state) files and stitch them together. That may be what you want. But if your tiff filenames are such that the neighbor images may not be next to each other in a directory listing, the script has the potential to generate some really big resulting tiffs.
For example, say you define a GROUP_SIZE of 2 and you have two geoTiffs tile1.tiff and tile2.tiff. If those two imagery tiles are separated by 100 miles, the resulting tiff will include those two tiles plus whitespace in between them. The resulting tiff file size is going to be much larger than the sum of the file sizes of tile1 and tile2.
Even one outlier tile has the potential to really drive up filesize if it is far enough away from the rest of the group.
It would be better to somehow parse the geographic information in the group and only stitch together adjacent tiles. Perhaps you define a max limit to the number of adjacent tiles in a single group. But it would be nice to avoid generating groups with outliers.
