gdalwarp with alpha band masking

[Jan Hartmann]

Thanks Chris (and Frank). This is exacly the scenario for my historical
cadastral mapping project. Thanks for sharing!

Jan


Jan Hartmann
Department of Geography
University of Amsterdam


Chris Hodgson wrote:

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


> 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.
>