Cheat sheet for GDAL/OGR command-line tools

Cheat sheet for GDAL/OGR command-line geodata tools.

Vector operations

Get vector information

ogrinfo -so input.shp layer-name

Or, for all layers

ogrinfo -al -so input.shp

Print vector extent

ogrinfo input.shp layer-name | grep Extent

List vector drivers

ogr2ogr --formats

Convert between vector formats

ogr2ogr -f "GeoJSON" output.json input.shp

Print count of features with attributes matching a given pattern

ogrinfo input.shp layer-name | grep "Search Pattern" | sort | uniq -c

Clip vectors by bounding box

ogr2ogr -f "ESRI Shapefile" output.shp input.shp -clipsrc <x_min> <y_min> <x_max> <y_max>

Clip one vector by another

ogr2ogr -clipsrc clipping_polygon.shp output.shp input.shp

Reproject vector:

ogr2ogr output.shp -t_srs "EPSG:4326" input.shp

Merge features in a vector file by attribute (“dissolve”)

ogr2ogr -f "ESRI Shapefile" dissolved.shp input.shp -dialect sqlite -sql "select ST_union(Geometry),common_attribute from input GROUP BY common_attribute"

Merge features (“dissolve”) using a buffer to avoid slivers

ogr2ogr -f "ESRI Shapefile" dissolved.shp input.shp -dialect sqlite \
-sql "select ST_union(ST_buffer(Geometry,0.001)),common_attribute from input GROUP BY common_attribute"

Merge vector files:

ogr2ogr merged.shp input1.shp
ogr2ogr -update -append merged.shp input2.shp -nln merged

Extract from a vector file based on query

To extract features with STATENAME ‘New York’,‘New Hampshire’, etc. from states.shp

ogr2ogr -where 'STATENAME like "New%"' states_subset.shp states.shp

To extract type ‘pond’ from water.shp

ogr2ogr -where "type = pond" ponds.shp water.shp

Subset & filter all shapefiles in a directory

Assumes that filename and name of layer of interest are the same…

basename -s.shp *.shp | xargs -n1 -I % ogr2ogr %-subset.shp %.shp -sql "SELECT field-one, field-two FROM '%' WHERE field-one='value-of-interest'"

Extract data from a PostGis database to a GeoJSON file

ogr2ogr -f "GeoJSON" file.geojson PG:"host=localhost dbname=database user=user password=password" \
-sql "SELECT * from table_name"

Get the difference between two vector files

Create a layers.vrt file that looks like:

	<OGRVRTLayer name="file1">
	<OGRVRTLayer name="file2">

Then run:

ogr2ogr -f "ESRI Shapefile" difference.shp layers.vrt -dialect sqlite \
-sql "SELECT ST_Difference(file1.geometry,file2.geometry) AS geometry FROM file1,file2"

This will produce a vector file with the part of file1.shp that doesn’t intersect with file2.shp.

Or, do it all as a one-liner:

echo '<OGRVRTDataSource><OGRVRTLayer name="'$SUBTRACT_FROM_SHP'"><SrcDataSource>'$SUBTRACT_FROM_SHP'.shp</SrcDataSource></OGRVRTLayer><OGRVRTLayer name="'$SUBTRACT_SHP'"><SrcDataSource>'$SUBTRACT_SHP'.shp</SrcDataSource></OGRVRTLayer></OGRVRTDataSource>' | \
ogr2ogr -f "ESRI Shapefile" difference.shp /vsistdin/ -dialect sqlite \

Spatial join:

A spatial join transfers properties from one vector layer to another based on a spatial relationship between the features. Fields from the join layer can be aggregated in the output.

Given a set of points (trees.shp) and a set of polygons (parks.shp) in the same directory, create a polygon layer with the geometries from parks.shp and summaries of some columns in trees.shp:

ogr2ogr -f 'ESRI Shapefile' output.shp parks.shp -dialect sqlite \
-sql "SELECT p.Geometry, id, SUM(t.field1) field1_sum, AVG(t.field2) field2_avg
FROM parks p, 'trees.shp'.trees t WHERE ST_Contains(p.Geometry, t.Geometry) GROUP BY"

Note that features that from parks.shp that don’t overlap with trees.shp won’t be in the new file.

Raster operations

Get raster information

gdalinfo input.tif

List raster drivers

gdal_translate --formats

Force creation of world file (requires libgeotiff)

listgeo -tfw  mappy.tif

Report PROJ.4 projection info, including bounding box (requires libgeotiff)

listgeo -proj4 mappy.tif

Convert between raster formats

gdal_translate -of "GTiff" input.grd output.tif

Convert 16-bit bands (Int16 or UInt16) to Byte type
(Useful for Landsat 8 imagery…)

