gdalwarp with alpha band masking

[Chris Hodgson]

I'd like to share a fun project I've been working on lately with
Mapserver and GDAL, for the interest and potential benefit of others.
BEWARE - Long post follows...

My task was to take a somewhat disorganized collection of orthophotos in
tiff format and setup a WMS to serve a mosaic of these photos for use as
a backdrop to support operators conflating road networks. This at first
seemed like it would be an exceptionally simple task, just get all the
images, run gdaltindex to make a tileindex for mapserver, setup the
mapfile, and it's done. However, there were several problems that needed
to be overcome in particular the quantity of data, and the mixed
projections of the images.

There is 2-4TB of orthophotos to work with. We were able to reduce this
somewhat because we only need orthophotos that overlap with both of the
road networks we are conflating, but it is still more data than we want
to deal with. We decided to solve this by compressing to JPEG-2000
format - at a rate of 10:1, I cannot tell the difference between the
original image and the compressed one. We are using the Kakadu library
and gdal_translate to do this conversion, however, since we are using a
lossy compression we needed to do all of the image processing before the
compression in order to minimize image degradation.

Next problem - not all of the images were in the same projection, most
were in UTM, but I wanted them all in BC Albers. gdalwarp solves the
reprojection problem quite nicely, however, the output images have empty
"collars" around the edges where data is missing due to the reprojected
image not filling the entire bounding box of the image in the new
projection. One way around this would be to put all the tiles together
beforehand with something like gdal_merge.py, then reproject them, then
slice them up into new tiles afterward. However, due to the quantity of
the data, this was not an ideal solution. In fact, you don't even need
to put them all together, all you need to build one tile is itself and
its eight neighbors - however, this requires far more processing than
our final solution for effectively the same result. The final solution
was to ask Frank to add an alpha band mask to the output of gdalwarp, so
that instead of a black collar, there is transparency. With an
additional change to mapserver to automatically use this alpha band
channel (it already did for colour images, just not for greyscale
images), a tiled layer of these reprojected images looks as though it
has been mosaicked together - in fact, mapserver is basically doing the
mosaicking on the fly. Note that the reprojected images all overlap
eachother, but the overlapping areas are mostly transparent.

A couple things which aren't addressed in this solution are color
balance/histogram equalization, and fixing seams with morphing/warping.
The histogram equalization might be easily solved with imageMagik, but
I'm not sure. The seaming is trickier stuff, but I suspect OSSIM might
be able to do this sort of thing. In fact, OSSIM might be able to do all
of this, but I didn't want try to learn a whole new piece of software on
a rush-project.

In the end, here is my processing stream:

- use imageMagik to convert to grayscale if not already so:
convert -colorspace GRAY $in_file $out_file

- use gdalwarp to convert to the right projection, and add transparency
to existing and resulting nodata areas:
gdalwarp -dstalpha -srcnodata 0 -co 'ALPHA=YES' -s_srs 'epsg:26910'
-t_srs 'epsg:42102' -rc $in_file $out_file

- use gdal_translate to convert the reprojected file to jpeg 2000 using
kakadu - options are from the kakadu examples site, they seem to work
well but then I haven't found options that don't work well:
gdal_translate -of JP2KAK -co 'Clayers=20' -co 'Clevels=8' -co
'Corder=RPCL'
                       -co 'Cprecincts={256,256},{256,256},{128,128}'
-co 'ORGgen_plt=yes'
                       -co 'ORGtparts=R' -co 'Cblk={32,32}' -co
'Creversible=no' -co 'quality=10' $in_file out_file

- build a tileindex of all of the files for mapserver:
gdaltindex tiles.shp *.jp2

- setup my layer in mapserver:
  LAYER
    NAME "orthophoto"
    METADATA
      "wms_title" "orthophoto"
    END
    STATUS ON
    TILEINDEX "tiles.shp"
    TILEITEM "Location"
    TYPE RASTER
    PROJECTION
      "init=epsg:42102"
    END
  END

Right now I have 1389 jp2 tiles totalling 45.5GB, and mapserver can
render a screen-size view of all of them in about 15 seconds - the
machine is an Athlon 64 3500+ with 2GB RAM and a pair of 7200RPM 250GB
SATA drives in software RAID 1 on Fedora Core 2 x86_64. The render time
seems to be pretty proportional down to around 1 second... ie. if you
look at half of it, it takes 7 seconds, a fifteenth of it takes 1
second, and basically any zoom level smaller than that takes about 1
second or less. JP2K files (especially with the options I specified
above) are created with multi-resolution and spatial indexing in mind,
so every jp2 file effectively is internally tiled and has many overview
layers, at little to no extra cost in file size - this contributes
greatly to the speed, I'm sure.

I'm still processing data right now, but once I'm all done I will use
gdal_merge.py to create an overview layer. It will be slightly degraded
since I will have to base it on the lossily-compressed jp2s, but the
quality of the overview isn't important. With an overview layer in place
(or maybe two), I expect to be able to provide sub-second response time
to render any bounding box at screen resolution.

Thanks to Frank Warmerdam for adding the new features into GDAL and
Mapserver that made this all possible.
Sorry, but I can't share the final result due to licensing agreements on
the orthophotos.

Hope this was worth the read for some of you :-)

Chris Hodgson
Refractions Research Inc.