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.

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.

Data conversion

Raster colorspaces conversion

rgb2pct.py

pct2rgb.py

Raster to Raster

Use gdal_translate

Raster to Polygon

gdal_polygonize

Vector to Raster

gdal_rasterize

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.

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

Contour

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:

ls -1 *.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