We have a collection of black-and-white map tiles as
*.tif images with accompanying
*.tfw files, that we'd like to merge into one single large map and cut with a polygon (turning the area outside the polygon transparent).
gdalinfo output for them looks like this:
Driver: GTiff/GeoTIFF Files: 702_232.tif 702_232.tfw Size is 4000, 4000 Coordinate System is `' Origin = (702000.000000000000000,233000.000000000000000) Pixel Size = (0.250000000000000,-0.250000000000000) Metadata: TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch) TIFFTAG_XRESOLUTION=96 TIFFTAG_YRESOLUTION=96 Image Structure Metadata: COMPRESSION=LZW INTERLEAVE=BAND MINISWHITE=YES Corner Coordinates: Upper Left ( 702000.000, 233000.000) Lower Left ( 702000.000, 232000.000) Upper Right ( 703000.000, 233000.000) Lower Right ( 703000.000, 232000.000) Center ( 702500.000, 232500.000) Band 1 Block=4000x8 Type=Byte, ColorInterp=Palette Image Structure Metadata: NBITS=1 Color Table (RGB with 2 entries) 0: 255,255,255,255 1: 0,0,0,255
So this seem to be 1-band TIFFs with a color palette.
To merge them, we're first building a VRT with
gdalbuildvrt -a_srs 'EPSG:21781' merged.vrt *.tif
(as you can see, we also set the source SRS, which is missing in the input data)
The individual input files (the tiles) and the VRT display fine (black and white, as intended) in QGIS 2.18.
To build an actual single image file (and at the same time cut the result to the polygon area1), we then use
gdalwarp \ -cutline swissBOUNDARIES3D/FGDB/swissBOUNDARIES3D_1_3_LV03_LN02.gdb \ -cl TLM_KANTONSGEBIET \ -cwhere "NAME='St. Gallen'" \ -dstalpha \ \ -co TILED=YES \ -co COMPRESS=LZW \ -co BIGTIFF=YES \ \ --config GDAL_CACHEMAX=40% \ -multi \ -et 0 \ -wo NUM_THREADS=ALL_CPUS \ -wm 5000 \ \ merged.vrt merged_and_cut.tif
The result is a 2-band TIFF (the 2nd band being alpha, required for the transparency after cutting) without a color palette. When displayed in QGIS, it appears in bright green instead of black and white.
If we add
-setci to the options for
gdalwarp we get a 2-band TIFF with a (sensible looking, even though much larger) color palette. But in QGIS 2.18, it also appears in bright green by default, at least until go to layer properties → style and manually set "Render type" to "Palette".
As the original files and the VRT work fine, I wondered whether we can use the "Raw mode" documented at http://gdal.org/frmt_gtiff.html to preserve the first band's color interpretation, but there's no indication on how to use with
GTIFF_RAW:merged.vrt in the the
gdalwarp invocation results in an error about the file not being found.)
How can I preserve the color interpretation with
gdalwarp? (So that it is correctly understood by QGIS and other software.)
For licensing reasons I think I cannot share the tiles I want to merge, but I've produced these two example tiles (with nonsensical content) that have the same properties and show the same behavior as the original tiles:
- Inputs (TIFFs and TFWs): ownCloud directory | download as ZIP archive
gdalbuildvrt -a_srs 'EPSG:21781' merged.vrt ../input/*.tifresults in
merged.vrt(2 KB) (direct download)
gdalwarp -dstalpha -co TILED=YES -co COMPRESS=LZW merged.vrt merged.tifresults in
merged.tif(17.6 MB) (direct download)
1 Thus why we use
gdalwarp. If we wouldn't cut to a polygon,
gdaltranslate instead of
gdalwarp would suffice to materialize the VRT to a single TIFF.