Advertisement

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

阅读量:

要获取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] + widthgt[4] + heightgt[5]

maxx = gt[0] + widthgt[1] + heightgt[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编程实现这一点,这就是获得相同信息的方式。

可悲的是,并非所有的经纬度对都是基于地球不同的球体模型而产生的。在本例中,我将转换为WGS84,这是GPSs中最受欢迎的地理坐标系,所有流行的web地图站点都使用它。坐标系由定义良好的字符串定义。它们的目录可从spatial ref获得,请参见例如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] + widthgt[4] + heightgt[5]

#get the coordinates in lat long

latlong = transform.TransformPoint(minx,miny)

希望这能满足你的要求。

全部评论 (0)

还没有任何评论哟~