Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

I have a set of orthophoto images which are each fully filled with data, but when placed together the coverage does not form a rectangle so there are nodata areas.

In MapGuide, I need to be able to display these nodata areas as white so that my users don't have to waste black ink when printing.

At large scale, where I display original full resolution images directly, this is no problem. I just set my map's background colour to white, and the areas where there are no orthophotos display the background.

For performance, I need to be able to merge all of these source images into a resampled composite overview image for display at smaller scales where more of the orthophotos are viewable at once.

I've been attempting to use GDAL to merge and resample the overview, but by default it appears to create the resampled composite GeoTIFF tile with black in the nodata areas, and MapGuide does not allow me to set black as transparent on colour rasters.

Is there a way for me to efficiently get what I want?

I have provided the answer I ended up with using GDAL, but would love to see solutions using other image processing utilities and GIS applications, both open source and proprietary.

share|improve this question
add comment

3 Answers

up vote 7 down vote accepted

The easiest way for me to deal with this problem was to use the GDAL Virtual Format. This format allowed me to treat the entire set of images as a single image object, and transform it in three relatively simple steps.

Creating the Virtual Dataset

GDAL (including Tamas Szekeres' GISInternals Windows binaries and recent versions of OSGeo4W) includes a utility called gdalbuildvrt which can be used to build an initial virtual dataset.

One simple way of using this is to add all of your images to a text file, and then use that text file as an input to gdalbuildvrt. Here's an example (you'll need to put the second command back on one line):

dir /b *.tif > my_images.txt
gdalbuildvrt 
  -hidenodata 
  -vrtnodata "255 255 255" 
  -resolution highest 
  -input_file_list my_images.txt 
  my_image.vrt

This will leave you with an XML file which you can treat as a single image for all GDAL operations. It also internally represents nodata as white, but hides the nodata definition from tools reading from it.

Creating the Resampled Overview

Next, you will perform the resampling and output of the overview image. You can do this with either gdal_translate or gdalwarp. For either of these, remember that the resultant size will be width * height * 3 (number of 8 bit bands) Bytes. If this will be larger than 4GB, you will want to look at the GeoTIFF options for the syntax to specify BigTIFF as your output (-co "BIGTIFF=YES").

For gdal_translate, you will need to determine the virtual image's dimensions using the handy gdalinfo command. Take these dimensions and divide each by some consistent factor to determine the output width and height of your file in pixels.

The command will looks something like (on one line):

gdal_translate
  -outsize 53120 14000
  -co "TILED=YES"
  -co "PROFILE=GEOTIFF"
  -co "BLOCKXSIZE=256"
  -co "BLOCKYSIZE=256"
  my_image.vrt
  my_image.tif

For gdalwarp, you'll need to know the resultant pixelsize; in this case I'm using .5 metre. You'll also want to make a call on resampling method. I prefer cubicspline for orthophoto overviews. It's softer, but you're not going to be using these down to full resolution and in my experience it creates a more compressible image if you're using something like JPEG or ECW.

gdalwarp 
  -r cubicspline 
  -of GTiff 
  -dstnodata "255 255 255" 
  -tr 0.5 0.5 
  -co "PROFILE=GEOTIFF" 
  -co "BIGTIFF=YES" 
  -co "TILED=YES" 
  -co "BLOCKXSIZE=256"
  -co "BLOCKYSIZE=256"
  my_image.vrt 
  my_image.tif

You may also want to consider using JPEG compression options for these resampled GeoTIFF overviews; it shrinks the output file considerably with (according to Frank) only a marginal performance penalty.

  -co "COMPRESS=JPEG" 
  -co "JPEG_QUALITY=80" 
  -co "PHOTOMETRIC=YCBCR"

Overviews

You will also want to run the handy gdaladdo command over the resultant image to build internal "pyramids", so that requests for lower resolutions than the full image dimensions can be met with a subset of data. The increase in performance is more than worth the disk space in most cases. You'll want to play around with the levels that you use here; for very large images you may be able to drop some. The gdaladdo command looks something like this:

gdaladdo 
  -r average 
  my_image.tif 
  2 4 8 16 32 64 128 256

I would suggest experimenting with these levels for optimal performance. You may find that a different resampling interval is better for your applicaition or, based on your image size, that you can drop some of the higher numbers (or that more are needed)

Also, if you are generating an external overview (using the -ro option) consider adding the JPEG compression configuration lines:

  --config COMPRESS_OVERVIEW JPEG 
  --config PHOTOMETRIC_OVERVIEW YCBCR 
  --config INTERLEAVE_OVERVIEW BAND 

(I believe that these are inherited from the parent GeoTIFF for embedded overviews)

Notes

When faced with this problem, I asked on the #gdal channel on freenode.irc.net. This is an amazing resource, and I am thoroughly indebted to Howard Butler, Frank Warmerdam, and Even Rouault for helping me out with this.

share|improve this answer
    
Thank you for being complete, for taking the time, and I'm sure it was time, to cover the whole workflow and not just the specific bit which answered the question. –  matt wilkie Jul 23 '10 at 4:40
1  
it would be interesting to see if gdalsetnull.py works on VRT's, then "edit the VRT file in a text editor, adding a <NoDataValue /> entity" would be unnecessary. –  matt wilkie Jul 23 '10 at 4:44
1  
@matt wilkie, it was definitely a bit of work rewriting and formatting for this site, but I also covered this previously as part of a larger workflow in my blog: jasonbirch.com/nodes/2009/08/11/290/fwtools-ftw-gdal That suggestion is definitely worth a try! –  JasonBirch Jul 23 '10 at 4:45
    
general windows cmd note: use caret ^ for a line break that will be joined when executed (e.g. add ^ to the end of each code line example above to keep both legibility and run-ability). Important caveat: never end a file or command line with a caret unless you want to invoke consume all memory –  matt wilkie Mar 12 at 19:43
add comment

Yes, but through trial and error i was able to detmerine that -vrtnodata 255 has the effect of flagging everything thats white as a nodata, not just off map, which gdal2tiles then treats with alpha transparency according to the -a flag.

So you end up with some of your source image set to transparent, in my case the white parts of dashed roads. This isnt terminal but it sure would be nice to be able to specify to gdal2tiles the "no src image" color, either as a result of original voids in the vrt, or as a result of warp. For my set a shade of pale blue would be just the thing.

After even more trial and error it seems hidenodata is the key. I dont know why these tools are so minimally documented. Heres what works for me, gdal 1.8.

gdalbuildvrt test.vrt -vrtnodata "209 231 245" -hidenodata BX*.tif BY*.tif
gdal2tiles -p raster -s nztm.prj test.vrt out
share|improve this answer
add comment

I'm not too familiar with GDAL, but I guess that a method or command is available where you can set a pixel of certain value with another value?

Nothing to do with it, but in plain SQL something like (just to illustrate - this is much more pseudo-code):

UPDATE raster SET pixel = 255 WHERE pixel = NoData;

I would like to know the answer!

share|improve this answer
    
GDAL does include several methods for replacing colours with other colours during processing in a way similar to your pseudocode, and a standalone utility called nearblack which allows you to clean up image edges containing compression artifacts and turn them white if desired. None of these worked well with the volume of data I was working with at the time. –  JasonBirch Jul 23 '10 at 4:16
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.