tiles - GDAL tiling with a lot of MrSID files

I'm trying to generate tiles for a region with a surface of 85k km2. The data source are MrSID files (2734 files). I'm using gdal to achieve this task.

It works fine with a small subset of MrSID, but with all the source files it would take more years than I could live :).

My process for generating tiles is the following:

# create the mosaic
gdalwarp -of GTiff -dstalpha -r cubic [MrSID Files] mosaic.tiff`

# generate tiles from mosaic
gdal2tiles.py -r cubic -s EPSG:23030 -z 3-17 -v mosaic.tiff

At first I've chosen a subset of 30 MrSID and it takes 10 hours and it works. Check here.

After that, I've tried with all the MrSID files and I can't reach to create the mosaic, after two days gdalwarp is at 1% processing the first MrSID.

As a workaround, I've tried to avoid the creation of the big mosaic, I've called gdal2tiles for each MrSID, and merging the shared tiles with imagemagic. The issue with this method is that I need a huge quantity of temporal files (30 MrSID files (90MB) give me around 30GB of temporal files).

Any ideas how can achieve this? Is it possible create such a big mosaic? Do you suggest me another strategy?

MrSID files are like this:

Driver: MrSID/Multi-resolution Seamless Image Database (MrSID)
Files: test/807-1-3.sid
Size is 7355, 4796
Coordinate System is `'

Origin = (309892.000000000000000,4291525.000000000000000)
Pixel Size = (1.000000000000000,-1.000000000000000)
Metadata:
  ICC__Profile=0,0,1,148,65,68,66,69,2,16,0,0,109,110,116,114,71,82,65,89,88,89,90,32,0,8,7,212,0,19,0,2,0,29,0,23,97,99,115,112,77,83,70,84,0,0,0,0,110,111,110,101,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,246,214,0,1,0,0,0,0,211,45,65,68,66,69,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,99,112,114,116,0,0,0,192,0,0,0,36,100,101,115,99,0,0,0,228,0,0,0,119,119,116,112,116,0,0,1,92,0,0,0,20,98,107,112,116,0,0,1,112,0,0,0,20,107,84,82,67,0,0,1,132,0,0,0,14,116,101,120,116,0,0,0,0,40,99,41,32,50,48,48,52,32,65,100,111,98,101,32,83,121,115,116,101,109,115,32,73,110,99,46,0,100,101,115,99,0,0,0,0,0,0,0,28,69,115,99,97,108,97,32,100,101,32,103,114,105,115,101,115,32,45,32,71,97,109,97,32,50,46,50,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,88,89,90,32,0,0,0,0,0,0,243,81,0,1,0,0,0,1,22,204,88,89,90,32,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,99,117,114,118,0,0,0,0,0,0,0,1,2,51,0,0
  IMAGE__COMPRESSION_BLOCK_SIZE=512
  IMAGE__COMPRESSION_GAMMA=2.000000
  IMAGE__COMPRESSION_NLEV=8
  IMAGE__COMPRESSION_VERSION=2,0,0
  IMAGE__COMPRESSION_WEIGHT=4.000000
  IMAGE__CREATION_DATE=Sat Jan 12 13:20:09 2002


  IMAGE__INPUT_FILE_SIZE=35277072.000000
  IMAGE__INPUT_FORMAT=TIFF w/ World File
  IMAGE__INPUT_NAME=H:\Imagenes\MrSid\Vuelo_Americano\comp\807_1_3.TIF
  IMAGE__STATISTICS_MAX=255
  IMAGE__STATISTICS_MEAN=231.585752
  IMAGE__STATISTICS_MIN=1
  IMAGE__STATISTICS_STANDARD_DEVIATION=58.039863
  IMAGE__TARGET_COMPRESSION_RATIO=10.000000
  PShop__ImageResources=

  VERSION=MG2
Corner Coordinates:
Upper Left  (  309892.000, 4291525.000) 
Lower Left  (  309892.000, 4286729.000) 
Upper Right (  317247.000, 4291525.000) 
Lower Right (  317247.000, 4286729.000) 
Center      (  313569.500, 4289127.000) 
Band 1 Block=1024x128 Type=Byte, ColorInterp=Gray
  Minimum=3.000, Maximum=255.000, Mean=232.406, StdDev=53.985
  Overviews: 3678x2398, 1839x1199, 920x600, 460x300, 230x150, 115x75, 58x38, 29x19

Answer

In order to avoid huge temporary files, you can use VRT as format of the mosaic file, whose main advantage is that it's just an XML file that will be created immediately:

:: create the mosaic (optionally setting the target extent and resolution)
gdalbuildvrt [-te xmin ymin xmax ymax] [-tr xres yres] mosaic.vrt *.sid

:: generate tiles from mosaic
gdal2tiles.py -r cubic -s EPSG:23030 -z 3-17 -v mosaic.vrt

In order to improve the time performance of the above workflow, I suggest to use a less expensive resampling method or, alternatively, to reduce the number of zoom levels to render in gdal2tiles.py.