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

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

要获取您的Geotiff角的坐标，请执行以下 *** 作：

from osgeo import gdalds = gdal.Open('path/to/file')width = ds.RasterXSizeheight = ds.RasterYSizegt = 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]

但是，它们可能不是纬度/经度格式。正如Justin指出的那样，您的Geotiff将使用某种坐标系进行存储。如果您不知道它是什么坐标系，可以通过运行

gdalinfo

以下命令进行查找：

gdalinfo ~/somedir/somefile.tif

哪个输出：

Driver: GTiff/GeoTIFFSize is 512, 512Coordinate 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中以编程方式执行此 *** 作，则可以通过此方法获取相同的信息。

如果坐标系

PROJCS

类似于上面的示例，则说明您正在处理投影坐标系。投影的Coordiante系统是球形地球表面的表示，但被展平并扭曲到一个平面上。如果需要纬度和经度，则需要将坐标转换为所需的地理坐标系。

遗憾的是，并非所有的纬度/经度对都是基于地球的不同球体模型而创建的。在此示例中，我将转换为WGS84，这是GPS偏爱的地理坐标系，并且所有流行的Web制图站点都使用该地理坐标系。坐标系由定义明确的字符串定义。它们的目录可从空间参考中获得，例如参见WGS84。

from osgeo import osr, gdal# get the existing coordinate systemds = gdal.Open('path/to/file')old_cs= osr.SpatialReference()old_cs.importFromWkt(ds.GetProjectionRef())# create the new coordinate systemwgs84_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 systemstransform = osr.CoordinateTransformation(old_cs,new_cs)#get the point to transform, pixel (0,0) in this casewidth = ds.RasterXSizeheight = ds.RasterYSizegt = ds.GetGeoTransform()minx = gt[0]miny = gt[3] + width*gt[4] + height*gt[5]#get the coordinates in lat longlatlong = transform.TransformPoint(minx,miny)

希望这会做您想要的。