从GeoTIFF文件获取经度和纬度 [英] Obtain Latitude and Longitude from a GeoTIFF File

问题描述

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

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

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

需要执行哪些操作？

我知道 GetGeoTransform（）方法对此很重要，但是我不知道如何处理它。

解决方案

要获得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] 
  

但是，这些可能不是纬度/经度格式。正如贾斯汀所说，你的geotiff将以某种坐标系统存储。如果你不知道它是什么坐标系，你可以运行 gdalinfo ：

  gdalinfo〜/ somedir / somefile.tif 
  

输出：

 驱动程序：GTiff / GeoTIFF 
尺寸为512,512 
坐标系为：
 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]，
单位[米，1]] 
原点=（440720.000000,3751320.000000）
像素大小=（60.000000，-60.000000）
角坐标：
左上角（440720.000 ， （117天38'28.21W，33天54'8.47N）
左下（440720.000，3720600.000）（117d38'20.79W，33d37'31.04N）
右上（471440.000，3751320.000） （117d18'32.07W，33d54'13.08N）
右下（471440.000,3720600.000）（117d18'28.50W，33d37'35.61N）
中心（456080.000,3735960.000）（117d28' 27.39W，33d45'52.46N）
 Band 1 Block = 512x16 Type = Byte，ColorInterp = Gray 
  

这个输出可能是你所需要的。如果你想在Python中以编程方式执行此操作，则可以通过这种方式获得相同的信息。

如果坐标系是 PROJCS 就像上面你正在处理投影坐标系的例子。一个预计的坐标系统是一个球面地球表面的表示，但在飞机上变平和变形。如果您想要经纬度，您需要将坐标转换为您想要的地理坐标系。

可悲的是，并非所有纬度/经度对都相同，基于不同的地球球体模型。在这个例子中，我将转换为 WGS84 ，这是地理坐标系统在全球定位系统中最受青睐并被所有人流行的网站映射网站。坐标系由定义良好的字符串定义。它们的目录可从空间参考中获得，请参阅 WGS84 。

  from osgeo import osr，gdal 
 
＃获取现有坐标系
 ds = gdal.Open（'path / to / file'）
 old_cs = osr.SpatialReference（）
 old_cs.ImportFromWkt（ ds.GetProjectionRef（））
 
＃创建新坐标系
 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）
 
 ＃创建一个转换对象在坐标系之间进行转换
 transform = osr.CoordinateTransformation（old_cs，new_cs）
 
＃获取要转换的点，像素（0,0）在这种情况下为
 width = ds。 RasterXSize 
 height = ds.RasterYSize 
 gt = ds.GetGeoTransform（）
 minx = gt [0] 
 miny = gt [3] + width * gt [4] + height * gt [5] 
 
＃获取纬度坐标
 latlong = transform.TransformPoint（x，y）
  

希望这会做你想做的。

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

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?

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

解决方案

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] 

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 

Which outputs:

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

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.

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.

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(x,y) 

Hopefully this will do what you want.

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