Python干货 | 遥感影像拼接
时间:2022-07-26
本文章向大家介绍Python干货 | 遥感影像拼接,主要内容包括其使用实例、应用技巧、基本知识点总结和需要注意事项,具有一定的参考价值,需要的朋友可以参考一下。
0.前言
因为没有喝上“秋天的第一份奶茶”,准备来更新一篇推送。
在上一篇推文中,我展示了如何使用Python
结合Landsat
制作遥感影像图(Python干货 | 制作遥感影像图)。
对于Landsat
数据来说,对某个区域的重访周期为16天,每个位置使用全球参考系(WRS)进行索引,即每一个位置都会对应一个Path和Row,相邻的影像之间会有部分区域是重叠的。
Fig.1 World Reference System
在某些遥感影像的应用场景中,如果我们关注的区域正好处于两景影像的交界处,如下图中的象山港,那我们就需要将影像拼接起来才可以使用。
单张影像是这样。
本文合并后是这样。
1.准备工作
相较于上一篇推送,我们这次为了实现遥感影像的镶嵌拼接,我们使用到了两个库, rasterio
和gdal
。
import rasterio as rio
import gdal
先介绍一下我们实现两组遥感影像拼接的思路,首先选取两景相邻的影像,分别得到他们的空间范围,再得到两景组合到一起之后的空间范围,使用gdal
新建一个tif
文件(数据中转用),分别得到原来两景影像在新建的tif
文件中的起始位置,将对应的数据写入新的tif
文件中,即实现镶嵌拼接。
上面说的是两景影像的拼接,如果是更多影像拼接同样适用,但是现阶段的方法如果拼接多的影像的话,需要的内存空间很大,容易导致内存溢出,感兴趣的朋友可以思考一下如何高效实现多景影像的拼接。
其中还有两处关键处理,一是如何去除重叠区域的无效信息,二是重叠区域的数据如何选择。希望各位看官能从代码里面找到答案。
2.动起手来
得到输入影像的四个角点。
def tiffileList2filename(tiffileList):
filename = []
prefix = []
for ifile in tiffileList:
file0 = ifile.split("\")[-1]
prefix.append(os.path.join(ifile, file0))
filename.append(os.path.join(ifile, file0) + "_B1.TIF")
return filename, prefix
def get_extent(tiffileList):
filename, prefix = tiffileList2filename(tiffileList)
rioData = rio.open(filename[0])
left = rioData.bounds[0]
bottom = rioData.bounds[1]
right = rioData.bounds[2]
top = rioData.bounds[3]
for ifile in filename[1:]:
rioData = rio.open(ifile)
left = min(left, rioData.bounds[0])
bottom = min(bottom, rioData.bounds[1])
right = max(right, rioData.bounds[2])
top = max(top, rioData.bounds[3])
return left, bottom, right, top, filename, prefix
得到新建tif文件的size,这里已知Landsat
空间分辨率为30m,如果是其他遥感数据,需对应进行修改。
def getRowCol(left, bottom, right, top):
cols = int((right - left) / 30.0)
rows = int((top - bottom) / 30.0)
return cols, rows
主程序,其中plot_rgb为上一篇推送中用到的函数。
if __name__ == '__main__':
tiffileList = [r'PathofLandsat8LC08_L1TP_118039_20160126_20170330_01_T1',
r'PathofLandsat8LC08_L1TP_118040_20160126_20170330_01_T1']
left, bottom, right, top, filename, prefix = get_extent(tiffileList)
cols, rows= getRowCol(left, bottom, right, top)
bands = ['B7', 'B5', 'B3']
n_bands = len(bands)
arr = np.zeros((n_bands, rows, cols), dtype=np.float)
# 打开一个tif文件
in_ds = gdal.Open(filename[0])
for i in range(len(bands)):
ibands = bands[i]
# 新建一个tif文件
driver = gdal.GetDriverByName('gtiff')
out_ds = driver.Create(ibands + 'mosaic.tif', cols, rows)
# 设置tif文件的投影
out_ds.SetProjection(in_ds.GetProjection())
out_band = out_ds.GetRasterBand(1)
# 设置新tif文件的地理变换
gt = list(in_ds.GetGeoTransform())
gt[0], gt[3] = left, top
out_ds.SetGeoTransform(gt)
# 对要拼接的影像进行循环读取
for ifile in prefix:
in_ds = gdal.Open(ifile + '_' + ibands + '.tif')
# 计算新建的tif文件及本次打开的tif文件之间的坐标漂移
trans = gdal.Transformer(in_ds, out_ds, [])
# 得到偏移起始点
success, xyz = trans.TransformPoint(False, 0, 0)
x, y, z = map(int, xyz)
# 读取波段信息
fnBand = in_ds.GetRasterBand(1)
data = fnBand.ReadAsArray()
# 写入tif文件之前,最大值设置为255,这一步很关键
data = data / 65535 * 255
data[np.where(data == 255)] = 0
# 影像重合部分处理,重合部分取最大值
xSize = fnBand.XSize
ySize = fnBand.YSize
outData = out_band.ReadAsArray(x, y, xSize, ySize)
data = np.maximum(data, outData)
out_band.WriteArray(data, x, y)
del out_band, out_ds
file2read = ibands + 'mosaic.tif'
arr[i, :, :] = tiff.imread(file2read)
os.remove(file2read)
plot_rgb(arr, rgb=(0, 1, 2))
3.小结
大功告成!
- [接口测试 - http.client篇] 14 源码初探及其工作机制分析
- 接口测试 | 25 requests + pytest测试实例
- 接口测试 | 22 requests基础入门
- 图的存储结构的实现(C/C++实现)
- 树和二叉树的存储结构的实现(C/C++实现)
- Selenium3源码之异常模块篇
- 移位密码原理及算法实现
- 排序算法的实现(C/C++实现)
- [开源] 分享导出博客园文章成本地 Markdown 文件存储的工具
- 单表代替密码原理及算法实现
- 【Android开发学习笔记之一】5大布局方式详解
- Selenium3源码之common下action_chains.py模块分析
- 图的简单应用(C/C++实现)
- 一个很easy的脚本--php获取服务器端的相关信息
- JavaScript 教程
- JavaScript 编辑工具
- JavaScript 与HTML
- JavaScript 与Java
- JavaScript 数据结构
- JavaScript 基本数据类型
- JavaScript 特殊数据类型
- JavaScript 运算符
- JavaScript typeof 运算符
- JavaScript 表达式
- JavaScript 类型转换
- JavaScript 基本语法
- JavaScript 注释
- Javascript 基本处理流程
- Javascript 选择结构
- Javascript if 语句
- Javascript if 语句的嵌套
- Javascript switch 语句
- Javascript 循环结构
- Javascript 循环结构实例
- Javascript 跳转语句
- Javascript 控制语句总结
- Javascript 函数介绍
- Javascript 函数的定义
- Javascript 函数调用
- Javascript 几种特殊的函数
- JavaScript 内置函数简介
- Javascript eval() 函数
- Javascript isFinite() 函数
- Javascript isNaN() 函数
- parseInt() 与 parseFloat()
- escape() 与 unescape()
- Javascript 字符串介绍
- Javascript length属性
- javascript 字符串函数
- Javascript 日期对象简介
- Javascript 日期对象用途
- Date 对象属性和方法
- Javascript 数组是什么
- Javascript 创建数组
- Javascript 数组赋值与取值
- Javascript 数组属性和方法
- PHP添加文字水印或图片水印的水印类完整源代码与使用示例
- Python 分布式缓存之Reids数据类型操作详解
- Pycharm打开已有项目配置python环境的方法
- python cv2.resize函数high和width注意事项说明
- pytorch SENet实现案例
- python如何安装下载后的模块
- Python爬虫如何应对Cloudflare邮箱加密
- 如何使用Python处理HDF格式数据及可视化问题
- tp5框架使用composer实现日志记录功能示例
- python 图像插值 最近邻、双线性、双三次实例
- tp5(thinkPHP5)框架实现多数据库查询的方法
- Python-openCV开运算实例
- php curl获取https页面内容,不直接输出返回结果的设置方法
- 详解php中curl返回false的解决办法
- Pytorch mask-rcnn 实现细节分享