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:
Reproject SSURGO mapunit data from geographic coordinates (Lat/Lon) to UTM zone 10:
Extract all polygons from a SSURGO shapefile where the mapunit symbol is 'ScA' :
Convert a shapefile into Google-friendly KML:
Raster Operations [raster tutorial]
Get information about a raster dataset
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):
-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:
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:
Comment out this line, and re-compile. Thanks to Mateusz Loskot on the GDAL IRC channel for pointing this out.
 ssurgo_utm.shp |
 ssurgo_ScA.shp |
 |
