2

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=2145x1 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()
  • Please list what you did in python, and what you got. – mdsumner Nov 5 '15 at 23:41
  • @mdsummer python code and result has been added. – apricity Nov 6 '15 at 0:21
  • new_cs should be the string for long-lat, here it's the same as the source – mdsumner Nov 6 '15 at 6:22
4

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?

  • I think the link you pointed me to is probably the solution, but I don't know what my source EPSG is. GDAL info says AUTHORITY["EPSG","9001"] but when I try going from 9001 to 4326 I get the error: TypeError: in method 'Geometry_Transform', argument 2 of type 'OSRCoordinateTransformationShadow *' – apricity Nov 6 '15 at 18:56
  • 9001 is the EPSG code for a metre (i.e. that is part of the UNIT), not a projected coordinate system. – BradHards Nov 6 '15 at 23:17
  • Your source CRS does not have an EPSG number (and does not have to, because you can build projections with almost any parameter set you want). – AndreJ Nov 7 '15 at 9:42
0

Thanks for the pointers all. I now understand that even though my source geotif has no EPSG I can get its projection from wkt and view it as proj.4:

ds = gdal.Open("mygeotiff")
in_srs = osr.SpatialReference()
in_srs.ImportFromWkt(ds.GetProjection())
print(in_srs.ExportToProj4()) : '+proj=lcc +lat_1=25 +lat_2=25 +lat_0=25 +lon_0=265 +x_0=0 +y_0=0 +a=6371229 +b=6371229 +units=m +no_defs '

The below code translates coordinates to long, lat in 4326 correctly without having a source EPSG:

# translates coord_x, coord_y from in_spatial_ref to coords in output_epsg
def get_point_as_long_lat(coord_x, coord_y, in_spatial_ref, output_epsg):
  # target Spatial Reference System
  out_spatial_ref = osr.SpatialReference()
  out_spatial_ref.ImportFromEPSG(output_epsg)

  coord_transform = osr.CoordinateTransformation(in_spatial_ref, out_spatial_ref)

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

  # transform point
  point.Transform(coord_transform)

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


ds = gdal.Open("mygeotiff")
gt = ds.GetGeoTransform()

in_srs = osr.SpatialReference()
in_srs.ImportFromWkt(ds.GetProjection())

origin_long_lat = get_point_as_long_lat(geo_transform[0], geo_transform[3], in_srs, 4326)

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.