Raster Processing Tutorial

This is a working copy of a tutorial to be presented at FOSS4G 2007 (L-01) that will hopefully also be useful as a web based tutorial in the future.

Objectives

A practical tutorial on processing and preparing raster data for visualization or futher analysis. Exercises will include:

• Exploring your image data
• Format translation
• Optimizing data for MapServer/MapGuide/etc
• Rescaling
• Mosaicing
• Reprojection
• Using virtual files.

Most exercises will be performed using GDAL ​command line utility programs, with a final exercise demonstrating ​Python scripting for specialized processing. It is intended that the workshop would be useful for anyone needing to prepare raster data for use in packages such as ​MapServer, ​MapGuide, ​QGIS, ​GRASS, ​OSSIM or proprietary GIS and imaging applications.

Audience

Suitable for those new to GDAL. Some familiarity with geospatial raster data and coordinate systems is assumed but not strictly required.

Packages

System Requirements: ​FWTools for Windows or Linux (recommended version 1.3.6 or later)

Datasets: Download from ​http://download.osgeo.org/gdal/workshop

Note to self, need to add geogratis and/or geobase license info on download site for cded and CanMatrix data.

Tutorial

format nicely for wiki later - keep plain ascii for now for easy printout at conference.

GDAL Overview
-----------------

5 minute introduction to GDAL. 


1 Exploring Your Image Data
---------------------------

1a) Use gdalinfo to get basic information about the dataset.

  gdalinfo world.png

  Note: 
   - Driver is "PNG"
   - Size is 2048x1024
   - 3 Bands of type Byte. 
   - No coordinate system and bounds are just pixel/line based.

1b) More with gdalinfo.
 
   gdalinfo -mm east.dem

   Note: 
    - Coordinate system is NAD83 (expressed as WKT). 
    - Origin, Pixel Size and bounds are all provided.
    - Band is now of type "Int16"
    - Nodata Value is -32767
    - Unit type is "m" 
    - Min/max elevation values shown (computed due to -mm switch)

1c) Even more with gdalinfo: 

   gdalinfo 092b05.tif

   Note: 
    - Extra bits of TIFF metadata. 
    - Notes the compression 
    - UTM corner coordinates also reported as lat/long. 
    - Color table

1d) Look at your data!

    openev world.png

  Note: No georeferenced coordinates 
   
    openev east.dem

  Note: You can point at pixels and see the elevations, and everything is
  georeferenced.

    openev 092b05.tif

  Right click on the layer in layer dialog and look at the "Coordinate
  System" and "Image Info" tabs for some info similar to gdalinfo.

1e) If you have a minute, try looking at some of the files with QGIS,
    uDig, etc 


2 Simple Format Translation
---------------------------

-- Know your Drivers --

2a) The --formats commandline switch of gdal_translate can be used to see a
    list of available format drivers.  Each format reports if it is read 
    only (ro), read/write (rw) or read/write/update (rw+).

 gdal_translate --formats

2b) The --format commandline switch can be used to query details about a 
    particular driver, including creation options, and permitted data types.

 gdalinfo --format jpeg
 gdal_translate --format png 

-- Translation --

2c) Translations are accomplished with the gdal_translate command.  The
    default output format is GeoTIFF:

 gdal_translate east.dem eastdem.tif

2d) The -of flag is used to select an output format and the -co flag is
    used to specify a creation option:

 gdal_translate -of JPEG -co QUALITY=40 world.png world.jpg

-- Output Datatype --

2e) The -ot switch can be used to alter the output data type.  

 gdal_translate -ot Byte east.dem bytedem.tif

2f) Use gdalinfo to verify data type, and OpenEV to observe that all 
    elevations have been clipped to 255.

-- Rescaling --

2g) The -scale switch can be used to rescale data.  

    gdal_translate -of JPEG -scale east.dem east.jpg 

2h) Explicit control of the input and output ranges is also available,
    and the gdalinfo -mm switch can be used to see pixel min/max values:

    gdalinfo -mm east.dem
    gdal_translate -of JPEG -scale -100 678 0 255 east.dem east.jpg


