Based on the advice in "GeoServer on Steroids" I would aim for a mosaic of GeoTiffs. Page 8 clearly states to choose a mosaic when:
- A single file gets too big (inefficient seeks, too much metadata
to read, etc..)
- Single granules can be large
- Use Tiling + Overviews + Compression on granules
but suggests moving to pyramids when
- Tremendously large dataset
- Too many files / too large files
- Need to serve at all scales
- Especially low resolution
So I would use GDAL to split them into a collection of GeoTiffs with tiles and overviews set and then drop them all into an Image Mosaic. Some experimentation could be useful to see what the best size of GeoTiff v. Number of granules is, but I'd err towards keeping it below a 100 if you are doing a lot of full views.
Update
To split up a GeoTiff you can use a script like:
#!/bin/bash
# If I was really bored I'd work out how to parse the gdalinfo output
#Upper Left ( -25.000, 1240025.000) ( 9d22'14.87"W, 60d50'27.75"N)
#Lower Left ( -25.0000000, -25.0000000) ( 7d33'24.37"W, 49d45'57.40"N)
#Upper Right ( 660025.000, 1240025.000) ( 2d48'16.84"E, 60d57'27.18"N)
#Lower Right ( 660025.000, -25.000) ( 1d37' 1.96"E, 49d50'34.37"N)
#Center ( 330000.000, 620000.000) ( 3d 6'26.52"W, 55d28' 7.68"N)
x=0
y=0
mx=700000
my=1300000
step=10000
for ((i=0; i < $mx; i+=$step)); do
for ((j=0; j < $my; j+=$step)); do
nx=$(($i + $step))
ny=$(($j + $step))
echo "gdal_translate -projwin" $nx $ny $i $j
done
done
You can try to use gdalbuildvrt (http://www.gdal.org/gdalbuildvrt.html) instead of gdal_merge.py. As the result the xml-like file describing virtual raster created. GDAL and QGIS works with this file.
gdalbuildvrt qqv.vrt -separate qq1.tif qq3.tif
The only problem is, that only first band from all files will be added. This simply fixed by editing vrt file.
$ gdalinfo qqv.vrt
Driver: VRT/Virtual Raster
Files: qqv.vrt
/home/bishop/tmp/qq1.tif
/home/bishop/tmp/qq3.tif
Size is 735, 547
Coordinate System is:
GEOGCS["Unknown datum based upon the WGS 84 ellipsoid",
DATUM["Not_specified_based_on_WGS_84_ellipsoid",
SPHEROID["WGS 84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6030"]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433],
AUTHORITY["EPSG","4030"]]
Origin = (513165.703783421427943,5016826.563432835973799)
Pixel Size = (59.969933184855272,-60.001568937759295)
Corner Coordinates:
Upper Left ( 513165.704, 5016826.563) (Invalid angle,Invalid angle)
Lower Left ( 513165.704, 4984005.705) (Invalid angle,Invalid angle)
Upper Right ( 557243.605, 5016826.563) (Invalid angle,Invalid angle)
Lower Right ( 557243.605, 4984005.705) (Invalid angle,Invalid angle)
Center ( 535204.654, 5000416.134) (Invalid angle,Invalid angle)
Band 1 Block=128x128 Type=Byte, ColorInterp=Undefined
Min=0.000 Max=255.000
NoData Value=nan
Band 2 Block=128x128 Type=Byte, ColorInterp=Undefined
Min=0.000 Max=255.000
NoData Value=nan
Band 3 Block=128x128 Type=Byte, ColorInterp=Undefined
Min=0.000 Max=255.000
NoData Value=nan
Band 4 Block=128x128 Type=Int16, ColorInterp=Undefined
Raster VRT format tutorial: http://www.gdal.org/gdal_vrttut.html
Best Answer
New raster tables must be added into existing GeoPackage as subdatasets as documented in the GDAL driver manual page http://www.gdal.org/drv_geopackage_raster.html. Have a look at the fifth example
gdal_translate -of GPKG new.tif existing.gpkg -co APPEND_SUBDATASET=YES -co RASTER_TABLE=new_table
There is no user interface in QGIS for giving the extra parameters so you must edit the automatically generated gdal_translate command manually. For example if you have an existing GeoPackage "test.gpkg" and you want to add image "Img_Sample.tif" as a new raster table named as "Img_Sample_2" edit the command to look like this:
If you have lots of images to convert i would suggest to forget QGIS and graphical user interface and use gdal_translate directly from command line with a little script.