[GIS] Getting origin as long, lat in GDAL without knowing EPSG

gdallatitude longitude

My goal is to get the long, lat coordinates for the origin of a particular geotif.

gdalinfo outputs:
Upper Left (-2764486.928, 3232110.510) (130d 7'22.42"W, 49d56'48.81"N)
But I don't understand the meaning of (-2764486.928, 3232110.510) – for instance what units are being represented here?

What I want is something like:
Upper Left (-125.0208333, 49.9375000) (125d 1'15.00"W, 49d56'15.00"N)
which gdalinfo yields on a different geotiff file that I have.

I'm actually using ds.GetGeoTransform() in python to get this same information gdalinfo outputs and am looking for a solution in python to convert (-2764486.928, 3232110.510) into a long, lat. I tried fmark's answer here, but it didn't convert (-2764486.928, 3232110.510) into a long,lat.

Here is the entire gdalinfo output

Driver: GTiff/GeoTIFF
Files: mask.tif
Size is 2145, 1377
Coordinate System is:
PROJCS["unnamed",
GEOGCS["Coordinate System imported from GRIB file",
DATUM["unknown",
SPHEROID["Sphere",6371229,0]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433]],
PROJECTION["Lambert_Conformal_Conic_2SP"],
PARAMETER["standard_parallel_1",25],
PARAMETER["standard_parallel_2",25],
PARAMETER["latitude_of_origin",25],
PARAMETER["central_meridian",265],
PARAMETER["false_easting",0],
PARAMETER["false_northing",0],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]]]
Origin = (-2764486.928100540300000,3232110.510093217700000)
Pixel Size = (2539.703000000000000,-2539.703000000000000)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
INTERLEAVE=BAND
Corner Coordinates:
Upper Left (-2764486.928, 3232110.510) (130d 7'22.42"W, 49d56'48.81"N)
Lower Left (-2764486.928, -265060.521) (121d33'48.72"W, 20d10'43.04"N)
Upper Right ( 2683176.007, 3232110.510) ( 60d51'59.22"W, 50d 6'46.25"N)
Lower Right ( 2683176.007, -265060.521) ( 69d11'55.57"W, 20d19' 6.55"N)
Center ( -40655.461, 1483524.995) ( 95d27' 9.15"W, 38d13' 5.62"N)
Band 1 Block=2145×1 Type=Float32, ColorInterp=Gray

Here's what I tried in python:

geo_transform = ds.GetGeoTransform()
old_cs= osr.SpatialReference()
old_cs.ImportFromWkt(ds.GetProjectionRef())
wkt = """
    GEOGCS["NAD83",
    DATUM["North_American_Datum_1983",
        SPHEROID["GRS 1980",6378137,298.2572221010002,
            AUTHORITY["EPSG","7019"]],
        TOWGS84[0,0,0,0,0,0,0],
        AUTHORITY["EPSG","6269"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4269"]]"""
new_cs = osr.SpatialReference()
new_cs.ImportFromWkt(wkt)
transform = osr.CoordinateTransformation(old_cs,new_cs)
width = ds.RasterXSize
height = ds.RasterYSize
gt = ds.GetGeoTransform()

#get the coordinates in lat long
latlong = transform.TransformPoint(gt[0], gt[3])

The output is latlong = (-2764486.928, 3232110.510, 0.0)

Using BradHards's link, I've got the code pasted below and believe I just need to figure out the source geotiffs EPSG. gdalinfo says it's 9001, but when I try using that the code blows up with 

TypeError: in method 'Geometry_Transform', argument 2 of type 'OSRCoordinateTransformationShadow *'

The geotiff originated as a GRIB file from nomads

ds = gdal.Open("mygeotiff")
gt = ds.GetGeoTransform()
long_lat = get_point_as_long_lat(gt[0], gt[3], ???what_epsg???)


def get_point_as_long_lat(coord_x, coord_y, input_epsg):
  # Spatial Reference System
  output_epsg = 4326

  # create a geometry from coordinates
  point = ogr.Geometry(ogr.wkbPoint)
  point.AddPoint(coord_x, coord_y)

  # create coordinate transformation
  in_spatial_ref = osr.SpatialReference()
  in_spatial_ref.ImportFromEPSG(input_epsg)

  out_spatial_ref = osr.SpatialReference()
  out_spatial_ref.ImportFromEPSG(output_epsg)

  coord_transform = osr.CoordinateTransformation(in_spatial_ref, out_spatial_ref)

  # transform point
  point.Transform(coord_transform)

  # return point in EPSG 4326
  return point.GetX(), point.GetY()

Best Answer

The coordinates are in the spatial reference system described in the file.

So -2764486.928, 3232110.510 is referenced to something (that gdal doesn't have a name for). The something is:

PROJCS["unnamed",
GEOGCS["Coordinate System imported from GRIB file",
DATUM["unknown",
SPHEROID["Sphere",6371229,0]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433]],
PROJECTION["Lambert_Conformal_Conic_2SP"],
PARAMETER["standard_parallel_1",25],
PARAMETER["standard_parallel_2",25],
PARAMETER["latitude_of_origin",25],
PARAMETER["central_meridian",265],
PARAMETER["false_easting",0],
PARAMETER["false_northing",0],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]]]

gdalinfo does show it in something more familiar - the part following: (130d 7'22.42"W, 49d56'48.81"N). You should be able to get that with python too. You haven't posted all of the code (e.g. how did you open the dataset), but there are plenty of examples, including How to convert projected coordinates to lat/lon using Python?

Related Question