GDAL重采样,可以通过写文件时改变图像尺寸和geo_transformes的分辨率信息实现。核心代码示例如下:
in_ds = gdal.Open(fi, gdal.GA_ReadOnly)
geotrans = in_ds.GetGeoTransform()
geotrans = list(geotrans)
geoproj = in_ds.GetProjection()
# 更新为采样后分辨率
geotrans[1] = tar_resolution #
geotrans[5] = -tar_resolution #
in_band = in_ds.GetRasterBand(2)
xsize = in_band.XSize
ysize = in_band.YSize
# 重新计算尺寸
ysize_t = int(ysize * img_resolution / tar_resolution)
xsize_t = int(xsize * img_resolution / tar_resolution)
# 按新的尺寸写
out_ds = in_ds.GetDriver().Create(fo_, xsize_t, ysize_t, 1,
in_band.DataType) # 创建一个构建重采样影像的句柄
out_ds.SetProjection(in_ds.GetProjection()) # 设置投影信息
# 更新geotrans
geotrans = tuple(geotrans)
GDAL裁切先计算索引,通过索引获取裁切后的图像,同时重新计算geo_transformes里的左上角点坐标。核心代码示例如下:
xsize_t = in_band.XSize
ysize_t = in_band.YSize
# 计算偏移
y_offset = ysize_t - height
x_offset = xsize_t - width
data = in_band.ReadAsArray()
if pro_mode == "r_b": # 假设裁剪右下部分
# 更新左上角点坐标
geotrans[3] = geotrans[3] + geotrans[5] * y_offset
geotrans[0] = geotrans[0] + geotrans[2] * x_offset
data = data[y_offset:, x_offset:]
完整代码示例:
from os import makedirs, remove
from os.path import exists, join
from osgeo import gdal
import numpy as np
def get_params(img_info, img_config):
return img_info["indir"], img_info["outdir"], img_info["img_list"], img_info["img_resolution"], \
img_info["pro_mode"], img_config["height"], img_config["width"], img_config["channel"], img_config["tar_resolution"]
def large2small(img_info, img_config):
indir, outdir, img_list, img_resolutions, pro_modes, height, width, channel, tar_resolution = get_params(img_info, img_config)
for f, img_resolution, pro_mode in zip(img_list, img_resolutions, pro_modes):
fi = join(indir, f)
fo_ = join(outdir, "_" + f)
fo = join(outdir, f)
# resize image
in_ds = gdal.Open(fi, gdal.GA_ReadOnly)
geotrans = in_ds.GetGeoTransform()
geotrans = list(geotrans)
geoproj = in_ds.GetProjection()
geotrans[1] = tar_resolution #
geotrans[5] = -tar_resolution #
in_band = in_ds.GetRasterBand(1)
xsize = in_band.XSize
ysize = in_band.YSize
ysize_t = int(ysize * img_resolution / tar_resolution)
xsize_t = int(xsize * img_resolution / tar_resolution)
if exists(fo_):
remove(fo_)
out_ds = in_ds.GetDriver().Create(fo_, xsize_t, ysize_t, 1,
in_band.DataType) # 创建一个构建重采样影像的句柄
out_ds.SetProjection(in_ds.GetProjection()) # 设置投影信息
geotrans = tuple(geotrans)
out_ds.SetGeoTransform(geotrans) # 设置地理变换信息
data = np.empty((ysize_t, xsize_t), np.uint8) # 设置一个与重采样影像行列号相等的矩阵去接受读取所得的像元值
in_band.ReadAsArray(buf_obj=data)
out_band = out_ds.GetRasterBand(1)
out_band.WriteArray(data)
out_band.FlushCache()
out_band.ComputeStatistics(False) # 计算统计信息
# out_ds.BuildOverviews('average', [1, 2, 4, 8, 16, 32]) # 构建金字塔
del in_ds # 删除句柄
del out_ds
# reload and clip
in_ds = gdal.Open(fo_, gdal.GA_ReadOnly)
geotrans = in_ds.GetGeoTransform()
geotrans = list(geotrans)
geoproj = in_ds.GetProjection()
in_band = in_ds.GetRasterBand(1)
xsize_t = in_band.XSize
ysize_t = in_band.YSize
y_offset = ysize_t - height
x_offset = xsize_t - width
data = in_band.ReadAsArray()
#
# pro_mode = "l_u"
if pro_mode == "l_b":
geotrans[3] = geotrans[3] + geotrans[5] * y_offset
data = data[y_offset:, :width]
elif pro_mode == "l_u":
data = data[:height, :width]
elif pro_mode == "r_u":
geotrans[0] = geotrans[0] + geotrans[2] * x_offset
data = data[:height, x_offset:]
elif pro_mode == "r_b":
geotrans[3] = geotrans[3] + geotrans[5] * y_offset
geotrans[0] = geotrans[0] + geotrans[2] * x_offset
data = data[y_offset:, x_offset:]
if exists(fo):
remove(fo)
options = ['COMPRESS=LZW']
out_ds = in_ds.GetDriver().Create(fo, width, height, 1,
in_band.DataType, options=options) # 创建一个构建重采样影像的句柄
out_ds.SetProjection(in_ds.GetProjection()) # 设置投影信息
geotrans = tuple(geotrans)
out_ds.SetGeoTransform(geotrans) # 设置地理变换信息
# data = np.empty((height, width), np.uint8) # 设置一个与重采样影像行列号相等的矩阵去接受读取所得的像元值
# in_band.ReadAsArray(buf_obj=data)
out_band = out_ds.GetRasterBand(1)
out_band.WriteArray(data)
out_band.FlushCache()
out_band.ComputeStatistics(False) # 计算统计信息
out_ds.BuildOverviews('average', [1, 2, 4, 8, 16, 32]) # 构建金字塔
del in_ds # 删除句柄
del out_ds
if exists(fo_):
remove(fo_)
本文内容由网友自发贡献,版权归原作者所有,本站不承担相应法律责任。如您发现有涉嫌抄袭侵权的内容,请联系:hwhale#tublm.com(使用前将#替换为@)