python+gdal+遥感图像拼接(mosaic)的实例
作为摄影测量与遥感的从业者，笔者最近开始深入研究gdal，为工作打基础！个人觉得gdal也是没有什么技术含量，调用别人的api。但是想想这也是算法应用的一个技能，多学无害！
关于遥感图像的镶嵌，主要分为6大步骤：
step1：
1）对于每一幅图像，计算其行与列；
2）获取左上角X,Y
3）获取像素宽和像素高
4）计算maxX和minY，切记像素高是负值
maxX1=minX1+(cols1*pixelWidth)
minY1=maxY1+(rows1*pixelHeight）
step2：计算输出图像的minX,maxX,minY,maxY
minX=min(minX1,minX2,…)
maxX=max(maxX1,maxX2,…)
y坐标同理
step3：计算输出图像的行与列
cols=int((maxX–minX)/pixelWidth)
rows=int((maxY–minY)/abs(pixelHeight)
step4：创建一个输出图像
driver.create()
step5：
1）计算每幅图像左上角坐标在新图像的偏移值
2）依次读入每幅图像的数据并利用1)计算的偏移值将其写入新图像中
step6：对于输出图像
1）刷新磁盘并计算统计值
2）设置输出图像的几何和投影信息
3）建立金字塔
下面附上笔者的代码：
#mosica两张图像
importos,sys,gdal
fromgdalconstimport*
os.chdir('c:/temp/****')#改变文件夹路径
#注册gdal(required)
gdal.AllRegister()
#读入第一幅图像
ds1=gdal.Open('**.img')
band1=ds1.GetRasterBand(1)
rows1=ds1.RasterYSize
cols1=ds1.RasterXSize
#获取图像角点坐标
transform1=ds1.GetGeoTransform()
minX1=transform1[0]
maxY1=transform1[3]
pixelWidth1=transform1[1]
pixelHeight1=transform1[5]#是负值（important）
maxX1=minX1+(cols1*pixelWidth1)
minY1=maxY1+(rows1*pixelHeight1)
#读入第二幅图像
ds2=gdal.Open('**.img')
band2=ds2.GetRasterBand(1)
rows2=ds2.RasterYSize
cols2=ds2.RasterXSize
#获取图像角点坐标
transform2=ds2.GetGeoTransform()
minX2=transform2[0]
maxY2=transform2[3]
pixelWidth2=transform2[1]
pixelHeight2=transform2[5]
maxX2=minX2+(cols2*pixelWidth2)
minY2=maxY2+(rows2*pixelHeight2)
#获取输出图像坐标
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=ds1.GetDriver()
dsOut=driver.Create('mosiac.img',cols,rows,1,band1.DataType)#1是bands，默认
bandOut=dsOut.GetRasterBand(1)
#读图1的数据并将其写到输出图像中
data1=band1.ReadAsArray(0,0,cols1,rows1)
bandOut.WriteArray(data1,xOffset1,yOffset1)
#读图2的数据并将其写到输出图像中
data2=band2.ReadAsArray(0,0,cols2,rows2)
bandOut.WriteArray(data2,xOffset2,yOffset2)
'''写图像步骤'''
#统计数据
bandOut.FlushCache()#刷新磁盘
stats=bandOut.GetStatistics(0,1)#第一个参数是1的话，是基于金字塔统计，第二个
#第二个参数是1的话：整幅图像重度，不需要统计
#设置输出图像的几何信息和投影信息
geotransform=[minX,pixelWidth1,0,maxY,0,pixelHeight1]
dsOut.SetGeoTransform(geotransform)
dsOut.SetProjection(ds1.GetProjection())
#建立输出图像的金字塔
gdal.SetConfigOption('HFA_USE_RRD','YES')
dsOut.BuildOverviews(overviewlist=[2,4,8,16])#4层
补充知识：运用Python的第三方库：GDAL进行遥感数据的读写
0背景及配置环境
0.1背景
GDAL(GeospatialDataAbstractionLibrary)是一个在X/MIT许可协议下的开源栅格空间数据转换库。它利用抽象数据模型来表达所支持的各种文件格式。它还有一系列命令行工具来进行数据转换和处理。
这个开源栅格空间数据转换库拥有许多和其他语言的接口，对于python,他有对应的第三方包GDAL，下载安装已在上篇文章中提到。
目的:可以使用Python的第三方包：GDAL进行遥感数据的读写，方便批处理。
0.2配置环境
电脑系统：win7x64
Python版本:3.6.4
GDAL版本:2.3.2
1读
1.1TIFF格式
标签图像文件格式（TagImageFileFormat，简写为TIFF）是一种灵活的位图格式，主要用来存储包括照片和艺术图在内的图像。它最初由Aldus公司与微软公司一起为PostScript打印开发。TIFF与JPEG和PNG一起成为流行的高位彩色图像格式。
TIFF文件以.tif为扩展名。
deftif_read(tifpath,bandnum):
"""
UseGDALtoreaddataandtransformthemintoarrays.
:paramtifpath:tif文件的路径
:parambandnum:需要读取的波段
:return:该波段的数据，narray格式。len(narray)是行数，len(narray[0])列数
"""
image=gdal.Open(tifpath)#打开该图像
ifimage==None:
print(tifpath+"该tif不能打开!")
return
lie=image.RasterXSize#栅格矩阵的列数
hang=image.RasterYSize#栅格矩阵的行数
im_bands=image.RasterCount#波段数
im_proj=image.GetProjection()#获取投影信息
im_geotrans=image.GetGeoTransform()#仿射矩阵
print('该tif:{}个行，{}个列，{}层波段,取出第{}层.'.format(hang,lie,im_bands,bandnum))
band=image.GetRasterBand(bandnum)#Gettheinformationofbandnum.
band_array=band.ReadAsArray(0,0,lie,hang)#Gettingdatafromzerothrowsand0columns
#band_df=pd.DataFrame(band_array)
delimage#减少冗余
returnband_array,im_proj,im_geotrans
2写
2.1TIFF格式
TIFF格式的数据格式有：Byete、int16、uint16、int32、uint32、float32、float64等7余种。
首先，要判断数据的格式，才能按需求写出。
deftif_write(self,filename,im_data,im_proj,im_geotrans):
"""
gdal数据类型包括
gdal.GDT_Byte,
gdal.GDT_UInt16,gdal.GDT_Int16,gdal.GDT_UInt32,gdal.GDT_Int32,
gdal.GDT_Float32,gdal.GDT_Float64
:paramfilename:存出文件名
:paramim_data:输入数据
:paramim_proj:投影信息
:paramim_geotrans:放射变换信息
:return:0
"""
if'int8'inim_data.dtype.name:#判断栅格数据的数据类型
datatype=gdal.GDT_Byte
elif'int16'inim_data.dtype.name:
datatype=gdal.GDT_UInt16
else:
datatype=gdal.GDT_Float32
#判读数组维数
iflen(im_data.shape)==3:
im_bands,im_height,im_width=im_data.shape
else:
im_bands,(im_height,im_width)=1,im_data.shape#多维或1.2维
#创建文件
driver=gdal.GetDriverByName("GTiff")#数据类型必须有，因为要计算需要多大内存空间
dataset=driver.Create(filename,im_width,im_height,im_bands,datatype)
dataset.SetGeoTransform(im_geotrans)#写入仿射变换参数
dataset.SetProjection(im_proj)#写入投影
ifim_bands==1:
dataset.GetRasterBand(1).WriteArray(im_data)#写入数组数据
else:
foriinrange(im_bands):
dataset.GetRasterBand(i+1).WriteArray(im_data[i])
deldataset
3展示
3.1TIFF格式
#这个展示的效果并不是太好，当做示意图用 deftif_display(self,im_data): """ :paramim_data:影像数据,narray :return:展出影像 """ #plt.imshow(im_data,'gray')#必须规定为显示的为什么图像 plt.imshow(im_data)#必须规定为显示的为什么图像 plt.xticks([]),plt.yticks([])#隐藏坐标线 plt.show()#显示出来，不要也可以，但是一般都要了
以上这篇python+gdal+遥感图像拼接(mosaic)的实例就是小编分享给大家的全部内容了，希望能给大家一个参考，也希望大家多多支持毛票票。
声明：本文内容来源于网络，版权归原作者所有，内容由互联网用户自发贡献自行上传，本网站不拥有所有权，未作人工编辑处理，也不承担相关法律责任。如果您发现有涉嫌版权的内容，欢迎发送邮件至：czq8825#qq.com（发邮件时，请将#更换为@）进行举报，并提供相关证据，一经查实，本站将立刻删除涉嫌侵权内容。