从GeoTIFF文件获取纬度和经度

Question

在Python中使用GDAL,如何获得GeoTIFF文件的纬度和经度？

GeoTIFF似乎不存储任何坐标信息.相反,它们存储XY原点坐标.但是,XY坐标不提供左上角和左下角的纬度和经度.

看来我需要做一些数学来解决这个问题,但我不知道从哪里开始.

执行此操作需要什么程序？

我知道这个GetGeoTransform()方法对此很重要,但是,我不知道该怎么做.

Answer 1

要获取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所说,你的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所支持的地理坐标系统,并被所有流行的Web地图站点使用.坐标系由定义良好的字符串定义.它们的目录可从空间参考中获得,例如参见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) 

希望这会做你想要的.

Answer 2

我不知道这是否是一个完整的答案,但这个网站说:

x/y地图尺寸称为东向和北向.对于地理坐标系中的数据集,这些数据集将保持经度和纬度.对于投影坐标系,它们通常是投影坐标系中的东向和北向.对于未经过引用的图像,东向和北向只是每个像素的像素/线偏移(如统一地理变换所暗示的).

所以他们实际上可能是经度和纬度.