GDAL and OGR: geodata conversion and re-projection tools

Submitted by dylan on Mon, 2005-09-12 21:53.

 
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 comparible 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.

ssurgo utm example
ssurgo_utm.shp
ssurgo ScA
ssurgo_ScA.shp

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

image CA SPC IV

884084.tif

image utm zone 10

884084-utm.tif



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; // <--- line 356

     
    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






( categories: )