一些基于python+gdal整理的小工具
#!/usr/bin/env python
# coding: utf-8
from osgeo import gdal
import os
import glob
import numpy as np
import math
def read_img( filename):
'''
读取影像为数组并返回信息
——————————————————————————
@param
filename ：输入的影像路径
@return:
影像的numpy数组格式，并显示影像的基本信息
'''
dataset = gdal.Open(filename) #打开文件
print('栅格矩阵的列数:' ,dataset.RasterXSize)
print('栅格矩阵的行数:',dataset.RasterYSize)
print('波段数:',dataset.RasterCount)
print('数据类型:',dataset.GetRasterBand(1).DataType)
print('仿射矩阵(左上角像素的大地坐标和像素分辨率)',dataset.GetGeoTransform())
print('地图投影信息:',dataset.GetProjection())
im_data = dataset.ReadAsArray()
del dataset
return im_data
def read2arr( filename): # 读取影像为数组
'''
读取影像为数组
——————————————————————————
@param
filenam ：输入的影像路径
@return:
影像的numpy数组格式
'''
dataset = gdal.Open(filename) #打开文件
im_data = dataset.ReadAsArray()
del dataset
return im_data
def ds2tif(dataset,out_fn):
"""
将GDAL dataset数据格式写入tif保存
——————————————————————————
@param：
dataset：输入的GDAL影像数据格式
out_fn：输出的文件路径
@return：
输出影像文件
"""
# 读取dataset信息
im_width = dataset.RasterXSize
im_height = dataset.RasterYSize
im_bands = dataset.RasterCount
im_geotrans = dataset.GetGeoTransform()
im_proj = dataset.GetProjection()
im_datatype = dataset.GetRasterBand(1).DataType
# 将dataset 写入 tif
driver = gdal.GetDriverByName("GTiff")
ds = driver.Create(out_fn ,im_width, im_height,im_bands,im_datatype)
ds.SetGeoTransform(im_geotrans)
ds.SetProjection(im_proj)
ds.GetRasterBand(1).WriteArray(dataset.ReadAsArray())
del ds
def arr2tif(arr ,out_fn, Transform=None,Projection=None,Band = 1,Datatype=6):
"""
将数组格式写入tif保存
——————————————————————————
@param：
arr：待保存的影像数组
out_fn：输出的文件路径
Transform：仿射矩阵六参数数组，默认为空,详细格式见GDAL。
Projection ：投影，默认为空,详细格式见GDA
Band ：波段数，默认为1
Datatype：保存数据格式（位深），默认为6，GDT_Float32
@return：
输出影像文件
"""
(x,y) = arr.shape
driver = gdal.GetDriverByName("GTiff")
ds = driver.Create(out_fn ,y, x,Band, Datatype)
if not Transform==None:
ds.SetGeoTransform(Transform)
if not Projection==None:
ds.SetProjection(Projection)
ds.GetRasterBand(Band).WriteArray(arr)
del ds
def GetExtent(in_fn):
'''
计算影像角点的地理坐标或投影坐标
——————————————————————————
@param：
影像路径
@return:
min_x： x方向最小值
max_y： y方向最大值
max_x： x方向最大值
min_y: y方向最小值
'''
ds=gdal.Open(in_fn)
geotrans=list(ds.GetGeoTransform())
xsize=ds.RasterXSize
ysize=ds.RasterYSize
min_x=geotrans[0]
max_y=geotrans[3]
max_x=geotrans[0]+xsize*geotrans[1]
min_y=geotrans[3]+ysize*geotrans[5]
ds=None
return min_x,max_y,max_x,min_y
def DsGetExtent(ds):
'''
读取dataset格式，计算影像角点的地理坐标或投影坐标
——————————————————————————
@param：
ds： GDAL dataset格式数据
@return:
min_x： x方向最小值
max_y： y方向最大值
max_x： x方向最大值
min_y: y方向最小值
'''
geotrans=list(ds.GetGeoTransform())
xsize=ds.RasterXSize
ysize=ds.RasterYSize
min_x=geotrans[0]
max_y=geotrans[3]
max_x=geotrans[0]+xsize*geotrans[1]
min_y=geotrans[3]+ysize*geotrans[5]
ds=None
return min_x,max_y,max_x,min_y
def pix2geo(Xpixel, Ypixel, GeoTransform):
'''
计算影像某一像素点的地理坐标或投影坐标
——————————————————————————
@param：
Xpixel ：像素坐标x
Ypixel： 像素坐标y
GeoTransform：仿射变换参数
@return:
XGeo： 地理坐标或投影坐标X
YGeo： 地理坐标或投影坐标Y
'''
XGeo = GeoTransform[0] + GeoTransform[1] * Xpixel + Ypixel * GeoTransform[2]
YGeo = GeoTransform[3] + GeoTransform[4] * Xpixel + Ypixel * GeoTransform[5]
return XGeo, YGeo
def geo2pix(dataset, x, y):
'''
根据GDAL的仿射变换参数模型将给定的投影或地理坐标转为影像图上坐标（行列号）
————————————————————————————————
@param
dataset: GDAL地理数据
x: 投影或地理坐标x
y: 投影或地理坐标y
@return:
影坐标或地理坐标(x, y)对应的影像图上行列号(row, col)
'''
trans = dataset.GetGeoTransform()
a = np.array([[trans[1], trans[2]], [trans[4], trans[5]]])
b = np.array([x - trans[0], y - trans[3]])
return np.linalg.solve(a, b)
def Mosaic_all (path):
'''
将指定路径文件夹下的tif影像全部镶嵌到一张影像上
————————————————————————————————
@param
path：tif影像存放路径
@return:
镶嵌合成的整体影像
'''
os.chdir(path)
#如果存在同名影像则先删除
if os.path.exists('mosaiced_image.tif'):
os.remove('mosaiced_image.tif')
in_files=glob.glob("*.TIF")
in_fn=in_files[0]
#获取待镶嵌栅格的最大最小的坐标值
min_x,max_y,max_x,min_y=GetExtent(in_fn)
for in_fn in in_files[1:]:
minx,maxy,maxx,miny=GetExtent(in_fn)
min_x=min(min_x,minx)
min_y=min(min_y,miny)
max_x=max(max_x,maxx)
max_y=max(max_y,maxy)
#计算镶嵌后影像的行列号
in_ds=gdal.Open(in_files[0])
geotrans=list(in_ds.GetGeoTransform())
width=geotrans[1]
height=geotrans[5]
columns=math.ceil((max_x-min_x)/width)
rows=math.ceil((max_y-min_y)/(-height))
in_band=in_ds.GetRasterBand(1)
driver=gdal.GetDriverByName('GTiff')
out_ds=driver.Create('mosaiced_image.tif',columns,rows,1,in_band.DataType)
out_ds.SetProjection(in_ds.GetProjection())
geotrans[0]=min_x
geotrans[3]=max_y
out_ds.SetGeoTransform(geotrans)
out_band=out_ds.GetRasterBand(1)
#定义仿射逆变换
inv_geotrans=gdal.InvGeoTransform(geotrans)
#开始逐渐写入
for in_fn in in_files:
in_ds=gdal.Open(in_fn)
in_gt=in_ds.GetGeoTransform()
#仿射逆变换
offset=gdal.ApplyGeoTransform(inv_geotrans,in_gt[0],in_gt[3])
x,y=map(int,offset)
print(x,y)
trans=gdal.Transformer(in_ds,out_ds,[]) #in_ds是源栅格，out_ds是目标栅格
success,xyz=trans.TransformPoint(False,0,0) #计算in_ds中左上角像元对应out_ds中的行列号
x,y,z=map(int,xyz)
print(x,y,z)
data=in_ds.GetRasterBand(1).ReadAsArray()
out_band.WriteArray(data,x,y) #x，y是开始写入时左上角像元行列号
del in_ds,out_band,out_ds
return 0
def Mosaic(ds1 , ds2, path):
'''
将两幅影像镶嵌至同一幅影像
————————————————————————————————
@param
ds1：镶嵌数据集1
ds2：镶嵌数据集1
@return:
镶嵌合成的整体影像
'''
band1 = ds1.GetRasterBand(1)
rows1 = ds1.RasterYSize
cols1 = ds1.RasterXSize
band2 = ds2.GetRasterBand(1)
rows2 = ds2.RasterYSize
cols2 = ds2.RasterXSize
(minX1,maxY1,maxX1,minY1) = DsGetExtent(ds1)
(minX2,maxY2,maxX2,minY2) = DsGetExtent(ds2)
transform1 = ds1.GetGeoTransform()
pixelWidth1 = transform1[1]
pixelHeight1 = transform1[5] #是负值（important）
transform2 = ds2.GetGeoTransform()
pixelWidth2 = transform1[1]
pixelHeight2 = transform1[5]
# 获取输出图像坐标
minX = min(minX1, minX2)
maxX = max(maxX1, maxX2)
minY = min(minY1, minY2)
maxY = max(maxY1, maxY2)
#获取输出图像的行与列
cols = int((maxX - minX) / pixelWidth1)
rows = int((maxY - minY) / abs(pixelHeight1))
# 计算图1左上角的偏移值（在输出图像中）
xOffset1 = int((minX1 - minX) / pixelWidth1)
yOffset1 = int((maxY1 - maxY) / pixelHeight1)
# 计算图2左上角的偏移值（在输出图像中）
xOffset2 = int((minX2 - minX) / pixelWidth1)
yOffset2 = int((maxY2 - maxY) / pixelHeight1)
# 创建一个输出图像
driver = gdal.GetDriverByName("GTiff")
out_ds = driver.Create( path, cols, rows, 1, band1.DataType)#1是bands，默认
out_band = out_ds.GetRasterBand(1)
# 读图1的数据并将其写到输出图像中
data1 = band1.ReadAsArray(0, 0, cols1, rows1)
out_band.WriteArray(data1, xOffset1, yOffset1)
#读图2的数据并将其写到输出图像中
data2 = band2.ReadAsArray(0, 0, cols2, rows2)
out_band.WriteArray(data2, xOffset2, yOffset2)
''' 写图像步骤'''
#第二个参数是1的话：整幅图像重度，不需要统计
# 设置输出图像的几何信息和投影信息
geotransform = [minX, pixelWidth1, 0, maxY, 0, pixelHeight1]
out_ds.SetGeoTransform(geotransform)
out_ds.SetProjection(ds1.GetProjection())
del ds1,ds2,out_band,out_ds,driver
return 0
def raster_overlap(ds1, ds2, nodata1=None, nodata2=None):
'''
两个栅格数据集取重叠区或求交集 （仅测试方形影像）
————————————————————————————————
@param:
ds1 (GDAL dataset) - GDAL dataset of an image
ds2 (GDAL dataset) - GDAL dataset of an image
nodata1 (number) - nodata value of image 1
nodata2 (number) - nodata value of image 2
@returns:
ds1c (GDAL dataset), ds2c (GDAL dataset): intersection dataset1 and intersection dataset2.
'''
##Setting nodata
nodata = 0
###Check if images NoData is set
if nodata2 is not None:
nodata = nodata2
ds2.GetRasterBand(1).SetNoDataValue(nodata)
else:
if ds2.GetRasterBand(1).GetNoDataValue() is None:
ds2.GetRasterBand(1).SetNoDataValue(nodata)
if nodata1 is not None:
nodata = nodata1
ds1.GetRasterBand(1).SetNoDataValue(nodata1)
else:
if ds1.GetRasterBand(1).GetNoDataValue() is None:
ds1.GetRasterBand(1).SetNoDataValue(nodata)
### Get extent from ds1
projection = ds1.GetProjection()
geoTransform = ds1.GetGeoTransform()
###Get minx and max y
[minx, maxy, maxx, miny] = DsGetExtent(ds1)
[minx_2, maxy_2, maxx_2, miny_2] = DsGetExtent(ds2)
min_x = sorted([maxx,minx_2,minx,maxx_2])[1] # 对边界值排序，第二三个为重叠区边界
max_y = sorted([maxy,miny_2,miny,maxy_2])[2]
max_x = sorted([maxx,minx_2,minx,maxx_2])[2]
min_y = sorted([maxy,miny_2,miny,maxy_2])[1]
###Warp to same spatial resolution
gdaloptions = {'format': 'MEM', 'xRes': geoTransform[1], 'yRes':
geoTransform[5], 'dstSRS': projection}
ds2w = gdal.Warp('', ds2, **gdaloptions)
ds2 = None
###Translate to same projection
ds2c = gdal.Translate('', ds2w, format='MEM', projWin=[min_x, max_y, max_x, min_y],
outputSRS=projection)
ds2w = None
ds1c = gdal.Translate('', ds1, format='MEM', projWin=[min_x, max_y, max_x, min_y],
outputSRS=projection)
ds1 = None
return ds1c,ds2c
本文章为转载内容，我们尊重原作者对文章享有的著作权。如有内容错误或侵权问题，欢迎原作者联系我们进行内容更正或删除文章。
