GDAL and OGR: geodata conversion and re-projection tools

Overview

Operations on geographic data are most efficient when the input files have identical spatial parameters: i.e. coordinates system and datum / ellispoid. The GDAL/OGR tools and libraries can be used to "reproject" downloaded geographic data from one coordinate system to another, both for vector and raster data. More information on this software can be found here: http://www.remotesensing.org/gdal/index.html The following examples illustrate some of the ways in which the GDAL/OGR tools can be used to manipulate the spatial parameters of some common types of geographic data: SSURGO (ESRI Shapefile format) and aerial imagery (GeoTiff format).

How Does it Compare with ArcToolbox?

A SSURGO boundary polygon was used in a quick test to see how similar the results of a vector re-projection done by ArcMap and OGR were. The input vector (GCS NAD83) was projected to UTM zone 11 NAD83 by both OGR and ArcMap 9.0, and then imported into GRASS. Both vectors were then exported as ascii, for simple access to the line vertex coordinates. A simple AWK script was used to numerically compare the differences in both X and Y coordinates:

paste arcmap.ascii ogr.ascii | \ awk ' /\./{print ($1-$3), ($2-$4) ; x+=($1-$3) ; y+=($2-$4); n++ } END{print "avg dx: "(x/n), "avg dy: "(y/n), "n: "n} '

The resulting average dx and dy values were: dx: 2.98838e-06 m and dy: -1.23948e-07 m with 5035 samples: looks good to me. Next time, a datum shift may be more telling of any subtle differences to be aware of. Update: NAD83 -> NAD27 Datum shift experiment looks like datum shifts are comparable as well.

Examples:
Command names are in boldface type, input files are red, and output files are in green. Figures are provided for additional context, and are referenced by input/output file name in the example code.

Vector Operations

Get information about a shapefile:

ogrinfo -al ssurgo_geo.shp | less

Reproject SSURGO mapunit data from geographic coordinates (Lat/Lon) to UTM zone 10:

ogr2ogr -t_srs "+proj=utm +zone=10 +datum=NAD83" ssurgo_utm.shp ssurgo_geo.shp

Extract all polygons from a SSURGO shapefile where the mapunit symbol is 'ScA' :

ogr2ogr -where "musym = 'ScA' " ssurgo_ScA.shp ssurgo_utm.shp

Convert a shapefile into Google-friendly KML:

ogr2ogr -f "KML" -t_srs EPSG:4326 trails.kml trails.shp

Raster Operations [raster tutorial]

Get information about a raster dataset

gdalinfo rasterfile.tiff

Reproject an aerial photo in CA State Plane Zone 4 (Lambert Conformal Conic projection, units = feet) to UTM Zone 10 (units = meters), and rescale to 1 meter output resolution, and use thin-plate-spline resampling (-tps):

gdalwarp -tps -t_srs '+proj=utm +zone=10 +datum=NAD83 +units=m' \
-s_srs '+proj=lcc +lat_1=36.0000 +lat_2=37.2500 +lat_0=35.333333333333336 \
+lon_0=-119.0000 +x_0=2000000 +y_0=500000 +datum=NAD83 +units=ft' \
-srcnodata 255 -dstnodata 255 -tr 1 1 884084.tif 884084-utm.tif

Note that it is usually a good idea to "optimise" the resulting image with gdal_translate.
Convert Multi-band GeoTiff file to JPEG:

gdal_translate -of JPEG 884084-utm.tif 884084-utm.jpg

 Reading MrSid format files details

Add the DSDK path to configure: --with-mrsid=/Geo_DSDK-6.0.7.1407/
There is a bug in Geo_DSDK-6.0.7.1407//include/base/lti_sceneBuffer.h line 356:

bool LTISceneBuffer::inWindow(lt_uint32 x, lt_uint32 y) const; //

Comment out this line, and re-compile. Thanks to Mateusz Loskot on the GDAL IRC channel for pointing this out.

gdal_translate file.sid file.tif
ssurgo_utm.preview.png
ssurgo_utm.shp
ssurgo_ScA.preview.png
ssurgo_ScA.shp
884084-utm.jpg