geoprocessing - How to convert GeoTIFF to grayscale and add gaussian blur?

I have a giant GeoTIFF I'd like to convert to grayscale and add a gaussian blur to. As far as I know, GDAL can strip two of the three color bands and call it grayscale, which might work, but is there a better way? Imagemagick has both grayscale and gaussian blur functionality, but can not process georeferenced TIFFs (as far as I know). Any suggestions how I should go about doing this?

Answer

GDAL has a wonderful file format called VRT, which is an XML wrapper around one or more raster files.

One feature of VRTs is their ability to encode square convolution kernels for any given band. It does involve playing around with XML in a text editor (or programatically), but if you're already used to the GDAL tools, it shouldn't be too hard.

To illustrate, I took this image of some Ordnance Survey data from around the old OS building in Southampton:

And ran gdalbuildvrt to generate an initial VRT file:

gdalbuildvrt shirley.vrt shirley.tif

Then I replaced the elements with elements, adding in the coefficients of a 5 x 5 Gaussian kernel:

  shirley.tif
  1
  
  
  
  
    5
    0.0036630037 0.0146520147 0.0256410256 0.0146520147 0.0036630037

           0.0146520147 0.0586080586 0.0952380952 0.0586080586 0.0146520147
           0.0256410256 0.0952380952 0.1501831502 0.0952380952 0.0256410256
           0.0146520147 0.0586080586 0.0952380952 0.0586080586 0.0146520147
           0.0036630037 0.0146520147 0.0256410256 0.0146520147 0.0036630037
    
  

Then ran gdal_translate to convert to a TIFF:

gdal_translate -co TILED=YES shirley_gauss.vrt shirley_gauss.tif

Which gives me this image:

With its georeferencing data intact.

For the greyscale part, I suggest you use Quantum GIS and its good (if currently slightly quirky) Raster Calculator. Simply load up your blurred image, select Raster | Raster Calculator... and use the following expression:

(shirley_gauss@1 * 0.3) + (shirley_gauss@2 * 0.59) + (shirley_gauss@3 * 0.11)

Loading that image into QGIS, gives me:

Other coefficients can be used for the greyscale conversion, but those are a good starting point.

Also, it shouldn't matter whether you blur first then reduce, or the other way round.