Prepping OS rasters for Geoserver

Prepping Ordnance Survey's raster products for use within Geoserver takes a bit of manipulation to get the best performance and quality out of the image tiles. Combine this with some tweaks to Geoserver and you can get fast, slippy, good looking tiles in your favourite mapping software.

I use GDAL and a batch file to process the tiles. I am also running this in a Windows environment and so installed GDAL using the OSGEO4W installer (the 64-bit version). The first part (1) of the batch file sets the GDAL environment for the script. Then we remove (2a) and create (2b) a directory to hold the output. I keep a list of the tiles to process in a text file. We uncompress (3) the tiles first and write out to a GeoTIFF with a SRS of EPSG:27700. We write out a world file (TFW) for each tile. For performance we create internally tiled files with tiles dimensions of 512px x 512px. We repeat the process (4) but create greyscale images instead.

Once those are done we create some overviews (5 and 6) to make rendering at different scales much quicker. The key here is the resampling method "-r average" and the compression for the overviews "--config COMPRESS_OVERVIEW DEFLATE". After much experimenting I got the best quality with these settings. You can also try using "-r nearest or -r gauss" and "--config COMPRESS_OVERVIEW JPEG".

Lastly, I update the statistics (7 and 8) for each tile using "gdalinfo -stats and then create a .prj file for each tile using a python script. Geoserver uses this information when building the Image Mosaic to correctly locate each tile.

rem 1. CREATE UNCOMPRESSED TILED GTIFF WITH OVERVIEWS

rem Set the root directory for OSGEO4W
set OSGEO4W_ROOT=C:\Apps\OSGeo4W64

rem Convert double backslashes to single
set OSGEO4W_ROOT=%OSGEO4W_ROOT:\\=\%
echo. & echo OSGEO4W home is %OSGEO4W_ROOT% & echo.
set PATH=%OSGEO4W_ROOT%\bin;%PATH%

rem Add application-specific environment settings
for %%f in ("C:\Apps\OSGeo4W64\etc\ini\*.bat") do call "%%f"

rem 2a. REMOVE EXISTING DIRECTORY
RD /Q /S C:\PROCESSING\Rasters\VML_COL_2uncomp

rem 2b. NOW MAKE A DIRECTORY FOR THIS PROCESS
mkdir C:\PROCESSING\Rasters\VML_COL_2uncomp

rem FOR ALL TILES IN VML_COL_NO.txt

rem 3. UNCOMPRESS COLOUR TILES FIRST
FOR /F %%I IN (C:\PROCESSING\RASTERS\VML_COL_process\VML_COL_NO.txt) DO gdal_translate -of GTiff -a_srs EPSG:27700 -co TFW=YES -co TILED=YES -co BLOCKXSIZE=512 -co BLOCKYSIZE=512 ..\VML_COL\%%I.tif ..\VML_COL_2uncomp\%%I.tif

rem 4. UNCOMPRESS GREYSCALE TILES SECOND
FOR /F %%I IN (C:\PROCESSING\RASTERS\VML_COL_process\VML_COL_NO.txt) DO gdal_translate -of GTiff -a_srs EPSG:27700 -expand gray -co TFW=YES -co TILED=YES -co BLOCKXSIZE=512 -co BLOCKYSIZE=512 ..\VML_COL\%%I.tif ..\VML_COL_2uncompGS\%%I.tif

rem 5. CREATE COLOUR OVERVIEWS
FOR /F %%I IN (C:\PROCESSING\RASTERS\VML_COL_process\VML_COL_NO.txt) DO gdaladdo -r average --config COMPRESS_OVERVIEW DEFLATE --config GDAL_TIFF_OVR_BLOCKSIZE 512 --config PHOTOMETRIC_OVERVIEW YCBCR --config INTERLEAVE_OVERVIEW PIXEL ..\VML_COL_2uncomp\%%I.tif 2 4 8 16 32

rem 6. CREATE GREYSCALE OVERVIEWS
FOR /F %%I IN (C:\PROCESSING\RASTERS\VML_COL_process\VML_COL_NO.txt) DO gdaladdo -r average --config COMPRESS_OVERVIEW DEFLATE --config GDAL_TIFF_OVR_BLOCKSIZE 512 --config PHOTOMETRIC_OVERVIEW YCBCR --config INTERLEAVE_OVERVIEW PIXEL ..\VML_COL_2uncompGS\%%I.tif 2 4 8 16 32

rem 7. ADD STATISTICS AND PRJ TO COLOUR RASTERS
cd C:\PROCESSING\RASTERS\VML_COL_2uncomp
FOR /F %%I IN ('dir /b *.tif') DO gdalinfo -stats %%I
python ..\VML_COL_process\create_prj.py

rem 8. ADD STATISTICS AND PRJ TO GREYSCALE RASTERS
cd C:\PROCESSING\RASTERS\VML_COL_2uncompGS
FOR /F %%I IN ('dir /b *.tif') DO gdalinfo -stats %%I
python ..\VML_COL_process\create_prj.py

The Python script to create the .prj file for each tile. Save as create_prj.py

import glob

content = 'PROJCS["OSGB 1936 / British National Grid",GEOGCS["OSGB 1936", DATUM["OSGB_1936",SPHEROID["Airy 1830",6377563.396,299.3249646,AUTHORITY["EPSG","7001"]],AUTHORITY["EPSG","6277"]],PRIMEM["Greenwich",0, AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328, AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4277"]], UNIT["metre",1, AUTHORITY["EPSG","9001"]],PROJECTION["Transverse_Mercator"], PARAMETER["latitude_of_origin",49],PARAMETER["central_meridian",-2],PARAMETER["scale_factor",0.9996012717],PARAMETER["false_easting",400000], PARAMETER["false_northing",-100000],AUTHORITY["EPSG","27700"], AXIS["Easting",EAST],AXIS["Northing",NORTH]]'

tifs = glob.glob('*.tif') 

for tif in tifs:  
prj = tif.split('.')[0] + '.prj'  
file = open(prj,'w')  
file.writelines(content)  
file.close()