gdalwarp with alpha band masking

[Frank Warmerdam]

Chris,

Thanks for walking folks through this.  I'm sure it helps others
get an idea of what is possible, and what some of the issues are. I'll
just add a few comments for everyone.

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

I believe Chris was using this to strip off greyscale colormaps on some
of his files.  The colormaps were just simple greyscale ramps, but sort
of messed up later processing.  It was preferable to just strip them
down to plain greyscale images.

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

The -dstalpha switch is new.  It should appear in GDAL 1.2.4. The
-co ALPHA=YES is a special clue to the GTIFF driver to mark the second
band as alpha at creation time for Greyscale+alpha images.  This is really
a hack to work around limits in GDAL for dynamically marking a band as
alpha.

> 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

I would like to stress that rendering an overview from so many files is
inherently "high overhead".  As Chris notes you should try and produce
at least one overview mosaic at high enough resolution that once you
go "below it" you will only be accessing a few of the full resolution
files (say 5-10).

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

I would add that it might be practical to extract the overview image just
by issuing a WMS request to the server, though some steps might be required
to assign georeferencing to the resulting image.

PS. the greyscale+alpha support in MapServer is in 4.3 now.  The RGB+alpha
support is already in 4.2.  I think generating and using alpha bands can
have many applications for masking rasters with MapServer.  Alpha blending
of rasters is only supported for files read via the MapServe GDAL linkage.

Best regards,
--
---------------------------------------+--------------------------------------
I set the clouds in motion - turn up   | Frank Warmerdam, warmerdam at pobox.com
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush    | Geospatial Programmer for Rent