GDAL reference

General notes

It can be helpful to run gdalwarp with the commandline option "--debug on" in which case the debug output should include the PROJ.4 rendering of the source and destination coordinate system.

At a very low level you can also define the PROJ_DEBUG environment variable to the value ON and PROJ.4 will report some details on what datum files are opened.

GDAL OGR Cookbook

A very nice resource: https://pcjericks.github.io/py-gdalogr-cookbook/

Clip or Extract

Geometric

If you want to clip a rectangle out of a larger raster, for example to remove the collar from a scanned map, you can use gdal_translate but only if the input raster is not in a spherical projection. (For example, if it's in WGS84 (EPSG 4326) this will work.)

If the data is not projected you can use a gdal_translate command with either projwin or srcwin specifying a window on the command line.

If it is projected it will fail.

Using a mask

This tutorial shows how to use a shapefile to clip a raster: http://linfiniti.com/2009/09/clipping-rasters-with-gdal-using-polygons/

I want to use the neatline from a GeoPDF to remove the collar from a USGS scanned map. To do this, I need to make a polygon in a shapefile from the points in the neatline.

Distance

gdal_proximity

Generalization

gdal_sieve

Raster polygons smaller than a provided threshold size (in pixels) and replaces replaces them with the pixel value of the largest neighbour polygon. The result can be written back to the existing raster band, or copied into a new file.

Map algebra

Data conversion

Raster colorspaces conversion

rgb2pct.py

pct2rgb.py

X,Y Table to Vector File

Import data from GPX file to File Geodatabase

This has to be better than using the standard ESRI tool! First, install the FileGDB support. See Building GDAL on a Mac for example.

Copy the file from the Garmin then run it through the filter.

EPSG:4152 is HARN Cal State III in feet

cp /Volumes/GARMIN/Garmin/GPX/Current/Current.gpx .
GPXNAME=`date +"gpx_%Y-%m-%d"`
PROJ="EPSG:4152"
FGDB="ws_gpx.gdb"
ogr2ogr -f FileGDB ${GPXNAME}_track_points Current.gpx tracks

Other useful options: -lco "FEATURE_DATASET=foo" would create a feature dataset called "foo" and put the feature class into it.

Import PDF to TIFF

Split the pages in the PDF apart

pdfseparate NAME.pdf "NAME_%d.pdf"

Convert the PDFs into TIFFs

If your version of gdal cannot handle PDF files you can use pdftoppm and ppm2tiff or use pdftoppm and gdal_translate.

for i in NAME_*.pdf; do
b=`basename $i '.pdf'`
# gdal_translate -of GTiff -co "TILED=YES" -co "TFW=YES" $i $b.tif
# or perhaps
pdftoppm "$i" | ppm2tiff "$b.tif"
done

Now you can georeference the TIFF file(s) if the PDFs were not GeoPDFs.

Data management

Projections

To reproject data, use gdalwarp

I had trouble going from NAD83HARN to NAD27. Using --debug on reveals no datum is defined, even though it is in the SID.

OGRCT: Source: +proj=lcc +lat_1=44.33333333333334
+lat_2=46 +lat_0=43.66666666666666 +lon_0=-120.5
+x_0=2500000 +y_0=0
+ellps=GRS80
+units=m +no_defs

I can override the source SRS on the command line, which gives me same results when I define it as HARN. When I override it and leave out the HARN, it finds the datum correctly. This is a bug.

OGRCT: Source: +proj=lcc +lat_1=44.33333333333334
+lat_2=46 +lat_0=43.66666666666666 +lon_0=-120.5
+x_0=2500000 +y_0=0
+datum=NAD83
+units=m +no_defs

Still no go, so I try referring to http://trac.osgeo.org/proj/wiki/FAQ

Mosaic

The gdal_merge command creates a single raster output from several inputs. It does this the old-fashioned way by copying pixels.

Build raster catalog

(Or, "what is a VRT")

In GDAL world this is called a VRT = Virtual Dataset (which in my mind would abbreviate to "VD" but no one asked me.)

To build one you can use gdalbuildvrt

I think it creates a small file referencing the members of the dataset, which can then be processed by GDAL commands that accept only one input file.

Example: takes a couple seconds to create the output file, hwy17.vrt

cd ~/GISData/CA/SantaClara/USGS_Ortho/hwy17/
gdabuildvrt hwy17.vrt *.jp2

and the resultant vrt file can be used in ArcGIS.

Surfaces

Contours

To build a vector contour layer from an elevation dataset, use gdal_contour

I got data for Vallejo area and built a 1' contour shapefile using this command:

find . -name *.tif > list.txt
gdalbuildvrt -input_file_list list.txt napa.vrt
gdal_contour -i .305 -a elevation napa.vrt napa_contour_ft.shp

The input TIFF is in meters so .305 gives an interval of 1 foot, but comes up with a messy elevation attribute that I recalculate in ArcMap.

ArcGIS 10.4 was also able to directly load the .vrt file, something that really surprised me. It could not display it in ArcMap but I was able to use Spatial Analyst to build a contour layer too, for comparison.

Raster surfaces

To generate raster surfaces from an elevation dataset, use gdaldem

The list of surfaces you can generate includes

hillshade
color relief
shaded relief
slope
aspect
terrain ruggedness index
topographic position index
roughness

GDAL reference

Contents