计算两幅影像的几何定位精度

---

1.准备文件

（1）准备好要比较的两幅影像
（2）在envi协助下选取两幅影像的同名点，并保存文件 .pts

2.matlab几何定位精度比较程序

（1）主程序

fclose all
clear all
clc

%% 
% 计算影像的几何定位精度
% 需要的文件： envi上对选取同名点的文件。
img_left_name='F:\GF_2022_3_14_camgen\mapproject_result\ba1_pm1_012_pinhole_map.tif';
img_right_name='F:\GF_2022_3_14_camgen\mapproject_result\ba1_pm1_006_pinhole_map.tif';
f_pts='F:\GF_2022_3_14_camgen\mapproject_result\pm1_012_006_pinhole_map.pts';

[img_left_px,img_left_py,img_right_px,img_right_py]=textread(f_pts,'%s%s%s%s','headerlines',5);

[img_left_data,img_left_geoinfo] = geotiffread(img_left_name);

[img_right_data,img_right_geoinfo] = geotiffread(img_right_name);

% 大概思路： 求出两影像同名点的经纬度坐标，再将这经纬度坐标转为UTM坐标，然后进行计算精度

% img_left_geoinfo.LongitudeLimits
% img_left_geoinfo.LatitudeLimits
% img_left_geoinfo.CellExtentInLatitude
% img_left_geoinfo.CellExtentInLongitude


img_left_lon=img_left_geoinfo.LongitudeLimits(1)+img_left_geoinfo.CellExtentInLongitude*str2num(char(img_left_px));
img_left_lat=img_left_geoinfo.LatitudeLimits(2)-img_left_geoinfo.CellExtentInLatitude*str2num(char(img_left_py));

img_right_lon=img_right_geoinfo.LongitudeLimits(1)+img_right_geoinfo.CellExtentInLongitude*str2num(char(img_right_px));
img_right_lat=img_right_geoinfo.LatitudeLimits(2)-img_right_geoinfo.CellExtentInLatitude*str2num(char(img_right_py));


[img_left_utm_x,img_left_utm_y,img_left_utm_zone]=ll2utm(img_left_lat,img_left_lon,'wgs84');
[img_right_utm_x,img_right_utm_y,img_right_utm_zone]=ll2utm(img_right_lat,img_right_lon,'wgs84');

x_utm_chazhi=img_left_utm_x-img_right_utm_x;   % 同名点在UTM坐标系下x方向差值，单位是米 

y_utm_chazhi=img_left_utm_y-img_right_utm_y;   % 同名点在UTM坐标系下y方向差值，单位是米 

（2）附加的经纬度坐标转UTM坐标程序

function [x,y,f]=ll2utm(lat,lon,datum)
%LL2UTM Lat/Lon to UTM coordinates precise conversion.
%   [X,Y]=LL2UTM(LAT,LON) or LL2UTM([LAT,LON]) converts coordinates 
%	LAT,LON (in degrees) to UTM X and Y (in meters). Default datum is WGS84.
%
%	LAT and LON can be scalars, vectors or matrix. Outputs X and Y will
%	have the same size as inputs.
%
%	LL2UTM(LAT,LON,DATUM) uses specific DATUM for conversion. DATUM can be
%	a string in the following list:
%		'wgs84': World Geodetic System 1984 (default)
%		'nad27': North American Datum 1927
%		'clk66': Clarke 1866
%		'nad83': North American Datum 1983
%		'grs80': Geodetic Reference System 1980
%		'int24': International 1924 / Hayford 1909
%	or DATUM can be a 2-element vector [A,F] where A is semimajor axis (in
%	meters)	and F is flattening of the user-defined ellipsoid.
%
%	[X,Y,ZONE]=LL2UTM(...) returns also computed UTM ZONE (negative value 
%	stands for southern hemisphere points).
%
%
%	XY = LL2UTM(...) or without any output argument returns a 2-column 
%	matrix [X,Y].
%
%	Notice:
%		- LL2UTM does not perform cross-datum conversion.
%		- precision is near a millimeter.
%
%
%	Reference:
%		I.G.N., Projection cartographique Mercator Transverse: Algorithmes,
%		   Notes Techniques NT/G 76, janvier 1995.
%
%   Author: Francois Beauducel, <beauducel@ipgp.fr>
%   Created: 2003-12-02
%   Updated: 2014-04-20


%	Copyright (c) 2001-2014, Fran鏾is Beauducel, covered by BSD License.
%	All rights reserved.
%
%	Redistribution and use in source and binary forms, with or without 
%	modification, are permitted provided that the following conditions are 
%	met:
%
%	   * Redistributions of source code must retain the above copyright 
%	     notice, this list of conditions and the following disclaimer.
%	   * Redistributions in binary form must reproduce the above copyright 
%	     notice, this list of conditions and the following disclaimer in 
%	     the documentation and/or other materials provided with the distribution
%	                           
%	THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 
%	AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 
%	IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 
%	ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE 
%	LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 
%	CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 
%	SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS 
%	INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 
%	CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 
%	ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 
%	POSSIBILITY OF SUCH DAMAGE.

