[GIS] Simple gdal coordinate conversion

gdal

Starting with an .img file how to I convert a specific (x,y) coordinate into (lon,lat) with gdal. I don't know what my input format is, but if I type gdalinfo topo.img it tells me the corner coordinates of the map, so it ought to be possible to get the coordinates for any of the grid points without further input. gdallocationinfo topo.img [x] [y] gives the elevation only.

The beginning of the gdalinfo output for the input file looks like this:

Driver: HFA/Erdas Imagine Images (.img)
Files: bigislandDTM.img
bigislandDTM.ige
bigislandDTM.rrd
Size is 32609, 35573
Coordinate System is:
PROJCS["NAD_1983_UTM_Zone_5N",
GEOGCS["NAD83",
DATUM["North_American_Datum_1983",
SPHEROID["GRS 1980",6378137,298.257222101,
AUTHORITY["EPSG","7019"]],
TOWGS84[0,0,0,0,0,0,0],
AUTHORITY["EPSG","6269"]],
PRIMEM["Greenwich",0,
AUTHORITY["EPSG","8901"]],
UNIT["degree",0.0174532925199433,
AUTHORITY["EPSG","9108"]],
AUTHORITY["EPSG","4269"]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",-153],
PARAMETER["scale_factor",0.9996],
PARAMETER["false_easting",500000],
PARAMETER["false_northing",0],
UNIT["Meter",1],
AUTHORITY["EPSG","26905"]]
Origin = (170717.059624216170050,2247879.339183096773922)
Pixel Size = (4.497066284309543,-4.497066284309542)
Corner Coordinates:
Upper Left ( 170717.060, 2247879.339) (156d 9' 9.63"W, 20d18' 1.76"N)
Lower Left ( 170717.060, 2087905.200) (156d 7'28.87"W, 18d51'25.08"N)
Upper Right ( 317361.894, 2247879.339) (154d44'57.59"W, 20d19'12.58"N)
Lower Right ( 317361.894, 2087905.200) (154d44' 1.61"W, 18d52'30.47"N)
Center ( 244039.477, 2167892.270) (155d26'24.29"W, 19d35'22.45"N)

Best Answer

Here is a very simple python script using GDAL (well, OGR and OSR which are part of the GDAL package) to project a coordinate from one spatial reference to another:

from osgeo import ogr
from osgeo import osr

InSR = osr.SpatialReference()
InSR.ImportFromEPSG(4326)       # WGS84/Geographic
OutSR = osr.SpatialReference()
OutSR.ImportFromEPSG(32756)     # WGS84 UTM Zone 56 South

Point = ogr.Geometry(ogr.wkbPoint)
Point.AddPoint(150.700532,-35.145618) # use your coordinates here
Point.AssignSpatialReference(InSR)    # tell the point what coordinates it's in
Point.TransformTo(OutSR)              # project it to the out spatial reference
print '{0},{1}'.format(Point.GetX(),Point.GetY()) # output projected X and Y coordinates

output: 290523.025428,6108387.72056

The values here are hard coded to show an example of how to convert from one to another but it shouldn't be too difficult to modify to suit your needs. You already have GDAL installed and probably have python.. if not GDAL/OGR/OSR is available in many more languages but is fundamentally the same process.

Related Question