3 Correcting Metadata
---------------------

-- Set Coordinate System and Bounds --

3a) Use gdal_translate to attach bounds, and coordinate system:

  gdal_translate -a_srs WGS84 -a_ullr -180 90 180 -90 world.png geoworld.tif

3b) Confirm reporting of coordinate system and georeferencing with gdalinfo:

  gdalinfo geoworld.tif

3c) Confirm reporting of georeferenced locations with OpenEV (optional)

  openev geoworld.tif

-- Reset NODATA values --

3d) Use gdal_translate to attach/change the nodata value of a file.

  gdal_translate -a_nodata 0 west.dem zeronodata.tif

3e) Use gdalinfo to verify nodata value differs between the two files.

  gdalinfo west.dem
  gdalinfo zeronodata.tif


Reprojecting
------------

For this process we assume that geoworld.tif has been properly created with
bounds and coordinate system info in the previous steps.  

a) The gdalwarp command can be used to reproject images.  Here we reproject
   the WGS84 geographic image to the Mercator projection:

   gdalwarp -t_srs '+proj=merc +datum=WGS84' geoworld.tif mercator.tif

   Compare the images with OpenEV

   openev mercator.tif
   openev geoworld.tif

b) Here we reproject to the Ortho projection.  

   gdalwarp -t_srs '+proj=ortho +datum=WGS84' geoworld.tif ortho.tif
  
   openev ortho.tif

c) Note how the poles are clipped?  This is because the edges at the pole
   can't be reprojected gdalwarp does not read all the data.  We can force
   gdalwarp to read a bunch of surplus data around chunks as one way to 
   resolve this. 

   gdalwarp -wo SOURCE_EXTRA=125 -t_srs '+proj=ortho +datum=WGS84' 
             geoworld.tif ortho.tif

d) gdalwarp can also be cause to treat particular values as nodata and to 
   produce alpha values in the output.  In this example we cause oceans to
   be treated as transparent, and generate alpha in the output. 

   gdalwarp -wo SOURCE_EXTRA=125 -srcnodata "11 10 50" -dstalpha
             -t_srs '+proj=ortho +datum=WGS84' 
             geoworld.tif ortho.tif

   Hmm, it seems the water isn't quite as uniformly blue as I hoped.  It
   must have lived as a lossily compressed image at one point!


Mosaicing
---------

a) gdal_merge.py is a python script that can be used for simple mosaicing
   tasks.  Mosaic the east.dem and west.dem into a single file:

   gdal_merge  east.dem west.dem -o mergeddem.tif

   For extra credit, open gdal_merge.py in notepad and review.  It's scary,
   but not as scary as you might think!

   notepad "C:\Program Files\FWTools-1.3.6\pymod\gdal_merge.py"

b) The same task can be accomplished with gdalwarp.  gdalwarp has a variety 
   of advantages over gdal_merge, but can be slow to merge many files:

   gdalwarp east.dem west.dem warpmerged.tif

c) Now lets mosaic two CanMatrix image for our area:

   gdalwarp 092b05.tif 092b06.tif mosaic1.tif

   Seems successful, but if we look at the image with OpenEV ... a mess!

   openev mosaic1.tif

d) The problem is that the images have different pseudocolor tables.  So 
   lets convert from pseudocolor to RGB and then mosaic.  To speed up 
   later steps we will also reduce the resolution by half.

   gdal_translate -outsize 50% 50% 092b05.tif small5.tif
   gdal_translate -outsize 50% 50% 092b06.tif small6.tif

   pct2rgb small5.tif rgb5.tif
   pct2rgb small6.tif rgb6.tif

   gdalwarp rgb5.tif rgb6.tif mosaic2.tif

   Hmm, better but clearly we have a problem with collar overlaps. 

