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.
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
- 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.
Suitable for those new to GDAL. Some familiarity with geospatial raster data and coordinate systems is assumed but not strictly required.
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.
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.