gdal_translate -of "GTiff" -co "COMPRESS=LZW" -scale 0 65535 0 255 -ot Byte input_uint16.tif output_byte.tif

You can change ‘0’ and ‘65535’ to your image’s actual min/max values to preserve more color variation or to apply the scaling to other band types - find that number with:

gdalinfo -mm input.tif | grep Min/Max

Convert a directory of raster files of the same format to another raster format

basename -s.img *.img | xargs -n1 -I % gdal_translate -of "GTiff" %.img %.tif

Reproject raster:

gdalwarp -t_srs "EPSG:102003" input.tif output.tif

Be sure to add -r bilinear if reprojecting elevation data to prevent funky banding artifacts.

Georeference an unprojected image with known bounding coordinates:

gdal_translate -of GTiff -a_ullr <top_left_lon> <top_left_lat> <bottom_right_lon> <bottom_right_lat> \
-a_srs EPSG:4269 input.png output.tif

Clip raster by bounding box

gdalwarp -te <x_min> <y_min> <x_max> <y_max> input.tif clipped_output.tif

Clip raster to SHP / NoData for pixels beyond polygon boundary

gdalwarp -dstnodata <nodata_value> -cutline input_polygon.shp input.tif clipped_output.tif

Crop raster dimensions to vector bounding box

gdalwarp -cutline cropper.shp -crop_to_cutline input.tif cropped_output.tif

Merge rasters -o merged.tif input1.tif input2.tif


gdalwarp input1.tif input2.tif merged.tif

Or, to preserve nodata values:

gdalwarp input1.tif input2.tif merged.tif -srcnodata <nodata_value> -dstnodata <merged_nodata_value>

Stack grayscale bands into a georeferenced RGB

Where LC81690372014137LGN00 is a Landsat 8 ID and B4, B3 and B2 correspond to R,G,B bands respectively: -co "PHOTOMETRIC=RGB" -separate LC81690372014137LGN00_B{4,3,2}.tif -o LC81690372014137LGN00_rgb.tif

Fix an RGB TIF whose bands don’t know they’re RGB -co "PHOTOMETRIC=RGB" input.tif -o output_rgb.tif

Export a raster for Google Earth

gdal_translate -of KMLSUPEROVERLAY input.tif output.kmz -co FORMAT=JPEG

Raster calculation (map algebra)

Average two rasters: -A input1.tif -B input2.tif --outfile=output.tif --calc="(A+B)/2"

Add two rasters: -A input1.tif -B input2.tif --outfile=output.tif --calc="A+B"


Create a hillshade from a DEM

gdaldem hillshade -of PNG input.tif hillshade.png

Change light direction:

gdaldem hillshade -of PNG -az 135 input.tif hillshade_az135.png 

Use correct vertical scaling in meters if input is projected in degrees

gdaldem hillshade -s 111120 -of PNG input_WGS1984.tif hillshade.png

Apply color ramp to a DEM
First, create a color-ramp.txt file:
(Height, Red, Green, Blue)

	0 110 220 110
	900 240 250 160
	1300 230 220 170
	1900 220 220 220
	2500 250 250 250

Then apply those colors to a DEM:

gdaldem color-relief input.tif color_ramp.txt color-relief.tif

Create slope-shading from a DEM
First, make a slope raster from DEM:

	gdaldem slope input.tif slope.tif 

Second, create a color-slope.txt file:
(Slope angle, Red, Green, Blue)

0 255 255 255
90 0 0 0  

Finally, color the slope raster based on angles in color-slope.txt:

gdaldem color-relief slope.tif color-slope.txt slopeshade.tif

Resample (resize) raster

gdalwarp -ts <width> <height> -r cubic dem.tif resampled_dem.tif

Entering 0 for either width or height guesses based on current dimensions.


gdal_translate -outsize 10% 10% -r cubic dem.tif resampled_dem.tif

For both of these, -r cubic specifies cubic interpolation: when resampling continuous data (like a DEM), the default nearest neighbor interpolation can result in “stair step” artifacts.

Burn vector into raster

gdal_rasterize -b 1 -i -burn -32678 -l layername input.shp input.tif

Extract polygons from raster input.tif -f "GeoJSON" output.json

Create contours from DEM

gdal_contour -a elev -i 50 input_dem.tif output_contours.shp

Get values for a specific location in a raster

gdallocationinfo -xml -wgs84 input.tif <lon> <lat>  

Convert GRIB band to .tif

Assumes data for entire globe in WGS84. Be sure to specify band, or you may end up with a nonsense RGB image composed of the first three.

gdal_translate input.grib -a_ullr -180 -90 180 90 -a_srs EPSG:4326 -b 1 output_band1.tif


