简介

参考文献¹

GDAL 是一个开源的操作栅格数据和矢量数据的库，本文记录下用 Python 中 GDAL 库操作 TIFF （GeoTIFF）的常见代码，包括读写、获取坐标系、获取指定位置像元值等。
TIFF 简单理解就是一种图像格式，类似于 jpg、png 等。
GeoTIFF 就是在普通 TIFF 文件上增加了地理位置、投影信息、坐标信息等，常用于遥感数据。

一、读取TIFF文件

from osgeo import gdal, osrfilename=r'F:\GOSIF_2001065.tif'
tiff_path=r'F:\GOSIF\8day'def GetTifInfo():dataset = gdal.Open(filename)  # 打开文件im_width = dataset.RasterXSize  # 栅格矩阵的列数im_height = dataset.RasterYSize  # 栅格矩阵的行数im_bands = dataset.RasterCount  # 波段数# print(im_bands)#栅格数据的六参数。# geoTransform[0]：左上角像素经度# geoTransform[1]：影像宽度方向上的分辨率(经度范围/像素个数)# geoTransform[2]：x像素旋转, 0表示上面为北方# geoTransform[3]：左上角像素纬度# geoTransform[4]：y像素旋转, 0表示上面为北方# geoTransform[5]：影像宽度方向上的分辨率(纬度范围/像素个数)extent = dataset.GetGeoTransform()  # 仿射矩阵，左上角像素的大地坐标和像素分辨率var = dataset.GetProjection()   # 栅格数据的投影# osr.SpatialReference 提供描绘和转换坐标系统的服务 地理坐标(用经纬度表示)；投影坐标(如 UTM ，用米等单位量度来定位)。pcs = osr.SpatialReference()# ImportFromWkt()函数可以把 WKT坐标系统设置到OGRSpatialReference中pcs.ImportFromWkt(dataset.GetProjection())gcs = pcs.CloneGeogCS()  #地理空间坐标系print(gcs)shape = (dataset.RasterYSize, dataset.RasterXSize)img = dataset.GetRasterBand(1).ReadAsArray()  # (height, width)# print(img)
# img(ndarray), gdal数据集、地理空间坐标系、投影坐标系、栅格影像大小return img, dataset, gcs, pcs, extent, shape

二、经纬度、行列号、投影坐标之间的转换

1.经纬度转投影坐标xy

参考文献²

地理坐标系（Geographic coordinate system），Geographic coordinate system直译为地理坐标系统，是以经纬度为地图的存储单位的。很明显，Geographic coordinate system是球面坐标系统。
Projection coordinate system（投影坐标系统），每一个投影坐标系统都必定会有Geographic Coordinate System的参数。投影坐标系统，实质上便是平面坐标系统，其地图单位通常为米。投影：将球面坐标转化为平面坐标的过程便称为投影。即，投影的条件：
a、球面坐标
b、转化过程（也就是算法）
要得到投影坐标就必须得有一个“拿来”投影的球面坐标，然后才能使用算法去投影

from osgeo import gdal, osrdef Lonlat2Xy(SourceGcs, TargetPcs, lon, lat):''':param SourceRef: 源地理坐标系统:param TargetRef: 目标投影:param lon: 待转换点的longitude值:param lat: 待转换点的latitude值:return:'''#创建目标空间参考spatialref_target=osr.SpatialReference()spatialref_target.ImportFromEPSG(TargetPcs) #2331为目标空间参考的ESPG编号，西安80 高斯可吕格投影#创建原始空间参考spatialref_source=osr.SpatialReference()spatialref_source.ImportFromEPSG(SourceGcs)  #4326 为原始空间参考的ESPG编号，WGS84#构建坐标转换对象，用以转换不同空间参考下的坐标trans=osr.CoordinateTransformation(spatialref_source,spatialref_target)# coordinate_after_trans 是一个Tuple类型的变量包含3个元素， [0]为y方向值，[1]为x方向值，[2]为高度coordinate_after_trans=trans.TransformPoint(lat,lon)# print(coordinate_after_trans)#以下为转换多个点（要使用list）coordinate_trans_points=trans.TransformPoints([(40,117),(36,120)])print(coordinate_trans_points)return coordinate_after_trans
if __name__=='__main__':point_trans = Lonlat2Xy(4326, 2331, 120, 36)# 4326 为原始空间参考的ESPG编号  2331为目标空间参考的ESPG编号# 120为经度，36为纬度print(point_trans)

运行结果：

[(4588003.080066519, 19041214.81787683, 0.0), (4196378.172222488, 19406098.190200843, 0.0)]
(4196378.172222488, 19406098.190200843, 0.0)

2.投影坐标xy转经纬度

代码如下（示例）：

