GDAL坐标转换六参的使用方法 .

本文详细介绍了二维坐标转换中的仿射变换原理及其在GDAL库中的实现方式,包括坐标平移、旋转和缩放的具体操作,并通过示例代码展示了如何在实际项目中应用这些概念。

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

GDAL六参坐标转换是一种二维坐标转换的参数,常在GDALDataset 类中的CPLErr SetGeoTransform( double* padfTransform )使用;下面让我们先来谈论一下二维的仿射变换。

仿射变换变换公式:

展开形式:

x’,y’为目标坐标,x,y为原始坐标,dx,dy为平移参数。A为存储用于绕(dx,dy)旋转,x,y方向的拉伸比例参数。

GDAL坐标转换六参的对应:

double *dfGeoTransform = new double[6];

dfGeoTransform[0] = dx;

dfGeoTransform[1] = a11;

dfGeoTransform[2] = a12;

dfGeoTransform[3] = dy;

dfGeoTransform[4] = a21;

dfGeoTransform[5] = a22;

初始化时,A2×2的单位矩阵,dx = dy = 0.0;

 即dfGeoTransform[] = {0.0,1.0,0.0,0.0,0.0,1.0};

坐标平移:

dx = dx + xtranslation;

dy = dy + xtranslation;

坐标旋转:

坐标缩放:

下面来个示例:

#include "ogrsf_frmts.h"
#include "stdio.h"
#include <iostream>
#include <string>
using namespace std;

// 定义坐标转换器
class GeoTransform
{
public:
	GeoTransform()
	{
		// 初始化六参数
		dfGeoTransform[0] = 0.0;
		dfGeoTransform[1] = 1.0;
		dfGeoTransform[2] = 0.0;
		dfGeoTransform[3] = 0.0;
		dfGeoTransform[4] = 0.0;
		dfGeoTransform[5] = 1.0;
	}

	~GeoTransform()
	{

	}

	// 设置平移
	void SetTranslate(const double dx, const double dy)
	{
		dfGeoTransform[0] += dx;
		dfGeoTransform[3] += dy;
	}

	// 设置旋转角度
	void SetRotate(const double dangle)
	{
		dfGeoTransform[1] = dfGeoTransform[1] * cos(dangle) 
			+ dfGeoTransform[2] * sin(dangle);
		dfGeoTransform[2] = - dfGeoTransform[1] * sin(dangle) 
			+ dfGeoTransform[2] * cos(dangle);
		dfGeoTransform[4] = dfGeoTransform[4] * cos(dangle) 
			+ dfGeoTransform[5] * sin(dangle);
		dfGeoTransform[5] = - dfGeoTransform[4] * sin(dangle) 
			+ dfGeoTransform[5] * cos(dangle);
	}
	// 设置x,y方向缩放比例
	void SetScale(const double xScale, const double yScale)
	{
		dfGeoTransform[1] = dfGeoTransform[1] * xScale;
		dfGeoTransform[2] = dfGeoTransform[2] * yScale;
		dfGeoTransform[4] = dfGeoTransform[4] * xScale;
		dfGeoTransform[5] = dfGeoTransform[5] * yScale;
	}

	// 转换坐标
	OGRPoint Transform(const double dx, const double dy)
	{
		double x, y;
		x = dfGeoTransform[0]
		+ dx * dfGeoTransform[1] 
		+ dy * dfGeoTransform[2];
		y = dfGeoTransform[3]
		+ dx * dfGeoTransform[4] 
		+ dy * dfGeoTransform[5];

		return OGRPoint(x,y);
	}

private:
	double dfGeoTransform[6];
};

