python：多波段遥感影像镶嵌与裁剪

Guyuecc0906

已于 2024-06-01 21:37:39 修改

阅读量743

点赞数 17

文章标签： python 开发语言

于 2024-04-15 14:39:01 首次发布

本文链接：https://blog.csdn.net/weixin_57058676/article/details/137779581

版权

本文介绍了如何使用Python和GDAL库在内存中处理多幅遥感影像，通过创建空的镶嵌结果栅格、写入数据和裁剪功能，实现影像的无缝拼接。详细步骤包括获取影像范围、确定镶嵌尺寸、写入数据并处理缺失值。

摘要由CSDN通过智能技术生成

两张待镶嵌影像

镶嵌和裁剪后的影像
（这里只用了两幅影像镶嵌，所以左边是平整的，完整来说需要六幅影像）

遥感影像镶嵌

影像镶嵌是指将多幅相邻的遥感影像拼接成一个连续的、无缝的大范围影像的过程。

在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()

Guyuecc0906

关注

17
点赞
踩
13

收藏

觉得还不错? 一键收藏
打赏
1
评论
python：多波段遥感影像镶嵌与裁剪

首先，我们需要在内存中创建一个“空”的镶嵌结果栅格，来接收我们的遥感数据。这个栅格有多大、几个波段、数据类型是什么都需要设置，以及给它添加坐标意义（投影、地理转换信息）。接下来，我们需要将每一幅输入影像的数据写入到我们创建的“空”栅格中。这个过程需要对每一幅影像进行遍历，并计算其在镶嵌结果栅格中的位置。影像镶嵌是指将多幅相邻的遥感影像拼接成一个连续的、无缝的大范围影像的过程。tiffiles_path_list是一个有多个遥感影像路径的列表。在Python中，我们可以使用GDAL的。
复制链接

扫一扫