赞
踩
在处理tif/tiff大文件的过程中,有一些业务需求将大图进行裁剪,arcgis的批量化裁剪操作比较繁琐,写个脚本记录并解决问题。
采用gdal实现。
// An highlighted block from osgeo import gdal import os from tqdm import tqdm # 打开大图 ds = gdal.Open(r'xxx.tif') # 获取大图的宽和高 width = ds.RasterXSize height = ds.RasterYSize # 设置裁剪大小 cut_size = 1024 #512 # 计算横向和纵向的裁剪数量 x_count = width // cut_size y_count = height // cut_size # 遍历所有裁剪后的小图 for y in tqdm(range(y_count)): for x in range(x_count): # 计算当前小图的左上角坐标 xoff = x * cut_size yoff = y * cut_size # 计算当前小图的宽和高 xsize = cut_size ysize = cut_size # 如果最后一行或最后一列,宽或高可能不足cut_size,需要重新计算宽和高 if x == x_count - 1: xsize = width - xoff if y == y_count - 1: ysize = height - yoff # 创建小图文件名,根据实际需要进行修改 out_tif = f'cut_image_{x}_{y}.tif' tif_path=os.path.join(r"save_path",out_tif) # 创建输出数据集 driver = gdal.GetDriverByName('GTiff') out_ds = driver.Create(tif_path, xsize, ysize, ds.RasterCount, ds.GetRasterBand(1).DataType) # 设置输出数据集的地理信息 out_ds.SetGeoTransform((ds.GetGeoTransform()[0] + xoff * ds.GetGeoTransform()[1], ds.GetGeoTransform()[1], 0, ds.GetGeoTransform()[3] + yoff * ds.GetGeoTransform()[5], 0, ds.GetGeoTransform()[5])) # 逐波段读取数据,并写入输出数据集 for band_num in range(ds.RasterCount): band = ds.GetRasterBand(band_num + 1) data = band.ReadAsArray(xoff, yoff, xsize, ysize) out_band = out_ds.GetRasterBand(band_num + 1) out_band.WriteArray(data) # 释放内存 out_ds.FlushCache() del out_ds
注释完整,修改路径即可。
Copyright © 2003-2013 www.wpsshop.cn 版权所有,并保留所有权利。