UTM坐标中的大圆

Question

我有一个关于笛卡尔坐标系下大圆计算的问题

是否有可能在UTM坐标中计算沿大圆的点，即给定两个点P1(easting1, northing1)和P2 (easting2, northing2)，是否有一个公式来计算沿连接这两个点的测地线上的路点(不转换为地理坐标并返回)?

任何帮助都是感激的!谢谢!

Answer 1

UTM坐标系是一种横向墨卡托投影，以经度带分隔，纬度范围限制，因此与投影相关的失真保持较小。同样，UTM坐标系是正形投影．因此，它保留了角度。这意味着在UTM区域内，UTM投影中的直线是在UTM区域内相对较短距离上的一个很好的近似大圆。如果你不需要很高的精度，你可以只画出P1和P2之间的直线。

下面的图表显示了直线与大圆的最大偏差。在本例中，P1为赤道和经度3°E上的红点，并计算显示框中所有6°x6°的点的最大偏差。

虽然在几百公里的距离内差异很小，但在更大的距离内，差异会更大，特别是当P1和P2位于不同的UTM区域时。

如果应用程序需要更大的距离，就不应该使用UTM投影。或者，如果您的应用程序要求在适合UTM投影的距离上有更高的精度，您可以计算大圆轨迹，但我不知道任何直接公式。

我所使用的，我知道这不是你想要的，是一些Matlab函数，用于转换UTM坐标和另一个计算地理坐标中的大圆轨道:

函数[lat,lon] = greatCircle(P1,P2,N) % greatCircle计算大圆圆弧沿点连接点P1 %和P2 % P1,P2格式:[lat,lon] % N:大圆圆弧沿点数lat1=P1(1);lon1 = P1 (2);lat2 = P2 (1);lon2 = P2 (2);lon12 = lon2-lon1;a1 = atan2d (cosd (lat2)。*信德(lon12), (cosd (lat1) *信德(lat2)信德(lat1) * cosd (lat1)。* cosd (lon12)));% a2 = atan2d (cosd (lat1)。*信德(lon12), (-cosd (lat2) *信德(lat1) +信德(lat2) * cosd (lat1)。* cosd (lon12)));a0 = atand(信德(a1)。* cosd (lat1)。/√((cosd (a1)。^ 2)+(信德(a1)。^ 2)。*(信德(lat1)。^ 2)));sigma01 = atan2d(罐内(lat1) cosd (a1));sigma12 = atan2d(√(cosd (lat1) *信德(lat2)信德(lat1) * cosd (lat2) * cosd (lon12)) ^ 2 + (cosd (lat2) *信德(lon12)) ^ 2),(信德(lat1) *信德(lat2) + cosd (lat1) * cosd (lat2) * cosd (lon12))); sigma02=sigma01+sigma12; lon01=atan2d(sind(a0)*sind(sigma01),cosd(sigma01)); lon0=lon1-lon01; sigma=linspace(sigma01,sigma02,N); lat=atan2d(cosd(a0)*sind(sigma),sqrt(cosd(sigma).^2+sind(a0)^2*sind(sigma).^2)); lon=wrapTo360(atan2d(sind(a0)*sind(sigma),cosd(sigma))+lon0); end

对于转换:

函数(x, y, utmzone) = deg2utm(纬度,经度,黑龙江 ) % ------------------------------------------------------------------------- % [ x, y, utmzone] = deg2utm(纬度、经度)% %描述:函数将纬度/经度向量转换为UTM坐标(WGS84)。已经从UTM中提取了一些代码。m函数由Gabriel Ruiz Martinez设计。% %输入:% Lat:纬度向量。度。+ ddd。Lon:经度向量。度。+ ddd。ddddd WGS84 % %输出:% x, y, utmzone。例1:% Lat=[40.3154333; 46.283900; 37.577833; 28.645650; 38.855550; 25.061783]; % Lon=[-3.4857166; 7.8012333; -119.95525; -17.759533; -94.7990166; 121.640266]; % [x,y,utmzone] = deg2utm(Lat,Lon); % fprintf('%7.0f ',x) % 458731 407653 239027 230253 343898 362850 % fprintf('%7.0f ',y) % 4462881 5126290 4163083 3171843 4302285 2772478 % utmzone = % 30 T % 32 T % 11 S % 28 R % 15 S % 51 R % % Example 2: If you have Lat/Lon coordinates in Degrees, Minutes and Seconds % LatDMS=[40 18 55.56; 46 17 2.04]; % LonDMS=[-3 29 8.58; 7 48 4.44]; % Lat=dms2deg(mat2dms(LatDMS)); %convert into degrees % Lon=dms2deg(mat2dms(LonDMS)); %convert into degrees % [x,y,utmzone] = deg2utm(Lat,Lon) % % Author: % Rafael Palacios % Universidad Pontificia Comillas % Madrid, Spain % Version: Apr/06, Jun/06, Aug/06, Aug/06 % Aug/06: fixed a problem (found by Rodolphe Dewarrat) related to southern % hemisphere coordinates. % Aug/06: corrected m-Lint warnings %------------------------------------------------------------------------- % Argument checking % error(nargchk(2, 3, nargin)); %2 arguments required CHuso=true; if nargin==3 CHuso=false; end n1=length(Lat); n2=length(Lon); if (n1~=n2) error('Lat and Lon vectors should have the same length'); end % Memory pre-allocation % x=zeros(n1,1); y=zeros(n1,1); utmzone(n1,:)='60 X'; % Main Loop % for i=1:n1 la=Lat(i); lo=Lon(i); sa = 6378137.000000 ; sb = 6356752.314245; %e = ( ( ( sa ^ 2 ) - ( sb ^ 2 ) ) ^ 0.5 ) / sa; e2 = ( ( ( sa ^ 2 ) - ( sb ^ 2 ) ) ^ 0.5 ) / sb; e2cuadrada = e2 ^ 2; c = ( sa ^ 2 ) / sb; %alpha = ( sa - sb ) / sa; %f %ablandamiento = 1 / alpha; % 1/f lat = la * ( pi / 180 ); lon = lo * ( pi / 180 ); if CHuso Huso = fix( ( lo / 6 ) + 31); end S = ( ( Huso * 6 ) - 183 ); deltaS = lon - ( S * ( pi / 180 ) ); if (la<-72), Letra='C'; elseif (la<-64), Letra='D'; elseif (la<-56), Letra='E'; elseif (la<-48), Letra='F'; elseif (la<-40), Letra='G'; elseif (la<-32), Letra='H'; elseif (la<-24), Letra='J'; elseif (la<-16), Letra='K'; elseif (la<-8), Letra='L'; elseif (la<0), Letra='M'; elseif (la<8), Letra='N'; elseif (la<16), Letra='P'; elseif (la<24), Letra='Q'; elseif (la<32), Letra='R'; elseif (la<40), Letra='S'; elseif (la<48), Letra='T'; elseif (la<56), Letra='U'; elseif (la<64), Letra='V'; elseif (la<72), Letra='W'; else Letra='X'; end a = cos(lat) * sin(deltaS); epsilon = 0.5 * log( ( 1 + a) / ( 1 - a ) ); nu = atan( tan(lat) / cos(deltaS) ) - lat; v = ( c / ( ( 1 + ( e2cuadrada * ( cos(lat) ) ^ 2 ) ) ) ^ 0.5 ) * 0.9996; ta = ( e2cuadrada / 2 ) * epsilon ^ 2 * ( cos(lat) ) ^ 2; a1 = sin( 2 * lat ); a2 = a1 * ( cos(lat) ) ^ 2; j2 = lat + ( a1 / 2 ); j4 = ( ( 3 * j2 ) + a2 ) / 4; j6 = ( ( 5 * j4 ) + ( a2 * ( cos(lat) ) ^ 2) ) / 3; alfa = ( 3 / 4 ) * e2cuadrada; beta = ( 5 / 3 ) * alfa ^ 2; gama = ( 35 / 27 ) * alfa ^ 3; Bm = 0.9996 * c * ( lat - alfa * j2 + beta * j4 - gama * j6 ); xx = epsilon * v * ( 1 + ( ta / 3 ) ) + 500000; yy = nu * v * ( 1 + ta ) + Bm; if (yy<0) yy=9999999+yy; end x(i)=xx; y(i)=yy; utmzone(i,:)=sprintf('%02d %c',Huso,Letra); end

和

