使用GDAL的MEM内存文件保存临时文件

本文探讨了在使用GDAL进行图像处理时如何避免使用临时文件,介绍了VRT虚拟文件和MEM内存文件的应用。通过实例展示了如何利用MEM文件进行数据处理,以此减少磁盘I/O操作并提高算法效率。

在使用GDAL编写算法的时候,经常会将计算的中间结果存在一个临时的图像文件中,然后使用完再将其删除,如果临时文件就一个的话,创建一个也无所谓,但是当一个复杂的算法中可能会出现很多个临时文件的时候(我在编写Hariss角点自动匹配算法的时候有4个临时文件),这种情况下总觉得临时文件很不爽,此外第一个不爽的地方;第二个图像太大的时候,临时文件也会占用很大的空间,假如空间不足或者给定的临时文件路径不可写等问题会让人头疼;第三,在创建临时文件的读写上会耗用比较多的时间,尤其在磁盘的IO时,耗时比较多。

基于上面的几个原因,总想不用临时文件,发现GDAL里面提供了一个VRT的虚拟文件(参考http://www.gdal.org/gdal_vrttut.html),查看GDAL的App中的gdalenhance.cpp文件中,发现VRT使用起来比较复杂,或者我不太熟悉而已,过段时间仔细研究一下这个东西。除了VRT,GDAL中还有一个叫MEM的内存文件,简介参考这个网页:http://gdal.org/frmt_mem.html。关于MEM文件的使用,和普通的文件使用一样,参考下面的代码:

#include "gdal_priv.h"

int main()
{
    GDALAllRegister();

	//打开原始数据
    GDALDataset  *poSrcDS;
    poSrcDS = (GDALDataset *) GDALOpen( "C:\\test.img", GA_ReadOnly );
    if( poSrcDS == NULL )
        return -1;
	
	//获取MEM的文件驱动
    GDALDriver *poDriver;
    poDriver = GetGDALDriverManager()->GetDriverByName("MEM");
    if( poDriver == NULL )
        return -1;
        
	//这里或者使用Create函数创建一个MEM文件,文件路径给个空串就可以了
    GDALDataset *poDstDS;
    poDstDS = poDriver->CreateCopy( "", poSrcDS, FALSE, NULL, NULL, NULL );
	
	{
	//do something
	}
	
    //使用完直接Close掉就好了,不用删除临时文件了
    GDALClose( (GDALDatasetH) poDstDS );
    
    GDALClose( (GDALDatasetH) poSrcDS );
	
	return 0;
}


from osgeo import gdal, osr import numpy as np import os import tempfile import shutil def merge_tiff_files_with_color(input_files, output_file, resample_alg='nearest', compress='LZW'): """ 合并多个TIFF文件(具有不同分辨率和波段数)为一个TIFF文件,保留原始颜色并生成12级金字塔 参数: input_files: 输入文件路径列表 output_file: 输出文件路径 resample_alg: 重采样方法,默认为'nearest'(最邻近) compress: 压缩方式,默认为'LZW' """ if not input_files: raise ValueError("输入文件列表为空") # 创建临时目录处理中间文件 temp_dir = tempfile.mkdtemp() try: # 步骤1: 分析输入文件并处理颜色表 rgb_input_files = [] color_table = None same_color_table = True first_file = True for i, file in enumerate(input_files): ds = gdal.Open(file) if ds is None: raise Exception(f"无法打开文件:{file}") # 检查波段数和颜色表 band = ds.GetRasterBand(1) ct = band.GetColorTable() num_bands = ds.RasterCount # 如果是单波段且有颜色表 if num_bands == 1 and ct is not None: if first_file: color_table = ct.Clone() # 保存第一个文件的颜色表 first_file = False else: # 检查颜色表是否一致 if not same_color_table or not ct.IsSame(color_table): same_color_table = False # 转换为RGB保留颜色 rgb_file = os.path.join(temp_dir, f"rgb_{i}.tif") gdal.Translate(rgb_file, ds, bandList=[1, 1, 1], outputType=gdal.GDT_Byte, creationOptions=["PHOTOMETRIC=RGB"]) rgb_input_files.append(rgb_file) else: # 多波段文件直接使用(假设已经是RGB) rgb_input_files.append(file) ds = None # 步骤2: 计算所有输入文件中的最高分辨率(最小像元大小)和目标范围(并集) min_pixel_size_x = None min_pixel_size_y = None all_extents = [] # 每个元素为(minX, maxX, minY, maxY) for file in rgb_input_files: ds = gdal.Open(file) geotransform = ds.GetGeoTransform() if geotransform is None: ds = None raise Exception(f"文件{file}没有地理变换信息。") # 计算分辨率(绝对值) pixel_size_x = abs(geotransform[1]) pixel_size_y = abs(geotransform[5]) # 更新最小分辨率 if min_pixel_size_x is None or pixel_size_x < min_pixel_size_x: min_pixel_size_x = pixel_size_x if min_pixel_size_y is None or pixel_size_y < min_pixel_size_y: min_pixel_size_y = pixel_size_y # 计算当前影像的范围 minX = geotransform[0] maxX = minX + geotransform[1] * ds.RasterXSize minY = geotransform[3] + geotransform[5] * ds.RasterYSize maxY = geotransform[3] if geotransform[5] > 0: # 处理罕见情况 minY, maxY = maxY, minY all_extents.append((minX, maxX, minY, maxY)) ds = None # 计算所有范围的最小外接矩形(并集) all_minX = min([extent[0] for extent in all_extents]) all_maxX = max([extent[1] for extent in all_extents]) all_minY = min([extent[2] for extent in all_extents]) all_maxY = max([extent[3] for extent in all_extents]) # 步骤3: 设置Warp选项 creationOptions = ['TILED=YES', 'BIGTIFF=YES'] if compress is not None: creationOptions.append(f'COMPRESS={compress}') warp_options = gdal.WarpOptions( format='GTiff', outputBounds=(all_minX, all_minY, all_maxX, all_maxY), xRes=min_pixel_size_x, yRes=min_pixel_size_y, resampleAlg=resample_alg, creationOptions=creationOptions, dstSRS=None, multithread=True, targetAlignedPixels=True, ) # 执行Warp合并 ds = gdal.Warp(output_file, rgb_input_files, options=warp_options) # 如果所有输入文件有相同的颜色表,应用到输出文件 if same_color_table and color_table is not None: output_band = ds.GetRasterBand(1) output_band.SetColorTable(color_table) ds = None # 关闭文件,写入磁盘 # 步骤4: 为合并后的TIFF生成12级金字塔 print("开始生成12级金字塔...") ds = gdal.Open(output_file, gdal.GA_Update) if ds is None: raise Exception(f"无法打开输出文件以生成金字塔: {output_file}") # 设置金字塔配置选项 gdal.SetConfigOption('BIGTIFF_OVERVIEW', 'IF_NEEDED') # 支持大文件 if compress: gdal.SetConfigOption('COMPRESS_OVERVIEW', compress) # 应用压缩 # 定义12级金字塔层级 (2^1 到 2^12) overview_levels = [2 ** i for i in range(1, 13)] # 生成金字塔(使用最近邻采样法) ds.BuildOverviews("NEAREST", overview_levels) # 关闭数据集 ds = None print("金字塔生成完成!") return output_file finally: # 清理临时目录 shutil.rmtree(temp_dir, ignore_errors=True) # 使用示例 if __name__ == "__main__": input_files = [ r"F:\曹妃甸水域正射\二维\冯晨2\image_pyramid.tif", r"F:\曹妃甸水域正射\25-08-07-大学城-水域2\8-7-管委会-image_pyramid.tif", r"F:\曹妃甸水域正射\二维\B河\B河_image_pyramid.tif", r"F:\曹妃甸水域正射\25-08-01-大学城-水域1\image_pyramid.tif", r"F:\曹妃甸水域正射\二维\冯晨3\image_pyramid.tif", r"F:\曹妃甸水域正射\二维\主河道\主河道_image_pyramid.tif" ] output_file = r"E:\游戏\merged_output_with_color.tif" merge_tiff_files_with_color(input_files, output_file) print(f"合并完成,结果保存至: {output_file}") 为合并后的Tif文件去除黑边,给我修改后的完整代码
最新发布
08-20
评论 15
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值