geotif 添加坐标_从GeoTIFF文件获取经度和纬度

使用Python中的GDAL库,通过GetGeoTransform()方法获取GeoTIFF文件的角落坐标。但这些坐标可能不是经纬度格式,需结合gdalinfo找出坐标系统,并进行转换,例如从NAD27转换到WGS84以得到正确的经纬度。
摘要由CSDN通过智能技术生成

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.

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值