% Available datums
datums = [ ...
	{ 'wgs84', 6378137.0, 298.257223563 };
	{ 'nad83', 6378137.0, 298.257222101 };
	{ 'grs80', 6378137.0, 298.257222101 };
	{ 'nad27', 6378206.4, 294.978698214 };
	{ 'int24', 6378388.0, 297.000000000 };
	{ 'clk66', 6378206.4, 294.978698214 };
];

if nargin == 1
	if size(lat,2) ~= 2
		error('Single input argument must be a 2-column matrix [LAT,LON].')
	end
	lon = lat(:,2);
	lat = lat(:,1);
end
	
if nargin < 1
	error('Not enough input arguments.')
end

if all([numel(lat),numel(lon)] > 1) && any(size(lat) ~= size(lon))
	error('LAT and LON must be the same size or scalars.')
end

if nargin < 3
	datum = 'wgs84';
end

if ischar(datum)
	if ~any(strcmpi(datum,datums(:,1)))
		error('Unkown DATUM name "%s"',datum);
	end
	k = find(strcmpi(datum,datums(:,1)));
	A1 = datums{k,2};
	F1 = datums{k,3};	
else
	if numel(datum) ~= 2
		error('User defined DATUM must be a vector [A,F].');
	end
	A1 = datum(1);
	F1 = datum(2);
end

% constants
D0 = 180/pi;	% conversion rad to deg

K0 = 0.9996;	% UTM scale factor
X0 = 500000;	% UTM false East (m)

B1 = A1*(1 - 1/F1);
E1 = sqrt((A1*A1 - B1*B1)/(A1*A1));

p1 = lat/D0;					% Phi = Latitude (rad)
l1 = lon/D0;					% Lambda = Longitude (rad)
F0 = round((l1*D0 + 183)/6);	% UTM zone
P0 = 0/D0;
L0 = (6*F0 - 183)/D0;			% UTM origin longitude (rad)
Y0 = 1e7*(p1 < 0);				% UTM false northern (m)
N = K0*A1;

C = coef(E1,0);
B = C(1)*P0 + C(2)*sin(2*P0) + C(3)*sin(4*P0) + C(4)*sin(6*P0) + C(5)*sin(8*P0);
YS = Y0 - N*B;

C = coef(E1,2);
L = log(tan(pi/4 + p1/2).*(((1 - E1*sin(p1))./(1 + E1*sin(p1))).^(E1/2)));
z = complex(atan(sinh(L)./cos(l1 - L0)),log(tan(pi/4 + asin(sin(l1 - L0)./cosh(L))/2)));
Z = N.*C(1).*z + N.*(C(2)*sin(2*z) + C(3)*sin(4*z) + C(4)*sin(6*z) + C(5)*sin(8*z));
xs = imag(Z) + X0;
ys = real(Z) + YS;

% outputs zone if needed: scalar value if unique, or vector/matrix of the
% same size as x/y in case of crossed zones
if nargout > 2
	fu = unique(F0.*sign(lat));
	if isscalar(fu)
		f = fu;
	else
		f = F0;
	end
end

if nargout < 2
	x = [xs(:),ys(:)];
else
	x = xs;
	y = ys;
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function c = coef(e,m)
%COEF Projection coefficients
%	COEF(E,M) returns a vector of 5 coefficients from:
%		E = first ellipsoid excentricity
%		M = 0 for transverse mercator
%		M = 1 for transverse mercator reverse coefficients
%		M = 2 for merdian arc


if nargin < 2
	m = 0;
end

switch m
	case 0
	c0 = [-175/16384, 0,   -5/256, 0,  -3/64, 0, -1/4, 0, 1;
           -105/4096, 0, -45/1024, 0,  -3/32, 0, -3/8, 0, 0;
           525/16384, 0,  45/1024, 0, 15/256, 0,    0, 0, 0;
          -175/12288, 0, -35/3072, 0,      0, 0,    0, 0, 0;
          315/131072, 0,        0, 0,      0, 0,    0, 0, 0];
	  
	case 1
	c0 = [-175/16384, 0,   -5/256, 0,  -3/64, 0, -1/4, 0, 1;
             1/61440, 0,   7/2048, 0,   1/48, 0,  1/8, 0, 0;
          559/368640, 0,   3/1280, 0,  1/768, 0,    0, 0, 0;
          283/430080, 0, 17/30720, 0,      0, 0,    0, 0, 0;
       4397/41287680, 0,        0, 0,      0, 0,    0, 0, 0];

	case 2
	c0 = [-175/16384, 0,   -5/256, 0,  -3/64, 0, -1/4, 0, 1;
         -901/184320, 0,  -9/1024, 0,  -1/96, 0,  1/8, 0, 0;
         -311/737280, 0,  17/5120, 0, 13/768, 0,    0, 0, 0;
          899/430080, 0, 61/15360, 0,      0, 0,    0, 0, 0;
      49561/41287680, 0,        0, 0,      0, 0,    0, 0, 0];
   
end
c = zeros(size(c0,1),1);

for i = 1:size(c0,1)
    c(i) = polyval(c0(i,:),e);
end