Convert KML points to CSV (simple)

ogr2ogr -f CSV output.csv input.kmz -lco GEOMETRY=AS_XY

Convert KML to CSV (WKT)
First list layers in the KML file

ogrinfo -so input.kml

Convert the desired KML layer to CSV

ogr2ogr -f CSV output.csv input.kml -sql "select *,OGR_GEOM_WKT from some_kml_layer"

CSV points to SHP
This section needs retooling
Given input.csv


Make a .dbf table for ogr2ogr to work with from input.csv

ogr2ogr -f "ESRI Shapefile" input.dbf input.csv

Use a text editor to create a .vrt file in the same directory as input.csv and input.dbf. This file holds the parameters for building a full shapefile based on values in the DBF you just made.

  <OGRVRTLayer name="output_file_name">
    <SrcDataSource relativeToVRT="1">./</SrcDataSource>
    <GeometryField encoding="PointFromColumns" x="lon_column" y="lat_column"/>

Create shapefile based on parameters listed in the .vrt

mkdir shp
ogr2ogr -f "ESRI Shapefile" shp/ inputfile.vrt

The VRT file can be modified to give a new output shapefile name, reference a different coordinate system (LayerSRS), or pull coordinates from different columns.

MODIS operations

First, download relevant .hdf tiles from the MODIS ftp site:; use the MODIS sinusoidal grid for reference.

List MODIS Subdatasets in a given HDF (conf. the MODIS products table)

gdalinfo longFileName.hdf | grep SUBDATASET

Make TIFs from each file in list; replace ‘MOD12Q1:Land_Cover_Type_1’ with desired Subdataset name

mkdir output
find . '*.hdf' -exec gdalwarp -of GTiff 'HDF4_EOS:EOS_GRID:"{}":MOD12Q1:Land_Cover_Type_1' output/{}.tif \;

Merge all .tifs in output directory into single file -o output/Merged_Landcover.tif output/*.tif

BASH functions
Size Functions
This size function echos the pixel dimensions of a given file in the format expected by gdalwarp.

function gdal_size() {
	SIZE=$(gdalinfo $1 |\
		grep 'Size is ' |\
		cut -d\   -f3-4 |\
		sed 's/,//g')
	echo -n "$SIZE"

This can be used to easily resample one raster to the dimensions of another:

gdalwarp -ts $(gdal_size bigraster.tif) -r cubicspline smallraster.tif resampled_smallraster.tif

Extent Functions
These extent functions echo the extent of the given file in the order/format expected by gdal_translate -projwin. (Originally from Linfiniti).

function gdal_extent() {
	if [ -z "$1" ]; then 
		echo "Missing arguments. Syntax:"
		echo "  gdal_extent <input_raster>"
	EXTENT=$(gdalinfo $1 |\
		grep "Upper Left\|Lower Right" |\
		sed "s/Upper Left  //g;s/Lower Right //g;s/).*//g" |\
		tr "\n" " " |\
		sed 's/ *$//g' |\
		tr -d "[(,]")
	echo -n "$EXTENT"

function ogr_extent() {
	if [ -z "$1" ]; then 
		echo "Missing arguments. Syntax:"
		echo "  ogr_extent <input_vector>"
	EXTENT=$(ogrinfo -al -so $1 |\
		grep Extent |\
		sed 's/Extent: //g' |\
		sed 's/(//g' |\
		sed 's/)//g' |\
		sed 's/ - /, /g')
	EXTENT=`echo $EXTENT | awk -F ',' '{print $1 " " $4 " " $3 " " $2}'`
	echo -n "$EXTENT"

function ogr_layer_extent() {
	if [ -z "$2" ]; then 
		echo "Missing arguments. Syntax:"
		echo "  ogr_extent <input_vector> <layer_name>"
	EXTENT=$(ogrinfo -so $1 $2 |\
		grep Extent |\
		sed 's/Extent: //g' |\
		sed 's/(//g' |\
		sed 's/)//g' |\
		sed 's/ - /, /g')
	EXTENT=`echo $EXTENT | awk -F ',' '{print $1 " " $4 " " $3 " " $2}'`
	echo -n "$EXTENT"

Extents can be passed directly into a gdal_translate command like so:

gdal_translate -projwin $(ogr_extent boundingbox.shp) input.tif clipped_output.tif


gdal_translate -projwin $(gdal_extent target_crop.tif) input.tif clipped_output.tif

This can be a useful way to quickly crop one raster to the same extent as another. Add these to your ~/.bash_profile file for easy terminal access.


Top Contributors

dwtkns fitnr veltman timwallace lyzidiamond almccon andrewxhill ermolaev willemarcel