Python处理 TIFF 格式数据
简介
参考文献1
- GDAL 是一个开源的操作栅格数据和矢量数据的库，本文记录下用 Python 中 GDAL 库操作 TIFF （GeoTIFF）的常见代码，包括读写、获取坐标系、获取指定位置像元值等。
- TIFF 简单理解就是一种图像格式，类似于 jpg、png 等。
- GeoTIFF 就是在普通 TIFF 文件上增加了地理位置、投影信息、坐标信息等，常用于遥感数据。
一、读取TIFF文件
from osgeo import gdal, osr
filename=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
参考文献2
- 地理坐标系（Geographic coordinate system），Geographic coordinate system直译为地理坐标系统，是以经纬度为地图的存储单位的。很明显，Geographic coordinate system是球面坐标系统。
- Projection coordinate system（投影坐标系统），每一个投影坐标系统都必定会有Geographic Coordinate System的参数。投影坐标系统，实质上便是平面坐标系统，其地图单位通常为米。投影：将球面坐标转化为平面坐标的过程便称为投影。即，投影的条件：
a、球面坐标
b、转化过程（也就是算法）
要得到投影坐标就必须得有一个“拿来”投影的球面坐标，然后才能使用算法去投影
from osgeo import gdal, osr
def 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, osr
filename=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) #WGS84
ct = osr.CoordinateTransformation(prosrs, geosrs)
coords = ct.TransformPoint(y, x)
return coords[:2] #(lat,lon)
3.行列号转投影或地理坐标
根据原始tiff的坐标系，确定是投影坐标还是地理坐标。
from osgeo import gdal, osr
filename=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,YGeo
elif pcs.IsProjected(): #是投影坐标
coords = Xy2Lonlat(pcs, XGeo, YGeo)
return coords[:2]
4.投影或地理坐标转行列号
根据原始tiff的坐标系，确定是投影坐标还是地理坐标。
from osgeo import gdal, osr
filename=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 , lat
elif 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 shapefile
filename=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 feature
dstNodata=0) # set nodata value
ds = None # close dataset
else:
print("No '.tif' file found ...")
五、保存TIFF文件
未完待续
文章出处登录后可见！
已经登录？立即刷新