函数(纬度、经度)= utm2deg (yy, xx utmzone ) % ------------------------------------------------------------------------- % [ 纬度,经度]= utm2deg (x, y, utmzone) % %描述:函数向量的UTM坐标转化为纬度/经度向量(WGS84)。已经从UTMIP中提取了一些代码。m函数由Gabriel Ruiz Martinez设计。% %输入:% x, y, utmzone。% %输出:% Lat:纬度向量。度。+ ddd。Lon:经度向量。度。+ ddd。例1:% x=[458731; 407653; 239027; 230253; 343898; 362850]; % y=[4462881; 5126290; 4163083; 3171843; 4302285; 2772478]; % utmzone=['30 T'; '32 T'; '11 S'; '28 R'; '15 S'; '51 R']; % [Lat, Lon]=utm2deg(x,y,utmzone); % fprintf('%11.6f ',lat) % 40.315430 46.283902 37.577834 28.645647 38.855552 25.061780 % fprintf('%11.6f ',lon) % -3.485713 7.801235 -119.955246 -17.759537 -94.799019 121.640266 % % Example 2: If you need Lat/Lon coordinates in Degrees, Minutes and Seconds % [Lat, Lon]=utm2deg(x,y,utmzone); % LatDMS=dms2mat(deg2dms(Lat)) %LatDMS = % 40.00 18.00 55.55 % 46.00 17.00 2.01 % 37.00 34.00 40.17 % 28.00 38.00 44.33 % 38.00 51.00 19.96 % 25.00 3.00 42.41 % LonDMS=dms2mat(deg2dms(Lon)) %LonDMS = % -3.00 29.00 8.61 % 7.00 48.00 4.40 % -119.00 57.00 18.93 % -17.00 45.00 34.33 % -94.00 47.00 56.47 % 121.00 38.00 24.96 % % Author: % Rafael Palacios % Universidad Pontificia Comillas % Madrid, Spain % Version: Apr/06, Jun/06, Aug/06 % Aug/06: corrected m-Lint warnings %------------------------------------------------------------------------- % Argument checking % error(nargchk(3, 3, nargin)); %3 arguments required n1=length(xx); n2=length(yy); n3=size(utmzone,1); if (n1~=n2 || n1~=n3) error('x,y and utmzone vectors should have the same number or rows'); end c=size(utmzone,2); if (c~=4) error('utmzone should be a vector of strings like "30 T"'); end % Memory pre-allocation % Lat=zeros(n1,1); Lon=zeros(n1,1); % Main Loop % for i=1:n1 if (utmzone(i,4)>'X' || utmzone(i,4)<'C') fprintf('utm2deg: Warning utmzone should be a vector of strings like "30 T", not "30 t"\n'); end if (utmzone(i,4)>'M') hemis='N'; % Northern hemisphere else hemis='S'; end x=xx(i); y=yy(i); zone=str2double(utmzone(i,1:2)); sa = 6378137.000000 ; sb = 6356752.314245; % e = ( ( ( sa ^ 2 ) - ( sb ^ 2 ) ) ^ 0.5 ) / sa; e2 = ( ( ( sa ^ 2 ) - ( sb ^ 2 ) ) ^ 0.5 ) / sb; e2cuadrada = e2 ^ 2; c = ( sa ^ 2 ) / sb; % alpha = ( sa - sb ) / sa; %f % ablandamiento = 1 / alpha; % 1/f X = x - 500000; if hemis == 'S' || hemis == 's' Y = y - 10000000; else Y = y; end S = ( ( zone * 6 ) - 183 ); lat = Y / ( 6366197.724 * 0.9996 ); v = ( c / ( ( 1 + ( e2cuadrada * ( cos(lat) ) ^ 2 ) ) ) ^ 0.5 ) * 0.9996; a = X / v; a1 = sin( 2 * lat ); a2 = a1 * ( cos(lat) ) ^ 2; j2 = lat + ( a1 / 2 ); j4 = ( ( 3 * j2 ) + a2 ) / 4; j6 = ( ( 5 * j4 ) + ( a2 * ( cos(lat) ) ^ 2) ) / 3; alfa = ( 3 / 4 ) * e2cuadrada; beta = ( 5 / 3 ) * alfa ^ 2; gama = ( 35 / 27 ) * alfa ^ 3; Bm = 0.9996 * c * ( lat - alfa * j2 + beta * j4 - gama * j6 ); b = ( Y - Bm ) / v; Epsi = ( ( e2cuadrada * a^ 2 ) / 2 ) * ( cos(lat) )^ 2; Eps = a * ( 1 - ( Epsi / 3 ) ); nab = ( b * ( 1 - Epsi ) ) + lat; senoheps = ( exp(Eps) - exp(-Eps) ) / 2; Delt = atan(senoheps / (cos(nab) ) ); TaO = atan(cos(Delt) * tan(nab)); longitude = (Delt *(180 / pi ) ) + S; latitude = ( lat + ( 1 + e2cuadrada* (cos(lat)^ 2) - ( 3 / 2 ) * e2cuadrada * sin(lat) * cos(lat) * ( TaO - lat ) ) * ( TaO - lat ) ) * ... (180 / pi); Lat(i)=latitude; Lon(i)=longitude; end

我希望这对你有所帮助。