Example n° 3: Serving a large number of Cartographic Black/White GeoTiff with Palette

In this example there is a group of Cartographic Black/White images from “CARTA TECNICA DELLA REGIONE TOSCANA”. The purpose of this example is to describe how the GDAL commands may be used for merging the input files provided.

Note

On a Windows machine you can set-up a shell with all GDAL Utilities, running directly the file OSGeo4W.bat under the %TRAINING_ROOT% folder.

Note

Data have the same pixel resolution and same Coordinate Reference System EPSG:25832. Also each pixel is stored into a single bit.

  1. Navigate to the workshop directory $TRAINING_ROOT/data/user_data/gdal_processing_data (on Windows %TRAINING_ROOT%\data\user_data\gdal_processing_data) and find the CTR_data directory.
  2. Navigate inside the CTR_data directory with the GDAL shell.

Note

The following operations must be executed from the shell inside the selected directory.

  1. Calling the gdalinfo command on an image for retrieving the associated information:

    gdalinfo 20E27_1994.TIF
    

    And the result is:

    Driver: GTiff/GeoTIFF
    Files: 20E27_1994.TIF
               20E27_1994.TFW
    Size is 16050, 14050
    Coordinate System is `'
    GeoTransform =
      600769.026848671, 0.1, 7.3789937e-007
      4863785.940434861, -8.172141e-008, -0.1
    Metadata:
      TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
      TIFFTAG_SOFTWARE=IrfanView
      TIFFTAG_XRESOLUTION=72
      TIFFTAG_YRESOLUTION=72
    Image Structure Metadata:
      COMPRESSION=PACKBITS
      INTERLEAVE=BAND
      MINISWHITE=YES
    Corner Coordinates:
    Upper Left  (  600769.027, 4863785.940)
    Lower Left  (  600769.037, 4862380.940)
    Upper Right (  602374.027, 4863785.939)
    Lower Right (  602374.037, 4862380.939)
    Center      (  601571.532, 4863083.440)
    Band 1 Block=16050x4 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
    

    From gdalinfo it is possible to note:

    • No CRS definition. An image without CRS cannot be displayed on GeoServer.
    • Color Interpretation as palette.
    • A GeoTransformation matrix is associated.
    • Tiles Striped (16050x4).
    • Packbits Compression.
  2. Executing the following commands on the tiff images:

    gdalwarp -t_srs EPSG:25832 -of vrt $f ${f:0:10}_temp.vrt
    
    gdal_translate -expand gray -of vrt ${f:0:10}_temp.vrt ${f:0:10}.vrt
    

    The first operation sets the CRS to each image and creates a VRT file, for reducing space occupation. Also the use of gdalwarp internally performs the GeoTransformation associated to the image.

    The second operation is used for changing the color interpretation from palette to gray. This operation is done because in the final steps other grey levels will be added with the interpolation. The expansion to the gray color interpretation leads to an expansion of the pixel dimension from 1 to 8 bits.

    Note

    The expand gray option does not create a multi banded image but only a single banded one.

    Note

    If the user wants to keep the palette, then can go directly to the Optional elaboration without expanding the Palette paragraph.

    These operations must be executed inside a script:

    Linux:

    #!/bin/bash
    FILES="*.TIF"
    echo start
    for f in $FILES
    do
            echo $f
            gdalwarp -t_srs EPSG:25832 -of vrt $f ${f:0:10}_temp.vrt
            gdal_translate -expand gray -of vrt ${f:0:10}_temp.vrt ${f:0:10}.vrt
    done
    echo stop
    

    Windows:

    for /R %%f in (*.tif) do (
            gdalwarp -t_srs EPSG:25832 -of vrt %%~f %%~f_temp.vrt
            gdal_translate -expand gray -of vrt %%~f_temp.vrt %%~f.vrt
    )
    
  3. Listing of all the VRT files into a single text list with the following command:

    ls *TIF.vrt > list.txt (Linux)
    
    or
    
    dir /b *TIF.vrt > list.txt (Windows)
    
  4. Merging of all the input files with the gdalbuildvrt command:

    gdalbuildvrt -srcnodata 255 -vrtnodata 255 -input_file_list list.txt merged_vrt.vrt
    

    Parameters used:

    • -srcnodata 255 -vrtnodata 255 : definition of the No Data associated with the file.
    • -input_file_list list.txt : definition of input files to elaborate.

    The gdalinfo output on the merged image is:

    Driver: VRT/Virtual Raster
    Files: merged_vrt.vrt
               20E27_1994.TIF.vrt
    
               ~
    
               20E60_1995.TIF.vrt
    Size is 50052, 62047
    Coordinate System is:
    PROJCS["ETRS89 / UTM zone 32N",
            GEOGCS["ETRS89",
                    DATUM["European_Terrestrial_Reference_System_1989",
                            SPHEROID["GRS 1980",6378137,298.257222101,
                                    AUTHORITY["EPSG","7019"]],
                            TOWGS84[0,0,0,0,0,0,0],
                            AUTHORITY["EPSG","6258"]],
                    PRIMEM["Greenwich",0,
                            AUTHORITY["EPSG","8901"]],
                    UNIT["degree",0.0174532925199433,
                            AUTHORITY["EPSG","9122"]],
                    AUTHORITY["EPSG","4258"]],
            PROJECTION["Transverse_Mercator"],
            PARAMETER["latitude_of_origin",0],
            PARAMETER["central_meridian",9],
            PARAMETER["scale_factor",0.9996],
            PARAMETER["false_easting",500000],
            PARAMETER["false_northing",0],
            UNIT["metre",1,
                    AUTHORITY["EPSG","9001"]],
            AXIS["Easting",EAST],
            AXIS["Northing",NORTH],
            AUTHORITY["EPSG","25832"]]
    Origin = (600768.734234663190000,4863785.940434861000000)
    Pixel Size = (0.100000372821407,-0.100000372821407)
    Corner Coordinates:
    Upper Left  (  600768.734, 4863785.940) ( 10d15'18.69"E, 43d55'13.06"N)
    Lower Left  (  600768.734, 4857581.217) ( 10d15'14.46"E, 43d51'51.99"N)
    Upper Right (  605773.953, 4863785.940) ( 10d19' 3.07"E, 43d55'10.54"N)
    Lower Right (  605773.953, 4857581.217) ( 10d18'58.64"E, 43d51'49.47"N)
    Center      (  603271.344, 4860683.579) ( 10d17' 8.72"E, 43d53'31.28"N)
    Band 1 Block=128x128 Type=Byte, ColorInterp=Gray
      NoData Value=255
    
  5. Transforming from VRT to GeoTiff with gdal_translate:

    gdal_translate -a_nodata none -co "BLOCKXSIZE=512" -co "BLOCKYSIZE=512" -co "TILED=YES" -co "BIGTIFF=YES" -co "COMPRESS=DEFLATE" merged_vrt.vrt merged_tif.tif
    

    The various input parameters are:

    • -a_nodata none : avoid setting 255 as No Data for a better image optimization.

    • -co “BLOCKXSIZE=512” -co “BLOCKYSIZE=512” -co “TILED=YES” : definition of the tile dimensions.

    • -co “BIGTIFF=YES” -co “COMPRESS=DEFLATE” : definition of the compression method.

      Note

      BIGTIFF=YES must be set for big images because when compression is used, by default gdal_translate is not able to check if the final image is a BigTiff or not.

    From gdalinfo:

    Driver: GTiff/GeoTIFF
    Files: merged_tif.tif
    Size is 50052, 62047
    Coordinate System is:
    PROJCS["ETRS89 / UTM zone 32N",
            GEOGCS["ETRS89",
                    DATUM["European_Terrestrial_Reference_System_1989",
                            SPHEROID["GRS 1980",6378137,298.2572221010002,
                                    AUTHORITY["EPSG","7019"]],
                            TOWGS84[0,0,0,0,0,0,0],
                            AUTHORITY["EPSG","6258"]],
                    PRIMEM["Greenwich",0],
                    UNIT["degree",0.0174532925199433],
                    AUTHORITY["EPSG","4258"]],
            PROJECTION["Transverse_Mercator"],
            PARAMETER["latitude_of_origin",0],
            PARAMETER["central_meridian",9],
            PARAMETER["scale_factor",0.9996],
            PARAMETER["false_easting",500000],
            PARAMETER["false_northing",0],
            UNIT["metre",1,
                    AUTHORITY["EPSG","9001"]],
            AUTHORITY["EPSG","25832"]]
    Origin = (600768.734234663190000,4863785.940434861000000)
    Pixel Size = (0.100000372821407,-0.100000372821407)
    Metadata:
      AREA_OR_POINT=Area
    Image Structure Metadata:
      COMPRESSION=DEFLATE
      INTERLEAVE=BAND
    Corner Coordinates:
    Upper Left  (  600768.734, 4863785.940) ( 10d15'18.69"E, 43d55'13.06"N)
    Lower Left  (  600768.734, 4857581.217) ( 10d15'14.46"E, 43d51'51.99"N)
    Upper Right (  605773.953, 4863785.940) ( 10d19' 3.07"E, 43d55'10.54"N)
    Lower Right (  605773.953, 4857581.217) ( 10d18'58.64"E, 43d51'49.47"N)
    Center      (  603271.344, 4860683.579) ( 10d17' 8.72"E, 43d53'31.28"N)
    Band 1 Block=512x512 Type=Byte, ColorInterp=Gray
    

    This image can be displayed on GeoServer but a further optimization step could bring to better performances. There could be two ways for optimizing the GeoServer performances:

    • building image overviews.
    • building a pyramid of the image.
  6. (Optional) Optimization methods. Here are described two possible optimizations each of them using a different interpolation type:

    • Creation of the overviews associated to the merged image for having better throughput:

      gdaladdo -r average merged_tif.tif 2 4 8 16 32 64 128
      

      Overviews are reduced views of the input image used by GeoServer for displaying the image at a lower resolutions.

      Parameters used:

      • -r average : definition of the interpolation method.
      • 2 ~ 128 : definition of the overviews level

      And with gdalinfo:

      Driver: GTiff/GeoTIFF
      Files: merged_tif.tif
      Size is 50052, 62047
      
      ~
      
      Band 1 Block=512x512 Type=Byte, ColorInterp=Gray
        Overviews: 25026x31024, 12513x15512, 6257x7756, 3129x3878, 1565x1939, 783x970, 392x485
      

      Then the result can be displayed in GeoServer by configuring the image as a GeoTiff (see Adding a GeoTiff section).

    • (Optional) Creation of a pyramid associated to the merged image and displaying the image on GeoServer with the ImagePyramid plugin (see Advanced Mosaic and Pyramid configuration section).

      For building a pyramid the gdalwarp command must be called several times. The operation to execute on the first image is:

      gdalwarp -r cubic -multi -tr 0.200000745642814 -0.200000745642814 -co BLOCKXSIZE=512 -co BLOCKYSIZE=512 -co TILED=YES -co COMPRESS=DEFLATE merged_tif.tif merged_tif_2.tif
      

      The parameters are:

      • -r cubic : definition of the interpolation method (average interpolation can be used only with GDAL 1.10).
      • -multi : forcing to use multithreading.
      • -tr 0,200000745642814 -0,200000745642814 : definition of the image resolution.

      From gdalinfo on the result image:

      Driver: GTiff/GeoTIFF
      Files: merged_tif_2.tif
      Size is 25026, 31024
      
      ~
      
      Band 1 Block=512x512 Type=Byte, ColorInterp=Gray
      

      Then the same operation, with another value for the resolution must be executed on the result image:

      gdalwarp -r cubic -multi -tr 0.400001491285628 -0.400001491285628 -co BLOCKXSIZE=512 -co BLOCKYSIZE=512 -co TILED=YES -co COMPRESS=DEFLATE merged_tif_2.tif merged_tif_4.tif
      

      These operation must be repeated until the final image has a resolution 128 times lower than that of the initial image.

      Note

      Each call of gdalwarp reduces by half the image resolution.

      After creating the various rasters, they must be saved inside a new directory. This directory must be internally divided into sub-directories numbered from 1 to 7, each of them containing a raster of smaller dimension(going from 1 to 7) and leaving the original raster in the super-directory.

      Then the user can configure the following structure with the ImagePyramid plugin.

  7. Displaying the result on GeoServer:

    ../../_images/ctr_screen_pyr.png

    Result as a pyramid (Zoom on the image for seeing the result).

    ../../_images/ctr_screen_ovr.png

    Result with overviews (Zoom on the image for seeing the result).

Optional elaboration without expanding the Palette

If the user wants to keep the palette the steps to achieve are quite similar.

  1. Executing the following commands on the tiff images:

    gdalwarp -t_srs EPSG:25832 -of vrt $f ${f:0:10}_temp.vrt
    
    gdal_translate -of vrt ${f:0:10}_temp.vrt ${f:0:10}.vrt
    

    These operations must be executed inside a script:

    Linux:

    #!/bin/bash
    FILES="*.TIF"
    echo start
    for f in $FILES
    do
            echo $f
            gdalwarp -t_srs EPSG:25832 -of vrt $f ${f:0:10}_temp.vrt
            gdal_translate -of vrt ${f:0:10}_temp.vrt ${f:0:10}.vrt
    done
    echo stop
    

    Windows:

    for /R %%f in (*.tif) do (
            gdalwarp -t_srs EPSG:25832 -of vrt %%~f %%~f_temp.vrt
            gdal_translate -of vrt %%~f_temp.vrt %%~f.vrt
    )
    
  2. Listing of all the VRT files into a single text list with the following command:

    ls *.vrt > list.txt (Linux)
    
    or
    
    dir /b *.vrt > list.txt (Windows)
    

    Warning

    Delete the _temp.vrt files from the list because they overlap with the final vrt files created.

  3. Merging of all the input files with the gdalbuildvrt command:

    gdalbuildvrt -srcnodata 0 -vrtnodata 0 -input_file_list list.txt merged_vrt.vrt
    

    Parameters used:

    • -srcnodata 0 -vrtnodata 0 : definition of the No Data associated with the file.
    • -input_file_list list.txt : definition of input files to elaborate.

    The gdalinfo output on the merged image is:

    Driver: VRT/Virtual Raster
    Files: merged_vrt.vrt
    20E27_1994.TIF.vrt
    
    ~
    
    20E60_1995.TIF.vrt
    Size is 50052, 62047
    Coordinate System is:
    PROJCS["ETRS89 / UTM zone 32N",
            GEOGCS["ETRS89",
                    DATUM["European_Terrestrial_Reference_System_1989",
                            SPHEROID["GRS 1980",6378137,298.257222101,
                                    AUTHORITY["EPSG","7019"]],
                            TOWGS84[0,0,0,0,0,0,0],
                            AUTHORITY["EPSG","6258"]],
                    PRIMEM["Greenwich",0,
                            AUTHORITY["EPSG","8901"]],
                    UNIT["degree",0.0174532925199433,
                            AUTHORITY["EPSG","9122"]],
                    AUTHORITY["EPSG","4258"]],
            PROJECTION["Transverse_Mercator"],
            PARAMETER["latitude_of_origin",0],
            PARAMETER["central_meridian",9],
            PARAMETER["scale_factor",0.9996],
            PARAMETER["false_easting",500000],
            PARAMETER["false_northing",0],
            UNIT["metre",1,
                    AUTHORITY["EPSG","9001"]],
            AXIS["Easting",EAST],
            AXIS["Northing",NORTH],
            AUTHORITY["EPSG","25832"]]
    Origin = (600768.734234663190000,4863785.940434861000000)
    Pixel Size = (0.100000372821407,-0.100000372821407)
    Corner Coordinates:
    Upper Left  (  600768.734, 4863785.940) ( 10d15'18.69"E, 43d55'13.06"N)
    Lower Left  (  600768.734, 4857581.217) ( 10d15'14.46"E, 43d51'51.99"N)
    Upper Right (  605773.953, 4863785.940) ( 10d19' 3.07"E, 43d55'10.54"N)
    Lower Right (  605773.953, 4857581.217) ( 10d18'58.64"E, 43d51'49.47"N)
    Center      (  603271.344, 4860683.579) ( 10d17' 8.72"E, 43d53'31.28"N)
    Band 1 Block=128x128 Type=Byte, ColorInterp=Palette
      NoData Value=0
      Color Table (RGB with 2 entries)
            0: 255,255,255,255
            1: 0,0,0,255
    
  4. Transforming from VRT to GeoTiff with gdal_translate:

    gdal_translate -co "NBITS=1" -co "BLOCKXSIZE=512" -co "BLOCKYSIZE=512" -co "TILED=YES" -co "BIGTIFF=YES" -co "COMPRESS=DEFLATE" merged_vrt.vrt merged_tif.tif
    

    The various input parameters are:

    • -co “NBITS=1” : sets the bits per pixel to 1, because the Palette contains only 0 or 1.

    • -co “BLOCKXSIZE=512” -co “BLOCKYSIZE=512” -co “TILED=YES” : definition of the tile dimensions.

    • -co “BIGTIFF=YES” -co “COMPRESS=DEFLATE” : definition of the compression method.

      Note

      BIGTIFF=YES must be set for big images because when compression is used, by default gdal_translate is not able to check if the final image is a BigTiff or not.

    From gdalinfo:

    Size is 50052, 62047
    Coordinate System is:
    PROJCS["ETRS89 / UTM zone 32N",
            GEOGCS["ETRS89",
                    DATUM["European_Terrestrial_Reference_System_1989",
                            SPHEROID["GRS 1980",6378137,298.2572221010002,
                                    AUTHORITY["EPSG","7019"]],
                            TOWGS84[0,0,0,0,0,0,0],
                            AUTHORITY["EPSG","6258"]],
                    PRIMEM["Greenwich",0],
                    UNIT["degree",0.0174532925199433],
                    AUTHORITY["EPSG","4258"]],
            PROJECTION["Transverse_Mercator"],
            PARAMETER["latitude_of_origin",0],
            PARAMETER["central_meridian",9],
            PARAMETER["scale_factor",0.9996],
            PARAMETER["false_easting",500000],
            PARAMETER["false_northing",0],
            UNIT["metre",1,
                    AUTHORITY["EPSG","9001"]],
            AUTHORITY["EPSG","25832"]]
    Origin = (600768.734234663190000,4863785.940434861000000)
    Pixel Size = (0.100000372821407,-0.100000372821407)
    Metadata:
      AREA_OR_POINT=Area
    Image Structure Metadata:
      COMPRESSION=DEFLATE
      INTERLEAVE=BAND
    Corner Coordinates:
    Upper Left  (  600768.734, 4863785.940) ( 10d15'18.69"E, 43d55'13.06"N)
    Lower Left  (  600768.734, 4857581.217) ( 10d15'14.46"E, 43d51'51.99"N)
    Upper Right (  605773.953, 4863785.940) ( 10d19' 3.07"E, 43d55'10.54"N)
    Lower Right (  605773.953, 4857581.217) ( 10d18'58.64"E, 43d51'49.47"N)
    Center      (  603271.344, 4860683.579) ( 10d17' 8.72"E, 43d53'31.28"N)
    Band 1 Block=512x512 Type=Byte, ColorInterp=Palette
      NoData Value=0
      Image Structure Metadata:
            NBITS=1
      Color Table (RGB with 2 entries)
            0: 255,255,255,255
            1: 0,0,0,255
    
  5. (Optional) Optimization methods described here are similar to that described above:

    • The overview creation method is equal to that described above.
    • For creating the pyramid the commands to use are the same as described above with the addition of the -co “NBITS=1” command.
  6. Displaying the result on GeoServer:

    ../../_images/ctr_palette_pyr.png

    Result as a pyramid (Zoom on the image for seeing the result).

    ../../_images/ctr_palette_ovr.png

    Result with overviews (Zoom on the image for seeing the result).