GDAL影像合并

本文介绍如何利用GDAL库进行多幅影像的高效合并,重点讲解自己构建掩膜的方法。

摘要生成于 C知道 ,由 DeepSeek-R1 满血版支持, 前往体验 >

// merge.cpp : 定义控制台应用程序的入口点。
//

#include "stdafx.h"

#include "../3dparty/include/gdal_priv.h"
#include "../3dparty/include/gdal_alg.h"
#include "../3dparty/include/gdalwarper.h"

#include <cassert>

#ifdef _DEBUG
#pragma comment(lib, "../3dparty/lib/gdal_i.lib")
#else
#pragma comment(lib, "../3dparty/lib/gdal_i.lib")
#endif

#define SIZEOFARRAY(array) sizeof(array)/sizeof(array[0])

static double dfMinX=0.0, dfMinY=0.0, dfMaxX=0.0, dfMaxY=0.0;
static double dfXRes=0.0, dfYRes=0.0;

static int bTargetAlignedPixels = FALSE;
static int nForcePixels=0, nForceLines=0, bQuiet = FALSE;
static int bEnableDstAlpha = FALSE, bEnableSrcAlpha = FALSE;
static int bVRT = FALSE;

static char *pszDstFileName = "../data/out.tif";
static char *apszSrcFiles[] = 
{
	/*"../data/H50G038013DOM.tif",
	"../data/H50G038014DOM.tif",
	"../data/H50G038015DOM.tif",
	"../data/H50G039012DOM.tif",
	"../data/H50G039013DOM.tif",
	"../data/H50G039014DOM.tif",*/
	"../data/H50G039015DOM.tif",
	"../data/H50G039016DOM.tif"
};

GDALDataType eWorkingType = GDT_Byte;

static GDALResampleAlg eResampleAlg = GRA_NearestNeighbour;
static GDALDatasetH GDALWarpCreateOutput( char **papszSrcFiles, 
										 int nFilesCount, 
										 const char *pszFilename, 
										 const char *pszFormat, 
										 GDALDataType eDT );

static void CreateMask(char* const* paFilePath, const int nFilesCount);

int _tmain(int argc, _TCHAR* argv[])
{
	GDALAllRegister();

	GDALDatasetH hSrcDS = NULL;
	GDALDatasetH hDstDS = NULL;

	int nFilesCount = SIZEOFARRAY(apszSrcFiles);
	hDstDS = GDALWarpCreateOutput(apszSrcFiles, nFilesCount,
		pszDstFileName, "GTiff", GDT_Byte);

	char **papszSrcFiles = apszSrcFiles;
	CreateMask(papszSrcFiles, nFilesCount);


	for( int iSrc = 0; iSrc < nFilesCount; iSrc++ )
	{
		GDALDatasetH hSrcDS;

		hSrcDS = GDALOpen( apszSrcFiles[iSrc], GA_ReadOnly );
		GDALDataset *pSrcDS = (GDALDataset *)hSrcDS;
		assert(pSrcDS->GetRasterBand(1)->GetMaskBand() != NULL);

		if( hSrcDS == NULL )
			exit( 2 );

		if ( GDALGetRasterCount(hSrcDS) == 0 )
		{
			fprintf(stderr, "Input file %s has no raster bands.\n", apszSrcFiles[iSrc] );
			exit( 1 );
		}

		if( !bQuiet )
			printf( "Processing input file %s.\n", apszSrcFiles[iSrc] );

		if ( eResampleAlg != GRA_NearestNeighbour &&
			GDALGetRasterColorTable(GDALGetRasterBand(hSrcDS, 1)) != NULL)
		{
			if( !bQuiet )
				fprintf( stderr, "Warning: Input file %s has a color table, which will likely lead to "
				"bad results when using a resampling method other than "
				"nearest neighbour. Converting the dataset prior to 24/32 bit "
				"is advised.\n", apszSrcFiles[iSrc] );
		}

		if( GDALGetRasterColorInterpretation( 
			GDALGetRasterBand(hSrcDS,GDALGetRasterCount(hSrcDS)) ) 
			== GCI_AlphaBand 
			&& !bEnableSrcAlpha )
		{
			bEnableSrcAlpha = TRUE;
			if( !bQuiet )
				printf( "Using band %d of source image as alpha.\n", 
				GDALGetRasterCount(hSrcDS) );
		}
		void *hTransformArg = NULL;

		hTransformArg = GDALCreateGenImgProjTransform
### 使用 GDAL 实现遥感影像拼接 #### 数据准备与环境配置 在使用 GDAL 进行遥感影像拼接之前,需确保安装了支持 Python 的 GDAL 库以及 NumPy 等依赖库。对于 JPG 格式的无人机数据,由于其可能缺乏地理坐标信息,建议先通过 EXIF 工具提取 GPS 坐标并将其写入到 GeoTIFF 文件中[^1]。 #### 读取影像数据 GDAL 提供了一种高效的方式用于读取栅格数据,并可将其转换为数组形式以便进一步操作。以下是基本的代码框架: ```python from osgeo import gdal, ogr, osr import numpy as np def read_image(file_path): dataset = gdal.Open(file_path) if not dataset: raise Exception(f"Failed to open {file_path}") geotransform = dataset.GetGeoTransform() projection = dataset.GetProjection() bands = [] for i in range(dataset.RasterCount): band = dataset.GetRasterBand(i + 1).ReadAsArray() # Bands are indexed from 1 bands.append(band) return { 'geotransform': geotransform, 'projection': projection, 'bands': bands } ``` 此函数返回影像的几何变换参数、投影信息及其波段数据[^2]。 #### 横向合并矩阵 当多个影像具有相同的行列数时,可以通过简单的 Numpy 数组连接完成水平方向上的拼接。如果涉及不同分辨率或大小,则需要额外考虑插值方法。 ```python def merge_horizontally(image_list): merged_bands = [] first_geotransform = image_list[0]['geotransform'] first_projection = image_list[0]['projection'] for b_idx in range(len(image_list[0]['bands'])): horizontal_band = np.hstack([img['bands'][b_idx] for img in image_list]) merged_bands.append(horizontal_band) new_width = sum(img['bands'][0].shape[1] for img in image_list) updated_geo_transform = list(first_geotransform) updated_geo_transform[0] -= (new_width * updated_geo_transform[1]) / 2 result = {'geotransform': tuple(updated_geo_transform), 'projection': first_projection, 'bands': merged_bands} return result ``` 上述逻辑假设所有输入图片共享一致的空间参考系和像素尺寸[^4]。 #### 计算仿射变换参数 为了使最终合成图保持正确的地理位置关系,在每一步都需要更新相应的仿射变换矩阵(即 `GetGeoTransform()` 返回的结果)。这通常涉及到调整原点位置(x,y) 和旋转角度等要素。 #### 输出新影像文件 最后一步是把处理后的多维数组保存回磁盘作为新的 TIFF 文件。 ```python def write_to_file(output_filename, data_dict): driver = gdal.GetDriverByName('GTiff') out_ds = driver.Create( output_filename, data_dict['bands'][0].shape[1], # Width data_dict['bands'][0].shape[0], # Height len(data_dict['bands']), # Number of bands gdal.GDT_UInt16 # Data type ) out_ds.SetGeoTransform(data_dict['geotransform']) out_ds.SetProjection(data_dict['projection']) for idx, band_data in enumerate(data_dict['bands']): out_band = out_ds.GetRasterBand(idx + 1) out_band.WriteArray(band_data) del out_ds ``` 以上过程涵盖了从原始数据加载直至生成完整镶嵌产品的整个工作流[^3]。
评论 9
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值