e) So, lets try stripping the collars off.  We can do this manually by
   digitizing a polygon representing the "valid" area of the map.  This
   can be done in OpenEV, QGIS, etc, but for our purposes we provide a 
   pre-created cutline.  Review it with:

   openev rgb06.tif cutline06.shp

   "Burn away" the inverse of the polygon using the command:
   
gdal_rasterize -i -burn 17 -b 1 -b 2 -b 3 cutline06.shp -l cutline06 rgb6.tif 

   Review the result, we should have "17,17,17" for all outside pixels.

f) Now mosaic again, but request that "17" be treated as a nodata value.

   gdalwarp -srcnodata "17 17 17" rgb5.tif rgb6.tif mosaic3.tif

   A pretty good result, especially if we were to also strip the collar 
   off the first image.  

g) For extra points, rerun the warp with the -dstalpha option to add an
   alpha (transparency) band for nodata areas. 

   gdalwarp -dstalpha -srcnodata "17 17 17" rgb5.tif rgb6.tif mosaic4.tif 


Optimizing  
----------

gdal_translate -co TILED=YES -co COMPRESS=DEFLATE mosaic4.tif tiled.tif

gdaladdo tiled.tif 2 4 8 16 32 64

Thats about it!


Virtual Files
-------------

a) Create an example VRT with gdal_translate.

   gdal_translate -of VRT east.dem east.vrt 

   Inspect east.vrt, for instance with notepad.

b) Skim VRT tutorial, available in Related Pages at www.gdal.org. 

   Compare the elements of the tutorial to the elements of the file.

c) Create VRT from geoworld.tif.
  
   gdal_translate -of VRT geoworld.tif 0to360.vrt

d) Edit the <GeoTransform> tag and change:

   -1.800000e+02 

   to
   
   0.0

   to set the origin to 0 longitude instead of -180 (the dateline).  Then
   check the .vrt file with gdalinfo and OpenEV. 

   gdalinfo 0to360.vrt
   openev 0to360.vrt

   Our extents are right, but the imagery is still wrong!

e) Duplicate and modify this section:

    <SimpleSource>
      <SourceFilename relativeToVRT="1">geoworld.tif</SourceFilename>
      <SourceBand>1</SourceBand>
      <SrcRect xOff="0" yOff="0" xSize="2048" ySize="1024"/>
      <DstRect xOff="0" yOff="0" xSize="2048" ySize="1024"/>
    </SimpleSource>

   to look like this:

    <SimpleSource>
      <SourceFilename relativeToVRT="1">geoworld.tif</SourceFilename>
      <SourceBand>1</SourceBand>
      <SrcRect xOff="0" yOff="0" xSize="1024" ySize="1024"/>
      <DstRect xOff="1024" yOff="0" xSize="1024" ySize="1024"/>
    </SimpleSource>
    <SimpleSource>
      <SourceFilename relativeToVRT="1">geoworld.tif</SourceFilename>
      <SourceBand>1</SourceBand>
      <SrcRect xOff="1024" yOff="0" xSize="1024" ySize="1024"/>
      <DstRect xOff="0" yOff="0" xSize="1024" ySize="1024"/>
    </SimpleSource>
   
  Check the result with openev ... now the hemispheres are appropriately
  swapped!

f) Now lets create an overlapping wrap world going from -180 to 360. 

   gdal_translate -of VRT geoworld.tif overlap.vrt

   Change rasterXSize from 2048 to 3072.

   Leave GeoTransform unchanged.

   Add this SimpleSource for the 180 to 360 portion.
 
    <SimpleSource>
      <SourceFilename relativeToVRT="1">geoworld.tif</SourceFilename>
      <SourceBand>1</SourceBand>
      <SrcRect xOff="0" yOff="0" xSize="1024" ySize="1024"/>
      <DstRect xOff="2048" yOff="0" xSize="1024" ySize="1024"/>
    </SimpleSource>
   

Python Scripting
----------------

Do a demonstration of reading east.dem into a  Numeric array, 
setting all pixel values less than 1 to -100, and saving to a new 
file.