python Arcpy/GDAL对遥感影像进行重采样
【代码】python arcpy对遥感影像进行重采样。
·
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}")
更多推荐
已为社区贡献1条内容
所有评论(0)