python 遙感影像分塊

在深度學習或其他一些應用中，往往需要對較大的遙感影像進行分塊輸出，并且保留原始的地理信息，以備后續(xù)使用。

本篇文章的目的即對遙感影像進行分塊輸出，按照網(wǎng)格順序進行重新命名，其中邊緣部分可能存在不足裁剪大小的情況，通過向前填補方式進行填充，故分塊結果可能導致邊緣部分存在一定重復率。

具體代碼如下：

#!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 Pathclass TifCrop:def __init__(self,infile,crop_size,save_path,repete_rate=0):'''遙感影像分塊函數(shù):param infile: 輸入tif文件:param crop_size: 分塊大小，單值或元祖，int型。eg:200表示以 200*200個像元大小的方形進行分塊，(100,200)表示以 100*200個像元大小的矩形進行分塊:param repete_size:重復率, float, 其中值的范圍為[0,1)之間, 默認值為0'''self.infile = infileself.crop_size = crop_sizeself.save_path = save_pathself.repete_rate = repete_rate# crop_size 參數(shù)判斷if not isinstance(crop_size, int):if not isinstance(crop_size, tuple):raise Exception('crop_size 輸入?yún)?shù)錯誤')else:if not (isinstance(crop_size[0], int) and isinstance(crop_size[1], int)):raise Exception('crop_size 輸入?yún)?shù)錯誤')# repete_rate 參數(shù)判斷if repete_rate >= 1 or repete_rate < 0:raise Exception('repete_rate 輸出參數(shù)錯誤')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_sizecrop_size_c = self.crop_sizerepete_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))# 循環(huán)讀取# 邊緣部分按照向前擴充原則進行提取# 當重復率較高或分塊尺寸較小時，遇到邊緣部分可能存在分割相同的情況，故用以下參數(shù)進行判斷避免該情況發(fā)生start_point = (-1,-1)for i in range(col_num):for j in range(row_num):row_s = repete_size_r * jrow_e = repete_size_r * j + crop_size_r# 是否超出邊界判斷if row_e > rows:row_s = rows - crop_size_rrow_e = rowscol_s = repete_size_c * icol_e = repete_size_c * i + crop_size_c# 是否超出邊界判斷if col_e > cols:col_s = cols - crop_size_ccol_e = colsdata_crop = data[:,row_s:row_e,col_s:col_e]# 判斷輸出內(nèi)容是否與之前存在重復情況，非完全重疊部分再進行分塊輸出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).suffixself.tif_write(data_crop,new_geotrans,out_file)def tif_write(self,data,trans,ofile):'''tif寫入:param data: 分塊后數(shù)組:param trans: 更新后的geotransform，包括六參數(shù):param ofile: 輸出全路徑:return: None'''# 數(shù)據(jù)類型獲取if 'int8' in data.dtype.name:datatype = gdal.GDT_Byteelif 'int16' in data.dtype.name:datatype = gdal.GDT_UInt16else:datatype = gdal.GDT_Float32# 輸出tif文件按照單波段或多波段劃分bands,height,width = data.shape# 創(chuàng)建文件driver = gdal.GetDriverByName("GTiff")dataset = driver.Create(ofile, int(width), int(height), int(bands), datatype)if (dataset != None):dataset.SetGeoTransform(trans) # 寫入仿射變換參數(shù)dataset.SetProjection(self.projection) # 寫入投影for i in range(bands):dataset.GetRasterBand(i + 1).WriteArray(data[i])del datasetif __name__ == '__main__':file = r'F:\input.tif'cropsize = 300savepath = r'F:\outpath'TifCrop(file,cropsize,savepath).crop_tif()

歡迎指正~

總結

以上是生活随笔為你收集整理的python——遥感影像分块的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。