![f859567fea75404eb8cb574453026e30.png](https://i-blog.csdnimg.cn/blog_migrate/929facc63e4f943bce4b23fc4773a416.png)
两张待镶嵌影像
![fbf4caca133049f3abd7146e0d8e3f1d.png](https://i-blog.csdnimg.cn/blog_migrate/33b15389f9688b4f076b73f5a2fe9a22.png)
镶嵌和裁剪后的影像
(这里只用了两幅影像镶嵌,所以左边是平整的,完整来说需要六幅影像)
遥感影像镶嵌
影像镶嵌是指将多幅相邻的遥感影像拼接成一个连续的、无缝的大范围影像的过程。
在Python中,我们可以通过以下步骤实现影像镶嵌:
1. 创建“空”的镶嵌结果栅格
首先,我们需要在内存中创建一个“空”的镶嵌结果栅格,来接收我们的遥感数据。这个栅格有多大、几个波段、数据类型是什么都需要设置,以及给它添加坐标意义(投影、地理转换信息)。
def GetExtent(in_fn):
'''
这个函数主要用于获取一幅遥感影像的经纬度范围
'''
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
'''----------------第一步:在内存中创建“空”的镶嵌结果的空栅格:有多大、几个波段、数据类型是是什么、添加坐标意义(投影、地理转换信息)-------------------------'''
global out_ds
# 获取待镶嵌栅格的最大最小的坐标值
in_fn = tiffiles_path_list[0]
min_x, max_y, max_x, min_y = GetExtent(in_fn)
for in_fn in tiffiles_path_list[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)
print(f'镶嵌结果图的经纬度范围是:{min_x} {min_y} {max_x} {max_y}')
# 计算镶嵌后影像的行列数
in_ds = gdal.Open(tiffiles_path_list[0])
geotrans = list(in_ds.GetGeoTransform())
width = geotrans[1]
height = geotrans[5]
print(f'像素东西方向分辨率为{width},南北方向分辨率为{height}')
columns = math.ceil((max_x - min_x) / width)
rows = math.ceil((max_y - min_y) / (-height))
print(f'镶嵌结果图的行列数是:行:{rows} 列:{columns}')
in_band = in_ds.GetRasterBand(1)
driver = gdal.GetDriverByName('MEM') #创建驱动程序:创建的数据集保存到内存中,而不是磁盘
num_bands = in_ds.RasterCount
print(f'波段数为:{num_bands}')
print(f'第一张影像的投影信息是{in_ds.GetProjection()}')
out_ds = driver.Create('', columns, rows, num_bands, in_band.DataType) #这里只能设置输出结果的形状大小,波段等,现在是没有坐标意义的
# 给结果添加坐标意义
out_ds.SetProjection(in_ds.GetProjection())
geotrans[0] = min_x
geotrans[3] = max_y
out_ds.SetGeoTransform(geotrans)
2. 往空栅格里面写数据
接下来,我们需要将每一幅输入影像的数据写入到我们创建的“空”栅格中。这个过程需要对每一幅影像进行遍历,并计算其在镶嵌结果栅格中的位置。
'''------------------------------------------第二步:往空栅格里面写数据----------------------------------------------'''
for b in tqdm(range(1, num_bands + 1), desc='正在镶嵌真彩色近红外波段影像文件'): # 从所有影像的第一个波段开始镶嵌
for tif in tiffiles_path_list:
in_tif = gdal.Open(tif)
# 计算in_tif中左上角像元对应out_ds中的行列号
trans = gdal.Transformer(in_tif, out_ds, []) # in_tif是输入栅格,out_ds是目标栅格,第三个参数设置转换器选项,这里不需要
# success接收是否转换成功,False表示计算从源栅格到目标栅格的偏移
success, xyz = trans.TransformPoint(False, 0, 0)
x, y, z = map(int, xyz)
data = in_tif.GetRasterBand(b).ReadAsArray() #将in_tif第b波段读成数组
print(f'data:{data}')
# 找到不包含NaN值的行的索引
valid_rows_index = np.all(~np.isnan(data), axis=1)
# 提取有效的像元
valid_data = data[valid_rows_index]
print(f'valid_data:{valid_data}')
row_difference = len(data) - len(valid_data)
print(f'行差:{row_difference}')
if row_difference > 1:
row_difference = 1
print(f'子影像的有效值在空栅格中的像素坐标为第{y + row_difference}行,第{x}列')
# 将删除NaN后的有效像元写入到输出数据集中
out_ds.GetRasterBand(b).WriteArray(valid_data, x, y + row_difference)
# out_ds.GetRasterBand(b).WriteArray(data, x, y) # 将data从空栅格的第{y}行,第{x}列开始写入
del in_tif
遥感影像裁剪
在Python中,我们可以使用GDAL的gdal.Warp
函数来实现影像裁剪,对warp不了解的可以看我的这篇文章:python:GDAL的gdal.Warp功能及应用示例:裁剪与坐标系转换
如下所示:
input_raster = out_ds
# gdal.PushErrorHandler('CPLQuietErrorHandler')
result = gdal.Warp(
tiff_out_path,
input_raster,
format='GTiff',
cutlineDSName=geojson,
cropToCutline=True,
dstNodata=0
)
result.FlushCache()
完整代码:
tiffiles_path_list是一个有多个遥感影像路径的列表
geojson是矢量边界的路径,也可以换成shp
tiff_out_path是输出路径
import math
import numpy as np
from osgeo import gdal, ogr
from tqdm import tqdm
tiffiles_path_list = [r"F:\text\dongchengqu\Sentinel2_2023_K50E024011.tif",
r"F:\text\dongchengqu\Sentinel2_2023_J50E001011.tif"]
geojson = r"F:\地调系统\Test_Data(Beijing)\shp\朝阳区.geojson"
tiff_out_path = r'F:\text\result\ceshi4.tif'
def GetExtent(in_fn):
'''
这个函数主要用于获取一幅遥感影像的经纬度范围
'''
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
'''----------------第一步:在内存中创建“空”的镶嵌结果的空栅格:有多大、几个波段、数据类型是是什么、添加坐标意义(投影、地理转换信息)-------------------------'''
global out_ds
# 获取待镶嵌栅格的最大最小的坐标值
in_fn = tiffiles_path_list[0]
min_x, max_y, max_x, min_y = GetExtent(in_fn)
for in_fn in tiffiles_path_list[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)
print(f'镶嵌结果图的经纬度范围是:{min_x} {min_y} {max_x} {max_y}')
# 计算镶嵌后影像的行列数
in_ds = gdal.Open(tiffiles_path_list[0])
geotrans = list(in_ds.GetGeoTransform())
width = geotrans[1]
height = geotrans[5]
print(f'像素东西方向分辨率为{width},南北方向分辨率为{height}')
columns = math.ceil((max_x - min_x) / width)
rows = math.ceil((max_y - min_y) / (-height))
print(f'镶嵌结果图的行列数是:行:{rows} 列:{columns}')
in_band = in_ds.GetRasterBand(1)
driver = gdal.GetDriverByName('MEM') #创建驱动程序:创建的数据集保存到内存中,而不是磁盘
num_bands = in_ds.RasterCount
print(f'波段数为:{num_bands}')
print(f'第一张影像的投影信息是{in_ds.GetProjection()}')
out_ds = driver.Create('', columns, rows, num_bands, in_band.DataType) #这里只能设置输出结果的形状大小,波段等,现在是没有坐标意义的
# 给结果添加坐标意义
out_ds.SetProjection(in_ds.GetProjection())
geotrans[0] = min_x
geotrans[3] = max_y
out_ds.SetGeoTransform(geotrans)
'''------------------------------------------第二步:往空栅格里面写数据----------------------------------------------'''
for b in tqdm(range(1, num_bands + 1), desc='正在镶嵌真彩色近红外波段影像文件'): # 从所有影像的第一个波段开始镶嵌
for tif in tiffiles_path_list:
in_tif = gdal.Open(tif)
# 计算in_tif中左上角像元对应out_ds中的行列号
trans = gdal.Transformer(in_tif, out_ds, []) # in_tif是输入栅格,out_ds是目标栅格,第三个参数设置转换器选项,这里不需要
# success接收是否转换成功,False表示计算从源栅格到目标栅格的偏移
success, xyz = trans.TransformPoint(False, 0, 0)
x, y, z = map(int, xyz)
data = in_tif.GetRasterBand(b).ReadAsArray() #将in_tif第b波段读成数组
print(f'data:{data}')
# 找到不包含NaN值的行的索引
valid_rows_index = np.all(~np.isnan(data), axis=1)
# 提取有效的像元
valid_data = data[valid_rows_index]
print(f'valid_data:{valid_data}')
row_difference = len(data) - len(valid_data)
print(f'行差:{row_difference}')
if row_difference > 1:
row_difference = 1
print(f'子影像的有效值在空栅格中的像素坐标为第{y + row_difference}行,第{x}列')
# 将删除NaN后的有效像元写入到输出数据集中
out_ds.GetRasterBand(b).WriteArray(valid_data, x, y + row_difference)
# out_ds.GetRasterBand(b).WriteArray(data, x, y) # 将data从空栅格的第{y}行,第{x}列开始写入
del in_tif
# cutline_ds = ogr.Open(geojson)
# cutline_layer = cutline_ds.GetLayer()
input_raster = out_ds
# gdal.PushErrorHandler('CPLQuietErrorHandler')
result = gdal.Warp(
tiff_out_path,
input_raster,
format='GTiff',
cutlineDSName=geojson,
cropToCutline=True,
dstNodata=0
)
result.FlushCache()