当前位置:   article > 正文

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

根据经纬度坐标使用geotools找到tiff文件中对应的点

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.

声明:本文内容由网友自发贡献,不代表【wpsshop博客】立场,版权归原作者所有,本站不承担相应法律责任。如您发现有侵权的内容,请联系我们。转载请注明出处:https://www.wpsshop.cn/w/Gausst松鼠会/article/detail/329908
推荐阅读
相关标签
  

闽ICP备14008679号