int main()
{
	const char *pszDriverName = "DXF";
	OGRSFDriver *poDriver;
	
	RegisterOGRDXF();

	CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
	CPLSetConfigOption("GDAL_DATA", "C:\\GDAL\\data");  //改了这句话就对了:OGR创建DXF格式需要data文件夹下的head.dxf,所以首先要设置GDAL_DATA的目录。
	poDriver = OGRSFDriverRegistrar::GetRegistrar()->GetDriverByName(
		pszDriverName );
	if( poDriver == NULL )
	{
		printf( "%s driver not available.\n", pszDriverName );
		exit( 1 );
	}

	OGRDataSource *poDS;
	const char* pszFileName="out.dxf";
	poDS = poDriver->CreateDataSource( pszFileName, NULL );
	if( poDS == NULL )
	{
		printf( "Creation of output file failed.\n" );
		exit( 1 );
	}

	OGRLayer *poLayer;

	poLayer = poDS->CreateLayer( "out", NULL, wkbUnknown, NULL );
	if( poLayer == NULL )
	{
		printf( "Layer creation failed.\n" );
		exit( 1 );
	}

	OGRFeature *poFeature = OGRFeature::CreateFeature( poLayer->GetLayerDefn() );

	// 原始矩形
	OGRLinearRing oSquare;
	oSquare.addPoint(2.0,2.0);
	oSquare.addPoint(4.0,2.0);
	oSquare.addPoint(4.0,3.0);
	oSquare.addPoint(2.0,3.0);
	oSquare.closeRings();
	poFeature->SetGeometry( &oSquare );
	if( poLayer->CreateFeature( poFeature ) != OGRERR_NONE )
	{
		printf( "Failed to create feature in dxffile.\n" );
		exit( 1 );
	}

	// 定义转换器
	GeoTransform geoTransform;
	// 平移
	geoTransform.SetTranslate(2.0,3.0);
	// 旋转
	geoTransform.SetRotate(2);
	// 缩放比例
	geoTransform.SetScale(2.0,3.0);

	// 转换坐标之后的矩形
	OGRLinearRing oSquare2;
	oSquare2.addPoint(&geoTransform.Transform(2.0,2.0));
	oSquare2.addPoint(&geoTransform.Transform(4.0,2.0));
	oSquare2.addPoint(&geoTransform.Transform(4.0,3.0));
	oSquare2.addPoint(&geoTransform.Transform(2.0,3.0));
	oSquare2.closeRings();
	poFeature->SetGeometry( &oSquare2 );
	if( poLayer->CreateFeature( poFeature ) != OGRERR_NONE )
	{
		printf( "Failed to create feature in dxffile.\n" );
		exit( 1 );
	}

	OGRFeature::DestroyFeature( poFeature );

	OGRDataSource::DestroyDataSource( poDS );
}

### 利用地理参数计算遥感影像四个角点的经纬度坐标 地理参数是一种常见的仿射变换模型,用于描述遥感影像与其对应的地理坐标之间的关系。这参数通常包括左上角像素的地理坐标(经度和纬度)、横向和纵向的空间分辨率以及可能存在的旋转角度。 #### 计算公式 假设给定的参数为 \(a, b, c, d, e, f\),其中: - \(a\): 左上角像素中心的横坐标偏移量; - \(b\): 表示列方向上的旋转/倾斜因子; - \(c\): 像素宽度(沿X轴的方向分辨率); - \(d\): 行方向上的旋转/倾斜因子; - \(e\): 像素高度(沿Y轴的方向分辨率),通常为负数表示向下移动时纵坐标减小; - \(f\): 左上角像素中心的纵坐标偏移量。 对于任意像素位置 \((col, row)\),其对应的地理坐标可以通过以下公式计算: \[ Longitude = a + col \cdot c + row \cdot b \] \[ Latitude = f + col \cdot e + row \cdot d \] 这里需要注意的是,\(row\) 和 \(col\) 是从零开始计数的像素索引[^1]。 #### 实现代码 以下是实现上述公式的 Python 代码片段: ```python def calculate_coordinates(a, b, c, d, e, f, width, height): """ 使用参数法计算遥感影像四角点的地理坐标。 参数: a (float): 地理坐标的初始横坐标偏移量 b (float): 列方向上的旋转/倾斜因子 c (float): 像素宽度(沿X轴) d (float): 行方向上的旋转/倾斜因子 e (float): 像素高度(沿Y轴) f (float): 地理坐标的初始纵坐标偏移量 width (int): 影像宽度(列数) height (int): 影像高度(行数) 返回: list of tuples: 四个角点的地理坐标 [(lon, lat), ...] """ corners = [ (0, 0), # 左上角 (width - 1, 0), # 右上角 (width - 1, height - 1), # 右下角 (0, height - 1) # 左下角 ] results = [] for col, row in corners: lon = a + col * c + row * b lat = f + col * e + row * d results.append((lon, lat)) return results # 示例调用 params = { 'a': 120.0, 'b': 0.0, 'c': 0.0001, 'd': 0.0, 'e': -0.0001, 'f': 30.0, } image_size = {'width': 1000, 'height': 1000} coordinates = calculate_coordinates( params['a'], params['b'], params['c'], params['d'], params['e'], params['f'], image_size['width'], image_size['height'] ) print(coordinates) ``` 此函数会返回一个列表,包含四个角点的地理坐标[(lon, lat)],分别对应于左上、右上、右下和左下的位置。 #### RPC 模型的应用扩展 虽然题目主要涉及参数模型,但在实际应用中,RPC(Rational Polynomial Coefficients)模型更为通用且精确。RPC 模型能够更灵活地适应复杂的成像条件,尤其适用于高空间分辨率卫星影像的数据处理[^5]。如果需要进一步提高精度,则可考虑引入 GDAL 库中的 RPC 支持功能来进行坐标转换
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值