1. Python Arcpy对遥感影像进行重采样

使用Arcpy对影像进行重采样,需要使用对应的python版本为2.7,下面为对单个波段的数据重采样示例:

# -*- coding: utf-8 -*-
import os
import arcpy

# 设置父文件夹路径
parent_folder = ur"D:\Sentinel2\10m"
# 设置输出文件夹路径
output_folder = ur"D:\Sentinel2\100m"
# 设置重采样分辨率
target_resolution = "100"
if os.path.exists(parent_folder):
    # print("Directory exists.")
    # print("Contents:", os.listdir(parent_folder))
    # 遍历父文件夹下的所有子文件夹
    for subdir, dirs, files in os.walk(parent_folder):
        for file in files:
            # 检查文件扩展名是否为.tif
            if file.endswith(".tif") or file.endswith(".tiff"):
                # 构建输入和输出路径
                input_raster = os.path.join(subdir, file)
                # 获取相对路径
                relative_path = os.path.relpath(input_raster, parent_folder)
                # 构建输出路径,保留相对路径结构
                output_raster = os.path.join(output_folder, relative_path)
                # 获取输入影像的空间参考
                desc = arcpy.Describe(input_raster)
                coordinate_system = desc.spatialReference
                # 如果输出文件夹中的子目录不存在,创建它
                output_subdir = os.path.dirname(output_raster)
                if not os.path.exists(output_subdir):
                    os.makedirs(output_subdir)
                # 设置环境参数,包括地理转换和分辨率
                arcpy.env.outputCoordinateSystem = coordinate_system
                arcpy.env.cellSize = target_resolution
                # 使用Resample_management工具进行重采样
                arcpy.Resample_management(
                    input_raster,
                    output_raster,
                    target_resolution,
                    "BILINEAR",  # 选择合适的插值方法,例如"BILINEAR"或"CUBIC"
                )
                print u"重采样完成:%s" % output_raster
else:
    print "Directory does not exist."

2. Python GDAL对遥感影像进行重采样

最近数据处理又用到重采样,和之前不同的是,这次先把数据进行了波段合成,参见此博客方法Python 根据经纬度批量裁剪TIFF影像 & ENVI 波段合成,一年的数据存放在同一个tiff中,这时对其进行重采样,使用GDAL实现,本次所用python版本为3.11

import os
from osgeo import gdal

def resample_tiff(input_file, output_file, x_res, y_res):
    # 打开输入文件
    input_dataset = gdal.Open(input_file)
    if not input_dataset:
        raise Exception(f"无法打开文件 {input_file}")

    # 获取输入文件的投影和仿射变换参数
    projection = input_dataset.GetProjection()
    geotransform = input_dataset.GetGeoTransform()

    # 创建输出文件
    driver = gdal.GetDriverByName('GTiff')
    output_dataset = driver.Create(output_file,
                                   xsize=int(input_dataset.RasterXSize * geotransform[1] / x_res),
                                   ysize=int(input_dataset.RasterYSize * (-geotransform[5]) / y_res),   #要获取y方向的大小,直接用geotransform[5]计算得到的为负值,因此需要加一个“-”。
                                   bands=input_dataset.RasterCount,
                                   eType=input_dataset.GetRasterBand(1).DataType)

    # 设置投影和仿射变换参数
    output_dataset.SetProjection(projection)
    new_geotransform = (geotransform[0], x_res, geotransform[2], geotransform[3], geotransform[4], -y_res)
    output_dataset.SetGeoTransform(new_geotransform)

    # 使用重采样方法,本次选择双线性内插
    gdal.ReprojectImage(input_dataset, output_dataset, projection, projection, gdal.GRA_Bilinear)

    # 清理资源
    del output_dataset
    del input_dataset

if __name__ == "__main__":
    input_file = r'E:\test.tif'
    # 设置重采样后的分辨率大小
    resolutions = [20, 40, 60, 80, 100]

    for res in resolutions:
        base, ext = os.path.splitext(input_file)
        output_file = f'{base}_{res}m{ext}'
        resample_tiff(input_file, output_file, res, res)
        print(f"已生成 {output_file}")

Logo

腾讯云面向开发者汇聚海量精品云计算使用和开发经验,营造开放的云计算技术生态圈。

更多推荐