从GeoTIFF文件获取纬度和经度

Question

在Python中使用GDAL，怎么才能获取一个GeoTIFF文件的经纬度呢？

GeoTIFF文件似乎并不直接存储坐标信息。相反，它们保存的是XY原点坐标。不过，这些XY坐标并不能告诉我们左上角和左下角的经纬度。

看起来我需要做一些数学计算来解决这个问题，但我不知道该从哪里开始。

需要什么步骤才能完成这个操作呢？

我知道GetGeoTransform()这个方法对这个问题很重要，但我不知道接下来该怎么用它。

Answer 1

我不知道这算不算完整的回答，但这个网站上说：

x/y地图的维度被称为“东偏移”和“北偏移”。对于地理坐标系统的数据集，这些值代表经度和纬度。而对于投影坐标系统，它们通常是投影坐标系统中的东偏移和北偏移。对于没有地理参考的图像，东偏移和北偏移只是每个像素的像素/行偏移（这可以通过单位地理变换来理解）。

所以它们实际上可能就是经度和纬度。

Answer 2

要获取你的地理信息图（geotiff）角落的坐标，可以按照以下步骤进行：

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] 

不过，这些坐标可能不是以经纬度的格式呈现。正如Justin提到的，你的地理信息图会使用某种坐标系统存储。如果你不知道具体是什么坐标系统，可以通过运行 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常用的地理坐标系统，也是所有流行的网络地图网站使用的系统。这个坐标系统是通过一个明确的字符串来定义的。你可以在 spatial ref 找到它们的目录，比如 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) 

希望这能满足你的需求。