from osgeo import gdal, osrfilename=r'F:\GOSIF_2001065.tif'def Xy2Lonlat(SourcePcs, x, y):'''将投影坐标转为经纬度坐标（具体的投影坐标系由给定数据确定）:param SourcePcs:源投影坐标系:param x: 投影坐标x:param y: 投影坐标y:return: 投影坐标(x, y)对应的经纬度坐标(lon, lat)'''prosrs = osr.SpatialReference()# prosrs投影参考系prosrs.ImportFromEPSG(SourcePcs)#投影坐标的EPSG编号#geosrs地理参考系geosrs = osr.SpatialReference()geosrs.ImportFromEPSG(4326)            #WGS84ct = osr.CoordinateTransformation(prosrs, geosrs)coords = ct.TransformPoint(y, x)return coords[:2]   #(lat,lon)

3.行列号转投影或地理坐标

根据原始tiff的坐标系，确定是投影坐标还是地理坐标。

from osgeo import gdal, osrfilename=r'F:\GOSIF_2001065.tif'#tiff是投影坐标的情况未实测
def Rowcol2Lonlat(Xpixel,Ypixel):dataset = gdal.Open(filename)  # 打开文件GeoTransform = dataset.GetGeoTransform()XGeo = GeoTransform[0]+GeoTransform[1]*Xpixel+Ypixel*GeoTransform[2];YGeo = GeoTransform[3]+GeoTransform[4]*Xpixel+Ypixel*GeoTransform[5];pcs = osr.SpatialReference()pcs.ImportFromWkt(dataset.GetProjection())   if pcs.IsGeographic():    # 是地理坐标return XGeo,YGeoelif pcs.IsProjected():   #是投影坐标coords = Xy2Lonlat(pcs, XGeo, YGeo)  return coords[:2]

4.投影或地理坐标转行列号

根据原始tiff的坐标系，确定是投影坐标还是地理坐标。

from osgeo import gdal, osrfilename=r'F:\GOSIF_2001065.tif'#tiff是投影坐标的情况未实测
def Lonlat2Rowcol(lon,lat):dataset = gdal.Open(filename)  # 打开文件tiff_geotrans = dataset.GetGeoTransform()pcs = osr.SpatialReference()pcs.ImportFromWkt(dataset.GetProjection())if pcs.IsGeographic():  # 是地理坐标XGeo, YGeo = lon , latelif pcs.IsProjected():  # 是投影坐标YGeo, XGeo = Lonlat2Xy(4326, pcs, lon, lat)A = [[tiff_geotrans[1], tiff_geotrans[2]],  # 根据公式Xgeo=tiff_geotrans[0]+Xpixel*tiff_geotrans[1]+Yline*tiff_geotrans[2][tiff_geotrans[4], tiff_geotrans[5]]]  # Ygeo=tiff_geotrans[3]+Xpixel*tiff_geotrans[4]+Yline*tiff_geotrans[5]s = [[XGeo - tiff_geotrans[0]],  # 运用矩阵解二元一次方程组求得行列号[YGeo - tiff_geotrans[3]]]r = np.linalg.solve(A, s)# Xpixel, Ypixel = r[0], r[1]Xpixel,Ypixel = int(r[0]),int(r[1])return Xpixel,Ypixel

三、获取指定位置的值

def GetValueByCoordinates(Cox,Coy, coordinate_type):''':param Cox: x轴方向坐标:param Coy: y轴方向坐标:param coordinate_type: 'rowcol'，'lonlat'，'xy'3个取值，表示3种坐标:return: '''img, dataset, gcs, pcs, extend, shape = GetTifInfo()if coordinate_type == 'rowcol':value = img[Coy, Cox]elif coordinate_type == 'lonlat':col, row = Lonlat2Rowcol(Cox,Coy)value = img[row, col]elif coordinate_type == 'xy':row, col =Lonlat2Rowcol(Cox,Coy)value = img[row, col]Exception("coordinated_type error")return value

四、使用SHP文件裁剪TIFF文件

说明待补充
注意shp的坐标系应该和tiff文件一致，可在arcGIS中修改shp的坐标系

from osgeo import gdal, osr, ogr
import numpy as np
import os
import shapefilefilename=r'F:\GOSIF_2001065.tif'
shpfile = r'F:\门源县.shp'
tiff_path=r'F:\\'
out_path = r'F:\haibei_GOSIF'def TiffCutByShp():file_list = os.listdir(tiff_path)for file in file_list:if file.endswith('.tif'):print('Processing >>> ' + file)in_raster = gdal.Open(os.path.join(tiff_path, file))out_raster = os.path.join(out_path, file)r = shapefile.Reader(shpfile)ds = gdal.Warp(out_raster,# 裁剪图像保存完整路径（包括文件名）in_raster,# 待裁剪的影像format='GTiff', # 保存图像的格式cutlineDSName=shpfile,# 矢量文件的完整路径outputBounds=r.bbox,#shp文件的外包矩阵# cutlineWhere="FIELD = 'whatever'", # clip specific featuredstNodata=0)  # set nodata valueds = None  # close datasetelse:print("No '.tif' file found ...")