3

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)

| improve this question | |
  • Are you programming or just trying to work out where on the image a coordinate is to go to it? I have already answered this question before for programming, if that is what you're after have a look here gis.stackexchange.com/questions/96050/… but it's important that your query point and your raster are in the same coordinate system if they're not that's a different matter but can be solved with proj lib that comes with GDAL. – Michael Stimson Oct 19 '15 at 23:17
  • I need to make a topo map. I have successfully cropped out and downsampled a region from a giant img file with gdal_translate -of xyz -srcwin 12000 11000 7000 2000 -outsize 22.49% 22.49% topo.img tmp.txt. If I can get the corner coordinates of what I have cropped out, then I know how to make a topo map with GMT. In the example you link, it is assumed that x and longitude are related linearly, which is not exactly true if I look at the corner coordinates of my rectangular grid. – Norbert Oct 19 '15 at 23:32
  • Are your image and query coordinates in a different spatial reference (for example you have Geographic lat/lon and your image is in metres UTM)? If you can tell me what programming/scripting languages you have available I can walk you through projecting coordinates otherwise use GDALwarp to transform your image from one spatial reference to another such that query and image are the same. – Michael Stimson Oct 19 '15 at 23:47
  • I am doing this on a Linux prompt, but could as well write a shell script. I added information about my input file above. I had also tried img2grd, but can't get it to produce any output. Thanks! – Norbert Oct 20 '15 at 0:02
  • 1
    Why don't you just crop out to something that will maintain the metadata (e.g. GeoTiff), then you can just get it from gdalinfo output_file.tif – BradHards Oct 20 '15 at 2:29
3

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.

| improve this answer | |
  • thank god lol...pyproj is apparently a nightmare and this worked just great – Derek Eden Feb 19 at 3:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.