CODEWITHSUNDEEP
pythongdalgeotiff
Alex
Alex

Reputation: 9

How I can read corner coordinates in degrees?

I have a geotiff, gdalinfo back me info like that: .....

Pixel Size = (0.740750480607563,-0.740750480607565)  Metadata:   AREA_OR_POINT=Area 
Corner Coordinates:  Upper Left  (  408834.576, 4521463.892) ( 13d55'7.11"E,40d50'20.93"N)  
Lower Left (  408834.576, 4497345.057) (13d55'19.75"E, 40d37'18.87"N)

..... I am using python, and read info from tiff like that:

ds = gdal.Open("name.tiff") 
xoff, a, b, yoff, d, e = ds.GetGeoTransform()

xoff, a, b, yoff, d, e will be used in conversion functions(pixel to geo coordinates), but GetGeoTransform() give me xoff and yoff in metrs(this is first number in corner coordinates), not in degrees. How I can read coordinates of Upeer Left corner in degrees? I understand that i can start gdalinfo and read a lot of info, may be in python there is some function for this?

Upvotes: 0

Views: 1816

Answers (1)

Rutger Kassies
Rutger Kassies

Reputation: 64443

You need to convert the coordinates to the other projection, with GDAL this can be done by using the osr module.

It's important to be aware that with since GDAL 3 and above the order of coordinates depends on the projection. The example below shows how you can get the "old" behavior (always x,y), but this of course depends on your usecase. Remove those lines if you don't want that.

Get the properties from the dataset

wkt_srs = ds.GetProjection()
gt = ds.GetGeoTransform()
xs = ds.RasterXSize
ys = ds.RasterYSize

Convert the geotransform to the corner coordinates:

ulx, uly = gdal.ApplyGeoTransform(gt, 0, 0)
lrx, lry = gdal.ApplyGeoTransform(gt, xs, ys)

Set up a transformation object:

src_srs = gdal.osr.SpatialReference()
src_srs.ImportFromWkt(wkt_srs)

tar_srs = gdal.osr.SpatialReference()
tar_srs.ImportFromEPSG(4326)

# with recent versions of GDAL the axis order (x,y vs y,x) depends
# on the projection. Force "x,y" with:
src_srs.SetAxisMappingStrategy(gdal.osr.OAMS_TRADITIONAL_GIS_ORDER)
tar_srs.SetAxisMappingStrategy(gdal.osr.OAMS_TRADITIONAL_GIS_ORDER)

ct = gdal.osr.CoordinateTransformation(src_srs, tar_srs)

Convert the coordinates:

ulx_deg, uly_deg = ct.TransformPoint(ulx, uly)
lrx_deg, lry_deg = ct.TransformPoint(lrx, lry)

Upvotes: 1

Related Questions