python基础教程:python+gdal+遥感图像拼接(mosaic)的实
这篇文章主要介绍了python+gdal+遥感图像拼接(mosaic)的实例,具有很好的参考价值,希望对大家有所帮助。一起跟随小编过来看看吧
作为摄影测量与遥感的从业者,笔者最近开始深入研究gdal,为工作打基础!个人觉得gdal也是没有什么技术含量,调用别人的api。但是想想这也是算法应用的一个技能,多学无害!
关于遥感图像的镶嵌,主要分为6大步骤:
step1:
1)对于每一幅图像,计算其行与列;
2)获取左上角X,Y
3)获取像素宽和像素高
4)计算max X 和 min Y,切记像素高是负值
maxX1 = minX1 + (cols1 * pixelWidth)
minY1 = maxY1 + (rows1 * pixelHeight)
step2 :计算输出图像的min X ,max X,min Y,max Y
minX = min(minX1, minX2, …)
maxX = max(maxX1, maxX2, …)
y坐标同理
step3:计算输出图像的行与列
cols = int((maxX – minX) / pixelWidth)
rows = int((maxY – minY) / abs(pixelHeight)
step 4:创建一个输出图像
driver.create()
step 5:
1)计算每幅图像左上角坐标在新图像的偏移值
2)依次读入每幅图像的数据并利用1)计算的偏移值将其写入新图像中
step6 :对于输出图像
1)刷新磁盘并计算统计值
2)设置输出图像的几何和投影信息
3)建立金字塔
下面附上笔者的代码:
#mosica 两张图像
import os, sys, gdal
from gdalconst import *
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(Geospatial Data Abstraction Library)是一个在X/MIT许可协议下的开源栅格空间数据转换库。它利用抽象数据模型来表达所支持的各种文件格式。它还有一系列命令行工具来进行数据转换和处理。
这个开源栅格空间数据转换库拥有许多和其他语言的接口,对于python,他有对应的第三方包GDAL,下载安装已在上篇文章中提到。
目的: 可以使用Python的第三方包:GDAL进行遥感数据的读写,方便批处理。
0.2 配置环境
电脑系统: win7x64
Python版本: 3.6.4
GDAL版本: 2.3.2
1 读
1.1 TIFF格式
标签图像文件格式(Tag Image File Format,简写为TIFF)是一种灵活的位图格式,主要用来存储包括照片和艺术图在内的图像。它最初由Aldus公司与微软公司一起为PostScript打印开发。TIFF与JPEG和PNG一起成为流行的高位彩色图像格式。
TIFF文件以.tif为扩展名。
def tif_read(tifpath, bandnum):"""Use GDAL to read data and transform them into arrays.:param tifpath:tif文件的路径:param bandnum:需要读取的波段:return:该波段的数据,narray格式。len(narray)是行数,len(narray[0])列数"""image = gdal.Open(tifpath) # 打开该图像if image == None:print(tifpath + "该tif不能打开!")returnlie = 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) # Get the information of band num.band_array = band.ReadAsArray(0,0,lie,hang) # Getting data from zeroth rows and 0 columns# band_df = pd.DataFrame(band_array)del image # 减少冗余return band_array, im_proj, im_geotrans
2 写
2.1 TIFF格式
TIFF格式的数据格式有:Byete、int16、uint16、int32、uint32、float32、float64等7余种。
首先,要判断数据的格式,才能按需求写出
def tif_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:param filename: 存出文件名:param im_data: 输入数据:param im_proj: 投影信息:param im_geotrans: 放射变换信息:return: 0 """if 'int8' in im_data.dtype.name: # 判断栅格数据的数据类型datatype = gdal.GDT_Byteelif 'int16' in im_data.dtype.name:datatype = gdal.GDT_UInt16else:datatype = gdal.GDT_Float32# 判读数组维数if len(im_data.shape) == 3:im_bands, im_height, im_width = im_data.shapeelse: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) #写入投影if im_bands == 1:dataset.GetRasterBand(1).WriteArray(im_data) #写入数组数据else:for i in range(im_bands):dataset.GetRasterBand(i+1).WriteArray(im_data[i])del dataset
3 展示
3.1 TIFF格式
# 这个展示的效果并不是太好,当做示意图用def tif_display(self,im_data):""":param im_data: 影像数据,narray:return: 展出影像"""# plt.imshow(im_data,'gray') # 必须规定为显示的为什么图像plt.imshow(im_data) # 必须规定为显示的为什么图像plt.xticks([]), plt.yticks([]) # 隐藏坐标线plt.show() # 显示出来,不要也可以,但是一般都要了
以上这篇python+gdal+遥感图像拼接(mosaic)的实例就是小编分享给大家的全部内容了
写到这里,给大家推荐一个资源很全的python学习聚集地,点击进入,这里有资深程序员分享以前学习
心得,学习笔记,还有一线企业的工作经验,且给大家精心整理一份python零基础到项目实战的资料,
每天给大家讲解python最新的技术,前景,学习需要留言的小细节
python基础教程:python+gdal+遥感图像拼接(mosaic)的实相关推荐
- python基础教程-Python入门教程完整版(懂中文就能学会)
提取码:sjfo 目录大纲: 本套教程15天 学前环境搭建 1-3 天内容为Linux基础命令 4-13 天内容为Python基础教程 14-15 天内容为 飞机大战项目演练 视频概括: 第一阶段(1 ...
- 什么是python基础教程-python基础教程之python是什么?概念解析
Python,是一种面向对象的解释型计算机程序设计语言,由荷兰人Guido van Rossum于1989年发明,第一个公开发行版发行于1991年. Python是纯粹的自由软件, 源代码和解释器CP ...
- python办公自动化知识点_Python自动化办公知识点整理汇总|python基础教程|python入门|python教程...
https://www.xin3721.com/eschool/pythonxin3721/ 知乎上有人提问:用python进行办公自动化都需要学习什么知识呢? 这可能是很多非IT职场人士面临的困 ...
- python下载电影天堂视频教程_一篇文章教会你利用Python网络爬虫获取电影天堂视频下载链接|python基础教程|python入门|python教程...
https://www.xin3721.com/eschool/pythonxin3721/ [一.项目背景] 相信大家都有一种头疼的体验,要下载电影特别费劲,对吧?要一部一部的下载,而且不能直观的知 ...
- python基础教程-Python基础教程,Python入门教程(非常详细)
Python 英文本意为"蟒蛇",直到 1989 年荷兰人 Guido van Rossum (简称 Guido)发明了一种面向对象的解释型编程语言(后续会介绍),并将其命名为 P ...
- python添加库详细教程_Python 中如何自动导入缺失的库?|python基础教程|python入门|python教程...
https://www.xin3721.com/eschool/pythonxin3721/ 在写 Python 项目的时候,我们可能经常会遇到导入模块失败的错误:ImportError: No mo ...
- python3.6.2下载教程_Windows下升级Python3.7.7后(原Python3.6.2版本)如何切换Python版本|python基础教程|python入门|python教程...
https://www.xin3721.com/eschool/pythonxin3721/ 笔者:风起怨江南 出处:https://www.cnblogs.com/mengjinxiang 笔者原创 ...
- Python基础教程:Python pass语句详解
2019独角兽企业重金招聘Python工程师标准>>> Python pass 语句 Python pass是空语句,是为了保持程序结构的完整性. pass 不做任何事情,一般用做占 ...
- python中ioerror怎么解决_Python IOError错误异常原因|python基础教程|python入门|python教程...
https://www.xin3721.com/eschool/pythonxin3721/ python语言IOError错误一般多发生在对文件操作报错时,表示要打开的文件不存在,当然能引发IOEr ...
最新文章
- playframework学习笔记1 -- 开发环境和第一个工程
- 传道、授业、解惑:俞士纶院长参加数据科学研究院第五届“院长接待日”
- tableau实战系列(二十八)-以可视化的方式打开关联分析算法购物篮分析(Market Basket Analysis)
- 浅谈5G机房配套那些事
- c语言常用指令翻译,c语言常见专业词汇带翻译
- 语音识别热词_出门问问 TWS 耳机语音交互解决方案
- 冬季,拿什么来温暖你的心情
- Go 学习笔记(11):切片
- mysql大数据量的全量备份_mysql备份神器——Xtrabackup全量备份还原
- shapenet各类数据(转载)
- CMMI与Agile敏捷开发比较之一:两者的本质区别
- Koa v2.x 中文文档 上下文(Context)
- Transform.GetComponentsInChildRen()
- 根据业务单生成时将描扫记录触发到临时表(SQL触发)
- 一个计算机专业女孩的求学之路——七年之痒,痒之感悟
- URL Schemes入门
- 坚果pro android版本,坚果Pro有几个版本 哪个版本好?坚果Pro各版本的区别
- 多多情报通:拼多多视频上传多久审核?如何发布新品?
- GeoHash算法获取附近店铺和距离
- Android 第三方应用接入微信平台(1)