FAQ - Raster

1. Why won't gdalwarp or gdal_merge write to most formats?
2. How to improve gdalwarp perfomance?
3. How to convert a raster to a layer of polygons?
4. Can I use gdal_rasterize to generate non-solid polygons?
5. How do I use gdal_translate to extract or clip a sub-section of a raster?

Why won't gdalwarp or gdal_merge write to most formats?

GDAL supports many raster formats for reading, but significantly less formats for writing. Of the ones supported for writing most are only supported in create copy mode. Essentially this means they have to be written sequentially from a provided input copy of the image to be written. Programs like gdal_merge.py or gdalwarp that write chunks of imagery non-sequentially cannot easily write to these sequential write formats. Generally speaking formats that are compressed, such as PNG, JPEG and GIF are sequential write. Also some formats require information such as the coordinate system and color table to be known at creation time and so these are also sequential write formats.

When you encounter this problem it is generally a good idea to first write the result to GeoTIFF format, and then translate to the desired target format.

To determine which formats support which capabilities, use the --formats switch with pretty much any GDAL utility. Each driver will include either r (read-only), rw (read or sequential write) or rw+ (read, sequential write or random write).

How to improve gdalwarp perfomance?

Briefly: use the warp memory and config cachemax settings. For example gdalwarp --config GDAL_CACHEMAX 500 -wm 500 uses 500MB of RAM for read/write caching, and 500MB of RAM for working buffers during the warp.

For more details see UserDocs/GdalWarp Will increasing RAM increase the speed of gdalwarp?

How to convert a raster to a layer of polygons?

TBD

Can I use gdal_rasterize to generate non-solid polygons?

See ​How gdal_rasterize works in gdal-dev archives.

As Chris Barker suggests, GDAL's rasterization capability is not sophisticated from a render styling point of view. Other tools may be more appropriate if you want to do anything more sophisticated than rasterize the polygons in a single solid color.

Examples of other tools: ​Quantum GIS, ​GRASS, ​MapServer, ​GMT, ​SAGA GIS.

However, if your raster format supports transparency in the Alpha band (RGBA), then you can use gdal_rasterize to burn in fully transparent areas into your image like :

$ gdal_rasterize -b 4 -burn 0 -where your_field=some_value -l your_layer your_vector_file.shp your_raster

How do I use gdal_translate to extract or clip a sub-section of a raster?

Gdal_translate is designed to convert to and from a variety of raster formats, but it can also perform helpful geoprocessing operations during conversion.

If you would like to extract a sub-section of a raster you can use the -srcwin or -projwin options. In gdal terminology these are "subsetting" operations that allow you to select a "subwindow" to copy from the source dataset into the destination dataset.

Here is an example of using gdal_translate on NAIP orthophotography in sid format to select a small subwindow that shows Blakely Island, WA:

$ gdal_translate -projwin 510286 5385025 518708 5373405 ortho_1-1_1n_s_wa055_2006_1.sid naip_ortho_blakely_island.tif

This example uses the -projwin option which accepts the bounding coordinates in the projected coordinates rather than in pixels (-srcwin). Gdal_translate -projwin needs the upper left x coordinate, the upper left y coordinate, the lower right x coordinate, and the lower right y coordinate. The naip imagery in this example is in NAD 83 Utm 10, so to get these bounding coordinates I simply loaded up the index shapefile that comes packaged with naip imagery in Quantum GIS and read the screen coordinates to form my extent.

Note: Currently, clipping raster images using vector extent polygons is not supported but is under discussion (see http://trac.osgeo.org/gdal/ticket/1599). However is is fairly easy to get the extents of a given shapefile and convert those coordinate pairs into the format needed by gdal_translate without manually reading the extents from another application like qgis. Say you have a shapefile named clipping_mask.shp use ogrinfo to get the extents:

• note the use of the pipe and grep command is optional(| grep Extent), but is a slick way to limit the info reported by ogrinfo to just what you need in this case

$ ogrinfo clipping_mask.shp -so -al | grep Extent
# which gives the extent as xMin,yMin, xMax, yMax:
Extent: (268596, 5362330) - (278396, 5376592)
# which is (xMin,yMin) - (xMax,yMax)

Then copy and paste that text to create your gdal_translate clipping command:

# -projwin's ulx uly lrx lry is equivalent to xMin, yMax, xMax, yMin so just switch the Y coordinates
# For the above Extent that would turn into:
$ gdal_translate -projwin 268596 5376592 278396 5362330 src_dataset dst_dataset

.