Python批量拼接遥感影像数据

该代码示例展示了如何利用GDAL库在Python中实现多幅地理栅格图像的拼接。首先获取每张图像的边界坐标,然后计算出最大边界以确定输出图像的大小。接着,创建一个新的GeoTIFF文件作为拼接结果,并将所有输入图像的数据逐波段写入到结果文件中。最后,对输出图像设置投影和地理变换信息。

摘要生成于 C知道 ,由 DeepSeek-R1 满血版支持, 前往体验 >

直接上代码:

import os
import glob
from osgeo import gdal
from math import ceil
from tqdm import tqdm


def GetExtent(infile):
    ds = gdal.Open(infile)
    geotrans = ds.GetGeoTransform()
    xsize = ds.RasterXSize
    ysize = ds.RasterYSize
    min_x, max_y = geotrans[0], geotrans[3]
    max_x, min_y = geotrans[0] + xsize * geotrans[1], geotrans[3] + ysize * geotrans[5]
    ds = None
    return min_x, max_y, max_x, min_y

def RasterMosaic(file_list, outpath):
    Open = gdal.Open
    min_x, max_y, max_x, min_y = GetExtent(file_list[0])
    for infile in file_list:
        minx, maxy, maxx, miny = GetExtent(infile)
        min_x, min_y = min(min_x, minx), min(min_y, miny)
        max_x, max_y = max(max_x, maxx), max(max_y, maxy)

    in_ds = Open(file_list[0])
    in_band = in_ds.GetRasterBand(3) # 这里的3是数据的波段数量。根据自己的数据来

    geotrans = list(in_ds.GetGeoTransform())
    width, height = geotrans[1], geotrans[5]
    columns = ceil((max_x - min_x) / width)  # 列数
    rows = ceil((max_y - min_y) / (-height))  # 行数

    driver = gdal.GetDriverByName('GTiff')
    out_ds = driver.Create(outpath, columns, rows, 3, in_band.DataType, options=["TILED=YES", "COMPRESS=LZW", "BIGTIFF=YES"]) # 这里的3同上
    out_ds.SetProjection(in_ds.GetProjection())
    geotrans[0] = min_x  # 更正左上角坐标
    geotrans[3] = max_y
    out_ds.SetGeoTransform(geotrans)
    inv_geotrans = gdal.InvGeoTransform(geotrans)
    for in_fn in tqdm(file_list):
        in_ds = Open(in_fn)
        in_gt = in_ds.GetGeoTransform()
        offset = gdal.ApplyGeoTransform(inv_geotrans, in_gt[0], in_gt[3])
        x, y = map(int, offset)
        for i in range(3):  #每个波段都要考虑。这里的3也是,同上
            data = in_ds.GetRasterBand(i+1).ReadAsArray()
            out_ds.GetRasterBand(i+1).WriteArray(data, x, y)  # x,y是开始写入时左上角像元行列号
    del in_ds, out_ds

if __name__ == '__main__':
    image_path = "image"    # 待拼接图片路径
    result_path = "result"  # 拼接结果路径
    imageList = glob.glob(image_path + "/*.tif")
    result = os.path.join(result_path, "result.tif")
    RasterMosaic(imageList, result)
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

清纯世纪

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值