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. How can I create a blank raster based on a vector files extents for use with gdal_rasterize?
5. Can I use gdal_rasterize to generate non-solid polygons?
6. 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

How can I create a blank raster based on a vector files extents for use with gdal_rasterize?

Currently the gdal_rasterize utility cannot create a blank raster based a vector datasource although this feature has been considered (Ticket #1599). Until that feature is added you'll need to find a way to generate a raster that matches the extent of your vector layer. You can then use that blank raster to burn in features from your vector layer.

One approach is to use gdal_translate to clip out and flatten an existing raster file that covers your area of interest. See this ​example of syntax.

Another approach is to use one of the language bindings to GDAL to create a blank raster from scratch. The tricky part is understanding the GeoTransform syntax. This python snippet shows how to read in a shapefile and output a tiff that matches the shapefile extents

• Note: you'll need to modify the px variable based on your desired resolution and raster dimensions

#!/usr/bin/env python

from osgeo import gdal
from osgeo import osr
from osgeo import ogr
import numpy

shp = 'test.shp'
tiff = 'test.tif'
px = .001
tiff_width = 7850
tiff_height = 3500

# Import vector shapefile
vector = ogr.GetDriverByName('ESRI Shapefile')
src_ds = vector.Open(shp)
src_lyr = src_ds.GetLayerByIndex(index=0)
src_extent = src_lyr.GetExtent()

# Create new raster layer with 4 bands
raster = gdal.GetDriverByName('GTiff')
dst_ds = raster.Create( tiff, tiff_width, tiff_height, 4, gdal.GDT_Byte)

# Create raster GeoTransform based on upper left corner and pixel resolution
raster_transform = [src_extent[0], px, 0.0, src_extent[3], 0.0, -px]
dst_ds.SetGeoTransform( raster_transform )

# Get projection of shapefile and assigned to raster
srs = osr.SpatialReference()
srs.ImportFromWkt(src_lyr.GetSpatialRef().__str__())
dst_ds.SetProjection( srs.ExportToWkt() )

# Create blank raster with fully opaque alpha band
zeros = numpy.zeros( (tiff_height, tiff_width), numpy.uint8 )
dst_ds.GetRasterBand(1).WriteArray( zeros )
dst_ds.GetRasterBand(2).WriteArray( zeros )
dst_ds.GetRasterBand(3).WriteArray( zeros )
opaque = numpy.ones((tiff_height,tiff_width), numpy.uint8 )*255
dst_ds.GetRasterBand(4).WriteArray( opaque )

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

.