57 lines
2.2 KiB
Python
57 lines
2.2 KiB
Python
def resampling_by_scale(input_path, target_file, refer_img_path):
|
|
"""
|
|
按照缩放比例对影像重采样
|
|
:param input_path: GDAL地理数据路径
|
|
:param target_file: 输出影像
|
|
:param refer_img_path:参考影像
|
|
:return: True or False
|
|
"""
|
|
ref_dataset = gdal.Open(refer_img_path)
|
|
ref_cols = ref_dataset.RasterXSize # 列数
|
|
ref_rows = ref_dataset.RasterYSize # 行数
|
|
|
|
dataset = gdal.Open(input_path)
|
|
if dataset is None:
|
|
logger.error('resampling_by_scale:dataset is None!')
|
|
return False
|
|
|
|
band_count = dataset.RasterCount # 波段数
|
|
if (band_count == 0) or (target_file == ""):
|
|
logger.error("resampling_by_scale:Parameters of the abnormal!")
|
|
return False
|
|
|
|
cols = dataset.RasterXSize # 列数
|
|
scale = ref_cols/cols
|
|
|
|
# rows = dataset.RasterYSize # 行数
|
|
# cols = int(cols * scale) # 计算新的行列数
|
|
# rows = int(rows * scale)
|
|
cols = ref_cols
|
|
rows = ref_rows
|
|
|
|
geotrans = list(dataset.GetGeoTransform())
|
|
geotrans[1] = geotrans[1] / scale # 像元宽度变为原来的scale倍
|
|
geotrans[5] = geotrans[5] / scale # 像元高度变为原来的scale倍
|
|
|
|
if os.path.exists(target_file) and os.path.isfile(target_file): # 如果已存在同名影像
|
|
os.remove(target_file) # 则删除之
|
|
if not os.path.exists(os.path.split(target_file)[0]):
|
|
os.makedirs(os.path.split(target_file)[0])
|
|
|
|
band1 = dataset.GetRasterBand(1)
|
|
data_type = band1.DataType
|
|
target = dataset.GetDriver().Create(target_file, xsize=cols, ysize=rows, bands=band_count,
|
|
eType=data_type)
|
|
target.SetProjection(dataset.GetProjection()) # 设置投影坐标
|
|
target.SetGeoTransform(geotrans) # 设置地理变换参数
|
|
total = band_count + 1
|
|
for index in range(1, total):
|
|
# 读取波段数据
|
|
data = dataset.GetRasterBand(index).ReadAsArray(buf_xsize=cols, buf_ysize=rows)
|
|
out_band = target.GetRasterBand(index)
|
|
|
|
no_data_value = dataset.GetRasterBand(index).GetNoDataValue() # 获取没有数据的点
|
|
if not (no_data_value is None):
|
|
out_band.SetNoDataValue(no_data_value)
|
|
|
|
out_band.WriteArray(data) # |