要获取地理学家角的坐标, 请按照以下步骤操作:

from osgeo import gdal

ds = gdal.Open('path/to/file')

width = ds.RasterXSize

height = ds.RasterYSize

gt = ds.GetGeoTransform()

minx = gt[0]

miny = gt[3] + width*gt[4] + height*gt[5]

maxx = gt[0] + width*gt[1] + height*gt[2]

maxy = gt[3]

但是, 它们可能不是纬度/经度格式。正如贾斯汀指出的那样, 您的Geotiff将存储在某种坐标系中。如果您不知道坐标系是什么, 可以通过运行

gdalinfo

来找出:

gdalinfo ~/somedir/somefile.tif

输出是什么:

Driver: GTiff/GeoTIFF

Size is 512, 512

Coordinate System is:

PROJCS["NAD27 / UTM zone 11N",

GEOGCS["NAD27",

DATUM["North_American_Datum_1927",

SPHEROID["Clarke 1866",6378206.4,294.978698213901]],

PRIMEM["Greenwich",0],

UNIT["degree",0.0174532925199433]],

PROJECTION["Transverse_Mercator"],

PARAMETER["latitude_of_origin",0],

PARAMETER["central_meridian",-117],

PARAMETER["scale_factor",0.9996],

PARAMETER["false_easting",500000],

PARAMETER["false_northing",0],

UNIT["metre",1]]

Origin = (440720.000000,3751320.000000)

Pixel Size = (60.000000,-60.000000)

Corner Coordinates:

Upper Left ( 440720.000, 3751320.000) (117d38'28.21"W, 33d54'8.47"N)

Lower Left ( 440720.000, 3720600.000) (117d38'20.79"W, 33d37'31.04"N)

Upper Right ( 471440.000, 3751320.000) (117d18'32.07"W, 33d54'13.08"N)

Lower Right ( 471440.000, 3720600.000) (117d18'28.50"W, 33d37'35.61"N)

Center ( 456080.000, 3735960.000) (117d28'27.39"W, 33d45'52.46"N)

Band 1 Block=512x16 Type=Byte, ColorInterp=Gray

此输出可能是您需要的。但是, 如果要使用Python进行编程操作, 将获得相同的信息。

如上例所示, 如果坐标系为

PROJCS

, 则您正在处理投影坐标系。投影坐标系是球形地球表面的

表示形式, 但它是平坦的并且在

平面上变形。如果需要纬度和经度, 则需要将坐标转换为正确的地理坐标系。

不幸的是, 并非所有纬度/经度对都是相同的, 因为它们基于不同的球形地球模型。在此示例中, 我将转换为

WGS84

, 这是GPS中首选的地理坐标系, 并且被所有流行的地图站点使用。坐标系由定义明确的字符串定义。该目录可以从

空间参考

获得, 例如, 请参见

WGS84

。

from osgeo import osr, gdal

# get the existing coordinate system

ds = gdal.Open('path/to/file')

old_cs= osr.SpatialReference()

old_cs.ImportFromWkt(ds.GetProjectionRef())

# create the new coordinate system

wgs84_wkt = """

GEOGCS["WGS 84",

DATUM["WGS_1984",

SPHEROID["WGS 84",6378137,298.257223563,

AUTHORITY["EPSG","7030"]],

AUTHORITY["EPSG","6326"]],

PRIMEM["Greenwich",0,

AUTHORITY["EPSG","8901"]],

UNIT["degree",0.01745329251994328,

AUTHORITY["EPSG","9122"]],

AUTHORITY["EPSG","4326"]]"""

new_cs = osr.SpatialReference()

new_cs .ImportFromWkt(wgs84_wkt)

# create a transform object to convert between coordinate systems

transform = osr.CoordinateTransformation(old_cs,new_cs)

#get the point to transform, pixel (0,0) in this case

width = ds.RasterXSize

height = ds.RasterYSize

gt = ds.GetGeoTransform()

minx = gt[0]

miny = gt[3] + width*gt[4] + height*gt[5]

#get the coordinates in lat long

latlong = transform.TransformPoint(minx,miny)

希望这可以满足您的需求。