从 GeoTIFF 文件中获取纬度和经度

问题描述

在Python中使用GDAL，如何获取GeoTIFF文件的经纬度?

Using GDAL in Python, how do you get the latitude and longitude of a GeoTIFF file?

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

GeoTIFF's do not appear to store any coordinate information. Instead, they store the XY Origin coordinates. However, the XY coordinates do not provide the latitude and longitude of the top left corner and bottom left corner.

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

It appears I will need to do some math to solve this problem, but I don't have a clue on where to start.

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

What procedure is required to have this performed?

我知道 GetGeoTransform() 方法对此很重要，但是，我不知道如何处理它.

I know that the GetGeoTransform() method is important for this, however, I don't know what to do with it from there.

推荐答案

要获取 geotiff 角的坐标，请执行以下操作:

To get the coordinates of the corners of your geotiff do the following:

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来查找:

However, these might not be in latitude/longitude format. As Justin noted, your geotiff will be stored with some kind of coordinate system. If you don't know what coordinate system it is, you can find out by running 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 中以编程方式执行此操作，这就是您获得相同信息的方式.

This output may be all you need. If you want to do this programmaticly in python however, this is how you get the same info.

如果坐标系是 PROJCS 像上面的例子，你正在处理一个投影坐标系.投影坐标系是球状地球表面的表示，但在平面上变平并变形.如果你想要经纬度，你需要将坐标转换成你想要的地理坐标系.

If the coordinate system is a PROJCS like the example above you are dealing with a projected coordinate system. A projected coordiante system is a representation of the spheroidal earth's surface, but flattened and distorted onto a plane. If you want the latitude and longitude, you need to convert the coordinates to the geographic coordinate system that you want.

遗憾的是，并非所有纬度/经度对都是平等的，这是基于地球的不同球体模型.在此示例中，我将转换为 WGS84，这是 GPS 中偏爱并被所有人使用的地理坐标系流行的网络地图网站.坐标系由定义明确的字符串定义.它们的目录可从 spatial ref 获得，例如参见 WGS84.

Sadly, not all latitude/longitude pairs are created equal, being based upon different spheroidal models of the earth. In this example, I am converting to WGS84, the geographic coordinate system favoured in GPSs and used by all the popular web mapping sites. The coordinate system is defined by a well defined string. A catalogue of them is available from spatial ref, see for example 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) 

希望这能满足您的需求.

Hopefully this will do what you want.

这篇关于从 GeoTIFF 文件中获取纬度和经度的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持html5模板网！