FAQ - Raster
- How to make GeoTIFF from non-geospatial raster?
- Why won't gdalwarp or gdal_merge write to most formats?
How to improve
- How to convert a raster to a layer of polygons?
- How can I create a blank raster based on a vector files extents for use with gdal_rasterize? (GDAL < 1.8.0)
- Can I use gdal_rasterize to generate non-solid polygons?
- How do I use gdal_translate to extract or clip a sub-section of a raster?
- How to create or modify an image color table ?
- How to generate cloud optimized GeoTIFF files
How to make GeoTIFF from non-geospatial raster?
gdal_translate -a_ullr ulx uly lrx lry src_dataset dst_dataset
or with SRS assignment:
gdal_translate -a_srs EPSG:4326 -a_ullr ulx uly lrx lry src_dataset dst_dataset
In fact, src_dataset and dst_dataset can be any raster datasets supported by GDAL with creation capability.
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
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?
The gdal_polygonize.py python script can be used to convert raster layers to polygon layers since GDAL 1.6. It will produce very detailed polygons around each area of equal pixel value.
How can I create a blank raster based on a vector files extents for use with gdal_rasterize? (GDAL < 1.8.0)
The method detailed in this paragraph is no longer necessary with GDAL >= 1.8.0. See ticket #3505
Before GDAL 1.8.0, 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 geo-tiff that matches the shapefile extents. It deals with two conceivable alternatives: a) the calculation of the output pixel resolution of the geo-tiff by using the desired dataset dimensions and b) the calculation of the output dimensions based on desired pixel resolutions:
#!/usr/bin/env python try: from osgeo import gdal from osgeo import osr from osgeo import ogr except: import gdal import osr import ogr # To be adjusted to your needs shp = 'test.shp' tiff = 'test.tif' # Alternative A: Set raster dataset dimensions tiff_width = 7850 tiff_height = 3500 # Alternative B: Set pixel resolution x_res = 0.5 y_res = 0.75 # Import vector shapefile vector = ogr.GetDriverByName('ESRI Shapefile') src_ds = vector.Open(shp) src_lyr = src_ds.GetLayer(0) src_extent = src_lyr.GetExtent() # Alternative A: Calculate pixel resolution from desired dataset dimensions x_res = (src_extent - src_extent) / tiff_width y_res = (src_extent - src_extent) / tiff_height # Alternative B: Calculate raster dataset dimensions from desired pixel resolution tiff_width = int(math.ceil(abs(src_extent - src_extent) / x_res)) tiff_height = int(math.ceil(abs(src_extent - src_extent) / y_res)) # 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 dimensions raster_transform = [src_extent, x_res, 0.0, src_extent, 0.0, -y_res] 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 dst_ds.GetRasterBand(1).Fill(0) dst_ds.GetRasterBand(2).Fill(0) dst_ds.GetRasterBand(3).Fill(0) dst_ds.GetRasterBand(4).Fill(255)
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.
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
How to create or modify an image color table ?
The easiest non-scripting way to create or modify a color palette is to translate your file to VRT (XML) format, and edit the color table in the XML file.
# create the vrt gdal_translate -of VRT your.tif your.vrt ## --- edit colour table here --- # create final fixed image gdal_translate your.vrt your_fixedup.tif
The .vrt file doesn't actually include the image data, it just references back to your.tif (or other supported format). This general approach can be used to modify metadata, color tables, gcps and other things for which there is no easy to do set them with gdal_translate. There are other mechanisms to do this via scripting or c/c++ programming.
<ColorInterp>Palette</ColorInterp> <ColorTable> <Entry c1="0" c2="0" c3="0" c4="255"/> <Entry c1="255" c2="255" c3="255" c4="255"/> <Entry c1="0" c2="151" c3="164" c4="255"/> <Entry c1="203" c2="0" c3="23" c4="255"/> <Entry c1="131" c2="66" c3="37" c4="255"/> <Entry c1="201" c2="234" c3="157" c4="255"/> <Entry c1="137" c2="51" c3="128" c4="255"/> <Entry c1="255" c2="234" c3="0" c4="255"/> <Entry c1="167" c2="226" c3="226" c4="255"/> <Entry c1="255" c2="184" c3="184" c4="255"/> <Entry c1="218" c2="179" c3="214" c4="255"/> <Entry c1="209" c2="209" c3="209" c4="255"/> <Entry c1="207" c2="164" c3="142" c4="255"/> </ColorTable>