python 遥感影像分块
在深度学习或其他一些应用中,往往需要对较大的遥感影像进行分块输出,并且保留原始的地理信息,以备后续使用。
本篇文章的目的即对遥感影像进行分块输出,按照网格顺序进行重新命名,其中边缘部分可能存在不足裁剪大小的情况,通过向前填补方式进行填充,故分块结果可能导致边缘部分存在一定重复率。
具体代码如下:
#!usr/bin/env python
# -*- coding: utf-8 -*-
"""
date:2022/9/2
author:甲戌_Tr
email: liu_xxxi@163.com
"""
import sys,os
from osgeo import gdal
import numpy as np
from pathlib import Path as Path
class TifCrop:
def __init__(self,infile,crop_size,save_path,repete_rate=0):
'''
遥感影像分块函数
:param infile: 输入tif文件
:param crop_size: 分块大小,单值或元祖,int型。eg:200表示以 200*200个像元大小的方形进行分块,(100,200)表示以 100*200个像元大小的矩形进行分块
:param repete_size:重复率, float, 其中值的范围为[0,1)之间, 默认值为0
'''
self.infile = infile
self.crop_size = crop_size
self.save_path = save_path
self.repete_rate = repete_rate
# crop_size 参数判断
if not isinstance(crop_size, int):
if not isinstance(crop_size, tuple):
raise Exception('crop_size 输入参数错误')
else:
if not (isinstance(crop_size[0], int) and isinstance(crop_size[1], int)):
raise Exception('crop_size 输入参数错误')
# repete_rate 参数判断
if repete_rate >= 1 or repete_rate < 0:
raise Exception('repete_rate 输出参数错误')
def crop_tif(self):
if isinstance(self.crop_size,tuple):
crop_size_r = self.crop_size[0]
crop_size_c = self.crop_size[1]
else:
crop_size_r = self.crop_size
crop_size_c = self.crop_size
repete_size_r = int(crop_size_r * (1 - self.repete_rate))
repete_size_c = int(crop_size_c * (1 - self.repete_rate))
ds = gdal.Open(self.infile)
data = ds.ReadAsArray()
geotrans = ds.GetGeoTransform()
self.projection = ds.GetProjection()
# 将单波段影像添加一个维度
if len(data.shape) == 2:
data = np.array([data])
channel, rows, cols = data.shape
# 向上取整
col_num = int(np.ceil(cols / repete_size_c))
row_num = int(np.ceil(rows / repete_size_r))
# 循环读取
# 边缘部分按照向前扩充原则进行提取
# 当重复率较高或分块尺寸较小时,遇到边缘部分可能存在分割相同的情况,故用以下参数进行判断避免该情况发生
start_point = (-1,-1)
for i in range(col_num):
for j in range(row_num):
row_s = repete_size_r * j
row_e = repete_size_r * j + crop_size_r
# 是否超出边界判断
if row_e > rows:
row_s = rows - crop_size_r
row_e = rows
col_s = repete_size_c * i
col_e = repete_size_c * i + crop_size_c
# 是否超出边界判断
if col_e > cols:
col_s = cols - crop_size_c
col_e = cols
data_crop = data[:,row_s:row_e,col_s:col_e]
# 判断输出内容是否与之前存在重复情况,非完全重叠部分再进行分块输出
if (row_s,col_s) != start_point:
start_point = (row_s,col_s)
# 地理信息存放
new_geotrans = (geotrans[0] + geotrans[1] * col_s, geotrans[1], geotrans[2], geotrans[3] + geotrans[5] * row_s,geotrans[4],geotrans[5])
# 输出文件名称
out_file = self.save_path + os.sep + Path(self.infile).stem + '_' + str(j) + '_' + str(i) + Path(self.infile).suffix
self.tif_write(data_crop,new_geotrans,out_file)
def tif_write(self,data,trans,ofile):
'''
tif写入
:param data: 分块后数组
:param trans: 更新后的geotransform,包括六参数
:param ofile: 输出全路径
:return: None
'''
# 数据类型获取
if 'int8' in data.dtype.name:
datatype = gdal.GDT_Byte
elif 'int16' in data.dtype.name:
datatype = gdal.GDT_UInt16
else:
datatype = gdal.GDT_Float32
# 输出tif文件按照单波段或多波段划分
bands,height,width = data.shape
# 创建文件
driver = gdal.GetDriverByName("GTiff")
dataset = driver.Create(ofile, int(width), int(height), int(bands), datatype)
if (dataset != None):
dataset.SetGeoTransform(trans) # 写入仿射变换参数
dataset.SetProjection(self.projection) # 写入投影
for i in range(bands):
dataset.GetRasterBand(i + 1).WriteArray(data[i])
del dataset
if __name__ == '__main__':
file = r'F:\input.tif'
cropsize = 300
savepath = r'F:\outpath'
TifCrop(file,cropsize,savepath).crop_tif()
欢迎指正~